# Influence of low orbit design and strategies for gravity field recovery of Europa

William Desprats, Stefano Bertone, Daniel Arnold, Adrian Jäggi, Michel Blanc

► To cite this version:

William Desprats, Stefano Bertone, Daniel Arnold, Adrian Jäggi, Michel Blanc. Influence of low orbit design and strategies for gravity field recovery of Europa. Planetary and Space Science, 2023, 226, 10.1016/j.pss.2022.105631. insu-04473219

HAL Id: insu-04473219

https://insu.hal.science/insu-04473219

Submitted on 22 Feb 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_1.png" /></figure>

# Influence of low orbit design and strategies for gravity field recovery of Europa

<figure><img src="image_2.png" /></figure>

William Desprats<sup>a, *</sup>, Stefano Bertone<sup>b, c</sup>, Daniel Arnold<sup>a</sup>, Adrian Jäggı<sup>a</sup>, Michel Blanc<sup>d</sup>

## ARTICLE INFO

Handling Editor: Dr Olivier Witasse

**Keywords:**

Europa

Orbit determination

Gravity field

Simulation study

## ABSTRACT

The gravity field of the Solar System celestial bodies is generally estimated within the orbit determination process of probes visiting them. Low altitude, high inclination and near circular orbits are best suited for this purpose. As part of a global characterization effort of the Jovian moon Europa, we studied the influence of multiple orbital parameters and configurations on the recovery of its gravity field parameters. We made use of Repetitive Ground Track Orbits (RGTO), allowing for a proper definition of the ground track coverage on the surface of Europa, which is essential to study the impact of this coverage. The results presented here rely on closed-loop simulations performed using a development version of the Bernese GNSS Software with planetary geodesy capabilities. We simulate realistic range-rate observables (2-way Doppler X-band, $\sigma = 0.10$ mm/s) from the different orbits considered over a total mission duration of 3 months. These observations are then used to reconstruct the orbit and estimate gravity field coefficients in terms of spherical harmonic coefficients and the $k_2$ Love number. We systematically compare solutions based on different input orbital parameters and we quantify their impact on the gravity field recovery process, which is of great importance for future mission designs. Our best case scenario shows that the gravity field can be estimated up to degree and order 72 after 3 months in a circular polar orbit at 100 km over Europa's surface. Different gravity field recovery strategies are also discussed when starting with a very poor a priori knowledge of the gravity field, as it is the case for the Galileans moons. We propose and evaluate two approaches by either using pseudo-stochastic pulses (i.e., instantaneous velocity changes) to cope with the large model deficiencies, or by co-estimating low-degree gravity field coefficients and orbit parameters to bootstrap the estimation process.

## 1. Introduction

Europa is a privileged destination for the upcoming phase of Solar System exploration, as clues for its habitability have been accumulating during the past two decades. Data analysis of the Galileo mission suggested that Europa's internal structure is differentiated (Anderson, 1998) into a metallic core, a silicate mantle and a water ice-liquid outer shell. A subsurface ocean of liquid water beneath the icy crust was subsequently inferred from magnetometer data of the Galileo mission (Kivelson et al., 2000). Furthermore, the geological activity of Europa was brought to light with evidence that water plumes emerge from the surface (Roth et al., 2014).

As part of a global characterization of Europa, the study of its gravity

field and tides is essential to better constrain internal structure models. Dedicated missions to Europa will help to acquire a better understanding of this complex system, and of various interactions of Europa within the Jovian system (Bignami et al., 2005; National Research Council, 2011).

In this sense, ESA's mission JUpiter ICy Moons Explorer (JUICE: Grasset et al., 2013) will perform a series of fly-bys of the three outer Galilean moons starting in 2030. However, Ganymede is the main focus of JUICE and only 2 flybys of Europa are planned. NASA's Europa Clipper is expected to be launched by 2024 and will perform a series of 46 flybys of Europa starting in 2030, which is expected to largely contribute to our knowledge of this moon (Phillips and Pappalardo, 2014). Several studies already investigated the recovery of the low-degrees gravity field coefficients of Europa from flyby missions (Park et al., 2015; Mazarico

<footnote><sup>a</sup> Astronomical Institute, University of Bern, Switzerland</footnote>

<footnote><sup>b</sup> Center for Research and Exploration in Space Science and Technology, USA</footnote>

<footnote><sup>c</sup> INAF, Astrophysical Observatory of Torino, Italy</footnote>

<footnote><sup>d</sup> IRAP, CNRS-Université Paul Sabatier, Toulouse, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: william.desprats@unibe.ch (W. Desprats).</footnote>
<figure><img src="image_2.png" /></figure>

# Influence of low orbit design and strategies for gravity field recovery of Europa

<figure><img src="image_3.png" /></figure>

Check for
updates

William Desprats<sup>a, *,</sup> , Stefano Bertone<sup>b, c</sup> , Daniel Arnold<sup>a</sup> , Adrian Jäggı<sup>a</sup> , Michel Blanc<sup>d</sup>

ARTICLE INFO

Handling Editor: Dr Olivier Witasse

Keywords:

Europa

Orbit determination

Gravity field

Simulation study

A B S T R A C T

The gravity field of the Solar System celestial bodies is generally estimated within the orbit determination process of probes visiting them. Low altitude, high inclination and near circular orbits are best suited for this purpose. As part of a global characterization effort of the Jovian moon Europa, we studied the influence of multiple orbital parameters and configurations on the recovery of its gravity field parameters. We made use of Repetitive Ground Track Orbits (RGTO), allowing for a proper definition of the ground track coverage on the surface of Europa, which is essential to study the impact of this coverage. The results presented here rely on closed-loop simulations performed using a development version of the Bernese GNSS Software with planetary geodesy capabilities. We simulate realistic range-rate observables (2-way Doppler X-band, σ = 0.10 mm/s) from the different orbits considered over a total mission duration of 3 months. These observations are then used to reconstruct the orbit and estimate gravity field coefficients in terms of spherical harmonic coefficients and the k₂ Love number. We systematically compare solutions based on different input orbital parameters and we quantify their impact on the gravity field recovery process, which is of great importance for future mission designs. Our best case scenario shows that the gravity field can be estimated up to degree and order 72 after 3 months in a circular polar orbit at 100 km over Europa's surface. Different gravity field recovery strategies are also discussed when starting with a very poor a priori knowledge of the gravity field, as it is the case for the Galileans moons. We propose and evaluate two approaches by either using pseudo-stochastic pulses (i.e., instantaneous velocity changes) to cope with the large model deficiencies, or by co-estimating low-degree gravity field coefficients and orbit parameters to bootstrap the estimation process.

1. Introduction

Europa is a privileged destination for the upcoming phase of Solar System exploration, as clues for its habitability have been accumulating during the past two decades. Data analysis of the Galileo mission suggested that Europa's internal structure is differentiated (Anderson, 1998) into a metallic core, a silicate mantle and a water ice-liquid outer shell. A subsurface ocean of liquid water beneath the icy crust was subsequently inferred from magnetometer data of the Galileo mission (Kivelson et al., 2000). Furthermore, the geological activity of Europa was brought to light with evidence that water plumes emerge from the surface (Roth et al., 2014).

As part of a global characterization of Europa, the study of its gravity

field and tides is essential to better constrain internal structure models. Dedicated missions to Europa will help to acquire a better understanding of this complex system, and of various interactions of Europa within the Jovian system (Bignami et al., 2005; National Research Council, 2011).

In this sense, ESA's mission JUpiter ICy Moons Explorer (JUICE: Grasset et al., 2013) will perform a series of fly-bys of the three outer Galilean moons starting in 2030. However, Ganymede is the main focus of JUICE and only 2 flybys of Europa are planned. NASA's Europa Clipper is expected to be launched by 2024 and will perform a series of 46 flybys of Europa starting in 2030, which is expected to largely contribute to our knowledge of this moon (Phillips and Pappalardo, 2014). Several studies already investigated the recovery of the low-degrees gravity field co- efficients of Europa from flyby missions (Park et al., 2015; Mazarico

<footnote><sup>a</sup> Astronomical Institute, University of Bern, Switzerland</footnote>

<footnote>b Center for Research and Exploration in Space Science and Technology, USA</footnote>

<footnote>c INAF, Astrophysical Observatory of Torino, Italy</footnote>

<footnote>d IRAP, CNRS-Université Paul Sabatier, Toulouse, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: william.desprats@unibe.ch (W. Desprats).</footnote>
et al., 2015; Verma and Margot, 2018).

However, a mission including an orbiter around Europa would enable global and uniform mapping of the Galilean moon, as well as recovering its gravity field and other geodetic parameters to much higher resolution (Wu et al., 2001; Wahr et al., 2006). An orbiter would thus achieve a more detailed characterization of the presence and the extent of the ocean and its relation to the deeper interior. There have been several proposed missions to send an orbiter around Europa, e.g., the Jupiter Europa Orbiter (Clark et al., 2009), HADES (Böttcher et al., 2009), and the Joint Europa Mission (Blanc et al., 2020). As it will be done for JUICE in the case of Ganymede (Cappuccio et al., 2020), a low altitude, near polar and quasi circular orbit would be in general most beneficial to improve our knowledge of Europa's gravity field by means of spacecraft Doppler tracking data. However, candidate orbits have to fit all mission science objectives and to satisfy multiple mission and instrument constraints. An orbit meeting all these additional requirements might be less suited for gravity field recovery. Quantifying the accuracy to which geodetic parameters can be recovered based on different orbits is then of great importance for a proper orbit selection.

The purpose if this study is to systematically characterize an extensive set of orbits around Europa with respect to their sensitivity to measuring the gravity field and the $k_2$ tidal Love number from the analysis of spacecraft Doppler tracking data. To that end, we perform closed-loop simulations in a development version of the Bernese (GNSS) Software (BSW; Dach et al., 2015) adapted for planetary probes (Arnold et al., 2015; Bertone et al., 2021). In particular, we produce synthetic range-rate (Doppler) measurements from a set of reference orbits. The data is then used for the determination of spacecraft orbit and geodetic parameters that are compared to the ground truth to quantify the sensitivity of different Europa-bound orbits for geodetic parameters recovery. This study follows the template of the Joint Europa Mission proposal (Blanc et al., 2020) (e.g., that the spacecraft will remain in orbit around Europa during 3 months due to the harsh environmental conditions) but our conclusions would apply to the planning of any future mission including an orbiter around Europa.

Section 2 defines the reference orbits and the different parameters considered in this study. The force model, observables adopted in our simulations, and the estimation process are described in Sec. 3. In Sec. 4, we explore several orbital configurations and their impact on gravity field determination by assuming a perfect knowledge of the a priori force model. Finally, in Sec. 5, we propose several gravity field recovery strategies, when starting from the current limited knowledge (up to degree-2 only) of Europa's gravity field.

## 2. Orbit design

The orbits considered in this study are designed in the Europa-Jupiter system. We derive a set of orbits according to parameters of interest associated with different orientations of the orbital plane with respect to other celestial bodies of interest, namely the Earth. In general, we consider low altitude, high inclination and near circular orbits to provide a global mapping of Europa to the extent possible.

### 2.1. Repetitive Ground Track Orbits

Repetitive Ground Track Orbits (RGTO) are beneficial for the observation of time varying phenomena on the ground, as repeated observations of a given point of the surface of the celestial body are ensured. RGTO are commonly used for Earth observation missions and are, e.g., also considered for JUICE around Ganymede (Ortore et al., 2015; Boutonnet and Varga, 2020).

The ground tracks of an m:R RGTO repeat after m European days (3.55 Earth days). Within this period the probe completes R revolutions around Europa. This translates into

$$mD_n = RT_n, \qquad (1)$$

where $D_n$ is the nodal day of Europa, and $T_n$ is the nodal period of the probe. m and R are two integers prime to each other, i.e., every m Europa days ground tracks return to the same position. The whole set of associated orbit tracks defines a grid in the Europa-fixed reference frame (see Fig. 1).

We limit the investigation to low altitude (100 – 200 km) and high inclination (80°–100°) orbits, because of their relevance for gravity recovery and thus for mission proposals such as Blanc et al. (2020). Since at first order, $T_n$ is given by

$$T_n = 2\pi\sqrt{a^3/\mu_E}, \qquad (2)$$

with a being the semi-major axis of the probe's orbit and $\mu_E$ the standard gravitational parameter of Europa, the nodal period $T_n$ is also bounded. For a considered revisit time m, there is a finite number of integers R satisfying the ground track repetition condition given by Eq. (1) (see Table 1).

These orbits provide a homogeneous ground track coverage. We can define the spatial resolution by the equatorial distance between the ground tracks. We call this constant distance “cycle intertrack”, and it is defined by

$$\delta_{eq} = \frac{2\pi R_E}{R}, \qquad (3)$$

with $R_E$ being the equatorial radius of Europa (1562.6 km, Archinal et al., 2018). This gap depends only on the number of orbit revolutions between repetitions, noted as R. Large cycle intertracks, i.e., a lower spatial resolution of the ground tracks, can be detrimental to the gravity field recovery. Table 1 shows the cycle intertrack for different RGTO.

It is important to note that we can still study non-repetitive orbits in the framework of RGTO if we set the repetition cycle duration to be larger or equal to the mission duration. We analyzed such RGTO with m = 26, meaning that ground tracks would only repeat after 26 Europa days (≈ 3 months), equal to our target mission duration.

Using the polynomial approach detailed in Cinelli et al. (2015), we compute a set of initial guess orbital elements (semi-major axis a, eccentricity e, inclination i and argument of periapsis ω), describing a given m:R RGTO. This approach is based on the ground track repetition condition given by Eq. (1). The semi-major axis a then satisfies

$$d_T a^7 + d_1 a^{5,5} + d_K a^4 + d_2 a^2 + d_4 = 0, \qquad (4)$$

where $d_T$ is related to third body perturbations, $d_1$ is related to orbit

<figure><img src="image_3.png" /><figcaption>Fig. 1. 3:118 RGTO (i = 89°, h = 135 km) in Europa-fixed frame. After 118 revolutions around Europa and 3 Europa days (10.65 Earth days), the probe will follow the same track.</figcaption></figure>
et al., 2015; Verma and Margot, 2018).

However, a mission including an orbiter around Europa would enable global and uniform mapping of the Galilean moon, as well as recovering its gravity field and other geodetic parameters to much higher resolution (Wu et al., 2001; Wahr et al., 2006). An orbiter would thus achieve a more detailed characterization of the presence and the extent of the ocean and its relation to the deeper interior. There have been several proposed missions to send an orbiter around Europa, e.g., the Jupiter Europa Orbiter (Clark et al., 2009), HADES (Böttcher et al., 2009), and the Joint Europa Mission (Blanc et al., 2020). As it will be done for JUICE in the case of Ganymede (Cappuccio et al., 2020), a low altitude, near polar and quasi circular orbit would be in general most beneficial to improve our knowledge of Europa's gravity field by means of spacecraft Doppler tracking data. However, candidate orbits have to fit all mission science objectives and to satisfy multiple mission and instrument constraints. An orbit meeting all these additional requirements might be less suited for gravity field recovery. Quantifying the accuracy to which geodetic parameters can be recovered based on different orbits is then of great importance for a proper orbit selection.

The purpose if this study is to systematically characterize an extensive set of orbits around Europa with respect to their sensitivity to measuring the gravity field and the $k_2$ tidal Love number from the analysis of spacecraft Doppler tracking data. To that end, we perform closed-loop simulations in a development version of the Bernese (GNSS) Software (BSW; Dach et al., 2015) adapted for planetary probes (Arnold et al., 2015; Bertone et al., 2021). In particular, we produce synthetic range-rate (Doppler) measurements from a set of reference orbits. The data is then used for the determination of spacecraft orbit and geodetic parameters that are compared to the ground truth to quantify the sensitivity of different Europa-bound orbits for geodetic parameters recovery. This study follows the template of the Joint Europa Mission proposal (Blanc et al., 2020) (e.g., that the spacecraft will remain in orbit around Europa during 3 months due to the harsh environmental conditions) but our conclusions would apply to the planning of any future mission including an orbiter around Europa.

Section 2 defines the reference orbits and the different parameters considered in this study. The force model, observables adopted in our simulations, and the estimation process are described in Sec. 3. In Sec. 4, we explore several orbital configurations and their impact on gravity field determination by assuming a perfect knowledge of the a priori force model. Finally, in Sec. 5, we propose several gravity field recovery strategies, when starting from the current limited knowledge (up to degree-2 only) of Europa's gravity field.

## 2. Orbit design

The orbits considered in this study are designed in the Europa-Jupiter system. We derive a set of orbits according to parameters of interest associated with different orientations of the orbital plane with respect to other celestial bodies of interest, namely the Earth. In general, we consider low altitude, high inclination and near circular orbits to provide a global mapping of Europa to the extent possible.

### 2.1. Repetitive Ground Track Orbits

Repetitive Ground Track Orbits (RGTO) are beneficial for the observation of time varying phenomena on the ground, as repeated observations of a given point of the surface of the celestial body are ensured. RGTO are commonly used for Earth observation missions and are, e.g., also considered for JUICE around Ganymede (Ortore et al., 2015; Boutonnet and Varga, 2020).

The ground tracks of an m:R RGTO repeat after m European days (3.55 Earth days). Within this period the probe completes R revolutions around Europa. This translates into

$$mD_n = RT_n, \qquad (1)$$

where $D_n$ is the nodal day of Europa, and $T_n$ is the nodal period of the probe. m and R are two integers prime to each other, i.e., every m Europa days ground tracks return to the same position. The whole set of associated orbit tracks defines a grid in the Europa-fixed reference frame (see Fig. 1).

We limit the investigation to low altitude (100 – 200 km) and high inclination (80°–100°) orbits, because of their relevance for gravity recovery and thus for mission proposals such as Blanc et al. (2020). Since at first order, $T_n$ is given by

$$T_n = 2\pi\sqrt{a^3/\mu_E}, \qquad (2)$$

with a being the semi-major axis of the probe's orbit and $\mu_E$ the standard gravitational parameter of Europa, the nodal period $T_n$ is also bounded. For a considered revisit time m, there is a finite number of integers R satisfying the ground track repetition condition given by Eq. (1) (see Table 1).

These orbits provide a homogeneous ground track coverage. We can define the spatial resolution by the equatorial distance between the ground tracks. We call this constant distance “cycle intertrack”, and it is defined by

$$\delta_{eq} = \frac{2\pi R_E}{R}, \qquad (3)$$

with $R_E$ being the equatorial radius of Europa (1562.6 km, Archinal et al., 2018). This gap depends only on the number of orbit revolutions between repetitions, noted as R. Large cycle intertracks, i.e., a lower spatial resolution of the ground tracks, can be detrimental to the gravity field recovery. Table 1 shows the cycle intertrack for different RGTO.

It is important to note that we can still study non-repetitive orbits in the framework of RGTO if we set the repetition cycle duration to be larger or equal to the mission duration. We analyzed such RGTO with m = 26, meaning that ground tracks would only repeat after 26 Europa days (≈ 3 months), equal to our target mission duration.

Using the polynomial approach detailed in Cinelli et al. (2015), we compute a set of initial guess orbital elements (semi-major axis a, eccentricity e, inclination i and argument of periapsis ω), describing a given m:R RGTO. This approach is based on the ground track repetition condition given by Eq. (1). The semi-major axis a then satisfies

$$d_T a^7 + d_1 a^{5.5} + d_K a^4 + d_2 a^2 + d_4 = 0, \qquad (4)$$

where $d_T$ is related to third body perturbations, $d_1$ is related to orbit

<figure><img src="image_4.png" /><figcaption>Fig. 1. 3:118 RGTO (i = 89°, h = 135 km) in Europa-fixed frame. After 118 revolutions around Europa and 3 Europa days (10.65 Earth days), the probe will follow the same track.</figcaption></figure>
<caption>**Table 1**</caption>

m:R RGTO with their altitude and cycle intertrack fulfilling the orbit requirements. The two numbers in the 2nd column denote the minimum and maximum possible values of R. A lower R corresponds to a higher altitude, and a larger cycle intertrack (see Eq. (3)). m = 1 corresponds to 3.55 Earth days duration, and m = 26 corresponds to 92.33 Earth days cycle duration.

<table><thead><tr><th>m</th><th>R</th><th>Altitude</th><th>Cycle intertrack (at equator)</th></tr></thead><tbody><tr><td rowspan="2">1</td><td>38</td><td>181 km</td><td>258 km (9.5°)</td></tr><tr><td>40</td><td>123 km</td><td>245 km (9.0°)</td></tr><tr><td rowspan="2">2</td><td>75</td><td>197 km</td><td>131 km (4.8°)</td></tr><tr><td>81</td><td>109 km</td><td>121 km (4.4°)</td></tr><tr><td rowspan="2">3</td><td>113</td><td>192 km</td><td>87 km (3.2°)</td></tr><tr><td>122</td><td>104 km</td><td>80 km (3.0°)</td></tr><tr><td rowspan="2">26</td><td>973</td><td>199 km</td><td>10 km (0.37°)</td></tr><tr><td>1061</td><td>101 km</td><td>9 km (0.34°)</td></tr></tbody></table>

characteristics, $d_K$ is related to the Keplerian motion, $d_2$ is related to $J_2$ effect and $d_4$ is related to $J_2^2$ effect. Equation (4) (which is Eq. (10) in Cinelli et al., 2015) is obtained by fixing the other orbital elements (e, i, ω), and is a function of the ratio $R/m$, of the physical parameters of Europa, and of the mean motion of Europa around Jupiter.

For each pair (m, R), we solve the degree-14 polynomial equation obtained by squaring Eq. (4) for the semi-major axis *a* while fixing *e*, *i* and *ω*. From the multiple solutions in our case only one falls within the range of studied altitudes. In this study, we do not consider elliptic orbits. We thus fix the *e* and *ω* to 0.

Equation (4) was derived by approximating third body perturbations, so that it only provides a first guess orbit. Such orbit then needs to be refined according to the full force model to ensure ground track repetition. In addition, the initial longitude of the ascending node, $\Omega$, has yet to be defined at this stage, which completes the characterization of the geometry of the orbital plane with respect to the third body Jupiter.

## 2.2. Orbit refinement

In order to precisely meet the condition for ground track repetition, we considered the Hill model (Lara and Russell, 2007) to refine our first guess orbits. This model takes into account the influence of Europa and Jupiter as point masses, plus the effect of the $J_2$ and $C_{22}$ gravity field coefficients of Europa.

We define the Jupiter-Europa rotating frame with the x-axis pointing from Jupiter to Europa, the y-axis in direction of the velocity vector of Europa around Jupiter, and the z-axis perpendicular to the orbital plane of Europa around Jupiter (see Fig. 2). In this frame centered around Europa, the spacecraft equations of motion are

$$ \ddot{x} + 2n_J \times \dot{x} = \nabla_x U_{tot}, \quad (5) $$

with

<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Orbital plane of the probe (in red) at $t_0$ and $t_s$ with respect to Jupiter and Earth. The red dots represent crossing points with respect to Europa equatorial plane ($z_0 = z_s = 0$).</figcaption></figure>

$$ U_{tot} = \frac{n_J^2}{2} (x^2 + y^2) + U_E + U_J, \quad (6) $$

and where $x = (x, y, z)$ is the spacecraft position in the rotating frame, a dot denotes time derivatives, a bold variable denotes a vector, $n_J$ is the angular velocity vector of Europa around Jupiter, $n_J$ its magnitude, and $U_E$ and $U_J$ are the gravitational potentials of Europa and Jupiter, respectively. The gravitational potential of Jupiter is then defined by

$$ U_J = \frac{n_J^2}{2} (3x^2 - r^2), \quad (7) $$

where $r = \|x\|$. Assuming hydrostatic equilibrium for Europa ($J_2 = \frac{10}{3}C_{22}$), its gravitational potential can be separated into two contributions: the potential due to central point mass and the degree-2 gravity field coefficients potential $U_p$ defined so that

$$ U_E = \frac{\mu_E}{r} + U_p = \frac{\mu_E}{r} + \frac{\mu_E R_E^2}{r^2} \frac{J_2}{5} \frac{7x^2 - 2y^2 - 5z^2}{r^2}, \quad (8) $$

where $\mu_E$ is the standard gravitational parameter of Europa (3202.72 km$^3$s$^{-2}$, Anderson, 1998).

For convenience, we decided in the following steps to fix $\Omega$ at the initial epoch $t_0$ so that the initial position vector of the probe $x(t_0)$ lies on the axis Europa-Jupiter, opposite to Jupiter, and that the initial velocity in x-direction vanishes. In the rotating frame, the initial state vector can then be written as $X(t_0) = (x_0, 0, 0, 0, \dot{y}_0, \dot{z}_0)$ (see Fig. 2). We subsequently performed a differential correction to refine the orbit to ensure its periodicity (Russell, 2006; Pavlak, 2013) in a two-steps iterative process described in Appendix A.

During the differential correction, the semi-major axis *a* can move by a few meters from the solution of Eq. (4). The orbits are also not precisely circular anymore, but their eccentricity remains below 0.003. One can note that the differential correction can be adapted to refine the orbit in more complex force models, for example by taking into account higher degrees and orders of the Europa gravity fields model.

Once the initial state vector $X(t_0)$ is corrected, the resulting propagated orbit is our reference orbit. It is important to note that these orbits do not have a reference epoch at this point: by fixing one, we get an orbital plane configuration with respect to the other bodies of interest in this study (Jupiter, Earth, Sun).

## 2.3. Choice of inclination

The inclination of the orbit plays a crucial role in the estimation of the gravity field. First, an induced polar gap in the ground coverage of Europa would limit the recovery of zonal and near-zonal spherical harmonic coefficients (van Gelderen and Koop, 1997). Additionally, the inclination influences the evolution of the probe's orbital plane geometry with respect to Earth due to orbit perturbation. This geometry is characterized by the angle between the orbital plane of the probe and the Earth direction ($\beta_{Earth}$).

The $\beta_{Earth}$ angle changes during a 3 months mission, but this change can be limited with a careful choice of the orbit inclination, in order to precisely study the impact of $\beta_{Earth}$ on the gravity field recovery.

Neglecting the obliquity of Europa (< 0.1°) and the inclination of Europa's orbital plane with respect to the ecliptic (< 2°), $\beta_{Earth}$ is related to the longitude of the Earth with respect to Europa in the ecliptic plane ($\Omega_{Earth}$), to the longitude of the ascending node of the probe around Europa ($\Omega$) and to its inclination $i$ (see Fig. 3) by

$$ \sin(\beta_{Earth}) = \sin(\Omega - \Omega_{Earth}) \cdot \sin(i), \quad (9) $$

where $\Omega_{Earth}$ depends on the Solar System configuration at the mission date. For a mission mid-2031, the time variation $\dot{\Omega}_{Earth}$ is approximately 0.1°/day. $\dot{\Omega}$ is caused by $J_2$ and by the influence of Jupiter as a third body and depends on the orbital characteristics (mainly the inclination, see
<caption>**Table 1**</caption>

m:R RGTO with their altitude and cycle intertrack fulfilling the orbit requirements. The two numbers in the 2nd column denote the minimum and maximum possible values of R. A lower R corresponds to a higher altitude, and a larger cycle intertrack (see Eq. (3)). m = 1 corresponds to 3.55 Earth days duration, and m = 26 corresponds to 92.33 Earth days cycle duration.

<table><thead><tr><th>m</th><th>R</th><th>Altitude</th><th>Cycle intertrack (at equator)</th></tr></thead><tbody><tr><td rowspan="2">1</td><td>38</td><td>181 km</td><td>258 km (9.5°)</td></tr><tr><td>40</td><td>123 km</td><td>245 km (9.0°)</td></tr><tr><td rowspan="2">2</td><td>75</td><td>197 km</td><td>131 km (4.8°)</td></tr><tr><td>81</td><td>109 km</td><td>121 km (4.4°)</td></tr><tr><td rowspan="2">3</td><td>113</td><td>192 km</td><td>87 km (3.2°)</td></tr><tr><td>122</td><td>104 km</td><td>80 km (3.0°)</td></tr><tr><td rowspan="2">26</td><td>973</td><td>199 km</td><td>10 km (0.37°)</td></tr><tr><td>1061</td><td>101 km</td><td>9 km (0.34°)</td></tr></tbody></table>

characteristics, $d_K$ is related to the Keplerian motion, $d_2$ is related to $J_2$ effect and $d_4$ is related to $J_2^2$ effect. Equation (4) (which is Eq. (10) in Cinelli et al., 2015) is obtained by fixing the other orbital elements (e, i, $\omega$), and is a function of the ratio $R/m$, of the physical parameters of Europa, and of the mean motion of Europa around Jupiter.

For each pair (m, R), we solve the degree-14 polynomial equation obtained by squaring Eq. (4) for the semi-major axis *a* while fixing *e*, *i* and *ω*. From the multiple solutions in our case only one falls within the range of studied altitudes. In this study, we do not consider elliptic orbits. We thus fix the *e* and *ω* to 0.

Equation (4) was derived by approximating third body perturbations, so that it only provides a first guess orbit. Such orbit then needs to be refined according to the full force model to ensure ground track repetition. In addition, the initial longitude of the ascending node, $\Omega$, has yet to be defined at this stage, which completes the characterization of the geometry of the orbital plane with respect to the third body Jupiter.

## 2.2. Orbit refinement

In order to precisely meet the condition for ground track repetition, we considered the Hill model (Lara and Russell, 2007) to refine our first guess orbits. This model takes into account the influence of Europa and Jupiter as point masses, plus the effect of the $J_2$ and $C_{22}$ gravity field coefficients of Europa.

We define the Jupiter-Europa rotating frame with the x-axis pointing from Jupiter to Europa, the y-axis in direction of the velocity vector of Europa around Jupiter, and the z-axis perpendicular to the orbital plane of Europa around Jupiter (see Fig. 2). In this frame centered around Europa, the spacecraft equations of motion are

$$ \ddot{x} + 2n_J \times \dot{x} = \nabla_x U_{tot}, \quad (5) $$

with

<figure><img src="image_5.png" /><figcaption>**Fig. 2.** Orbital plane of the probe (in red) at $t_0$ and $t_s$ with respect to Jupiter and Earth. The red dots represent crossing points with respect to Europa equatorial plane ($z_0 = z_s = 0$).</figcaption></figure>

$$ U_{tot} = \frac{n_J^2}{2} (x^2 + y^2) + U_E + U_J, \quad (6) $$

and where $x = (x, y, z)$ is the spacecraft position in the rotating frame, a dot denotes time derivatives, a bold variable denotes a vector, $n_J$ is the angular velocity vector of Europa around Jupiter, $n_J$ its magnitude, and $U_E$ and $U_J$ are the gravitational potentials of Europa and Jupiter, respectively. The gravitational potential of Jupiter is then defined by

$$ U_J = \frac{n_J^2}{2} (3x^2 - r^2), \quad (7) $$

where $r = \|x\|$. Assuming hydrostatic equilibrium for Europa ($J_2 = \frac{10}{3}C_{22}$), its gravitational potential can be separated into two contributions: the potential due to central point mass and the degree-2 gravity field coefficients potential $U_p$ defined so that

$$ U_E = \frac{\mu_E}{r} + U_p = \frac{\mu_E}{r} + \frac{\mu_E R_E^2}{r^2} \frac{J_2}{5} \frac{7x^2 - 2y^2 - 5z^2}{r^2}, \quad (8) $$

where $\mu_E$ is the standard gravitational parameter of Europa (3202.72 km$^3$s$^{-2}$, Anderson, 1998).

For convenience, we decided in the following steps to fix $\Omega$ at the initial epoch $t_0$ so that the initial position vector of the probe $x(t_0)$ lies on the axis Europa-Jupiter, opposite to Jupiter, and that the initial velocity in x-direction vanishes. In the rotating frame, the initial state vector can then be written as $X(t_0) = (x_0, 0, 0, 0, \dot{y}_0, \dot{z}_0)$ (see Fig. 2). We subsequently performed a differential correction to refine the orbit to ensure its periodicity (Russell, 2006; Pavlak, 2013) in a two-steps iterative process described in Appendix A.

During the differential correction, the semi-major axis *a* can move by a few meters from the solution of Eq. (4). The orbits are also not precisely circular anymore, but their eccentricity remains below 0.003. One can note that the differential correction can be adapted to refine the orbit in more complex force models, for example by taking into account higher degrees and orders of the Europa gravity fields model.

Once the initial state vector $X(t_0)$ is corrected, the resulting propagated orbit is our reference orbit. It is important to note that these orbits do not have a reference epoch at this point: by fixing one, we get an orbital plane configuration with respect to the other bodies of interest in this study (Jupiter, Earth, Sun).

## 2.3. Choice of inclination

The inclination of the orbit plays a crucial role in the estimation of the gravity field. First, an induced polar gap in the ground coverage of Europa would limit the recovery of zonal and near-zonal spherical harmonic coefficients (van Gelderen and Koop, 1997). Additionally, the inclination influences the evolution of the probe's orbital plane geometry with respect to Earth due to orbit perturbation. This geometry is characterized by the angle between the orbital plane of the probe and the Earth direction ($\beta_{Earth}$).

The $\beta_{Earth}$ angle changes during a 3 months mission, but this change can be limited with a careful choice of the orbit inclination, in order to precisely study the impact of $\beta_{Earth}$ on the gravity field recovery.

Neglecting the obliquity of Europa (< 0.1°) and the inclination of Europa's orbital plane with respect to the ecliptic (< 2°), $\beta_{Earth}$ is related to the longitude of the Earth with respect to Europa in the ecliptic plane ($\Omega_{Earth}$), to the longitude of the ascending node of the probe around Europa ($\Omega$) and to its inclination $i$ (see Fig. 3) by

$$ \sin(\beta_{Earth}) = \sin(\Omega - \Omega_{Earth}) \cdot \sin(i), \quad (9) $$

where $\Omega_{Earth}$ depends on the Solar System configuration at the mission date. For a mission mid-2031, the time variation $\dot{\Omega}_{Earth}$ is approximately 0.1°/day. $\dot{\Omega}$ is caused by $J_2$ and by the influence of Jupiter as a third body and depends on the orbital characteristics (mainly the inclination, see
<figure><img src="image_5.png" /><figcaption>Fig. 3. Geometry of the orbital plane of the probe with respect to the Earth characterised by the β<sub>Earth</sub> angle.</figcaption></figure>

Cinelli et al., 2015). It can be derived from the Gaussian perturbation equation (e.g., Beutler, 2005),

$$ \dot{\Omega} = \frac{r \sin(u)}{h \sin(i)} N, \qquad (10) $$

where r is the distance of the probe to Europa, u is the argument of latitude of the probe, h = ||x × ẋ|| is the norm of the angular momentum and N the cross-track component of the perturbing accelerations. In the Hill model, these accelerations are the influence of Jupiter and of the degree 2 gravity field of Europa, characterised by their potentials U<sub>J</sub> and U<sub>P</sub> (see Eqs. (7) and (8)) so that

$$ N = \nabla_x(U_J + U_P) \cdot e_N. \qquad (11) $$

For near polar and quasi circular orbits the time derivative of the mean longitude of the ascending node can be expressed as

$$ \dot{\Omega} = -\frac{3}{2n} \left( J_2 n^2 \left( \frac{R_E}{a} \right)^2 + \frac{n_J^2}{2} \right) \cos(i), \qquad (12) $$

with n the mean motion of the probe around Europa, a the semi-major axis of the probe's orbit and R<sub>E</sub> the radius of Europa (see details in Appendix B). The orbits investigated in this paper result in $\dot{\Omega} \in [-0.74, 0.74]^{\circ}$/day. More specifically, we obtain $\dot{\Omega} = -0.1^{\circ}$/day $\approx -\dot{\Omega}_E$ for an 88.6°-inclined orbit. This implies $\dot{\beta}_{Earth} = 0^{\circ}$/day, i.e., a quasi-fixed geometry of the orbital probe plane with respect to Earth. As mentioned in Sec. 2.2, the inclination i<sub>0</sub> can be fixed at t<sub>0</sub>, which will be different from the mean inclination i<sub>m</sub>. The relation between i<sub>0</sub> and i<sub>m</sub> can be derived from the Gaussian perturbation equation for the inclination i

$$ \frac{di}{dt} = \frac{r \cos(u)}{h} N. \qquad (13) $$

Similarly to $\dot{\Omega}$ (see Appendix B), the time derivative of the inclination for near polar and quasi circular orbits can be expressed as

$$ \frac{di}{dt} = -4K_i \cos^2(u) \sin(2u_J), \qquad (14) $$

with

$$ K_i = \frac{3}{4} \frac{\sin(i)}{h} \left( \frac{1}{2} n_J^2 a^2 + \frac{3}{5} J_2 n^2 R_E^2 \right), \qquad (15) $$

and u<sub>J</sub> = n<sub>J</sub>t - Ω the argument of latitude of Jupiter around Europa.

Considering dot<sub>J</sub> = n<sub>J</sub> - $\dot{\Omega}$ ≈ n<sub>J</sub> ≪ n, Eq. (14) can be integrated via the first order perturbation calculation method (i.e., by setting K<sub>i</sub> = K<sub>im</sub>), which yields

$$ i(t) = i_m + K_{im} \frac{\cos(2u_J)}{n_J}, \qquad (16) $$

with i<sub>m</sub> the mean inclination. Because the initial conditions of the probe's orbit are fixed at (x<sub>0</sub>, 0, 0, 0, 0, 0, 0, 0) in the rotating frame (see Fig. 2), we also have u<sub>J</sub>(t<sub>0</sub>) = 180°, providing a simple relation between the mean inclination i<sub>m</sub> and the inclination i<sub>0</sub> at the initial time t<sub>0</sub>, such that

$$ i_0 = i_m + \frac{K_{im}}{n_J}. \qquad (17) $$

Eq. (17) allows us to choose the mean inclination during the refinement of the orbit detailed in Sec. 2.2. For example, in order to get a 117 km altitude orbit with a mean inclination of 89°, the orbit design should target an orbit with an initial inclination of 89.9°.

## 2.4. Scenario definition

In order to consistently compare the scientific value of each orbit, we considered several “scenarios” from which the orbit can be propagated in the simulation environment of the BSW using an detailed force model.

A scenario is defined as the combination of a given reference orbit (see Sec. 2.2) and an initial location on the cycle grid at a reference starting epoch. For practical reasons, we restrict the possible initial locations to the equatorial plane of Europa (see Fig. 2). Since an m:R RGTO completes R orbit revolutions before overlapping, there are R distinct ascending crossing points on the equatorial plane, as depicted in Fig. 4.

Therefore, after the orbit refinement (see Sec. 2.2) we have a list of state vectors X(t<sub>k</sub>) = (x<sub>k</sub>, y<sub>k</sub>, 0, 0, 0, 0, 0, 0, 0) in the rotating frame, with k = 0, ..., R-1. Defining a scenario amounts to selecting one of these state vectors with k = s and assigning a reference epoch to t<sub>s</sub>. This then fixes the initial longitude of the ascending node Ω, as well as β<sub>Earth</sub> (see Fig. 3).

Once the scenario is defined, the state vector X(t<sub>s</sub>) is converted from the rotating frame, where the force field of the Hill model is defined, to the inertial frame. The orbit is propagated from this vector to a full ephemeris and force model using the planetary extension of the BSW (see Sec. 3.1). Because the orbits are designed in a simpler force model, the ground track repetition condition is no longer precisely met. Moreover, for longer integration times, the considerable perturbations of Jupiter would cause a short lifetime of the mission i.e., the eccentricity will increase until the orbiter crashes on the surface of Europa.

To prevent this, the spacecraft will have to undergo more or less frequent orbit maintenance maneuvers. Our way of taking into account

<figure><img src="image_6.png" /><figcaption>Fig. 4. Distribution of the ascending orbit crossing points with respect to the equatorial plane for a 1:40 RGTO (top) and a 3:118 RGTO (bottom). The larger red dots denote the locations of the probe c<sub>k</sub>(k = 0, 1, 2) at the start of the first three propagation arcs in the case of c<sub>0</sub> = 1.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Fig. 3. Geometry of the orbital plane of the probe with respect to the Earth characterised by the β<sub>Earth</sub> angle.</figcaption></figure>

Cinelli et al., 2015). It can be derived from the Gaussian perturbation equation (e.g., Beutler, 2005),

$$ \dot{\Omega} = \frac{r \sin(u)}{h \sin(i)} N, \qquad (10) $$

where r is the distance of the probe to Europa, u is the argument of latitude of the probe, h = ||x × ˙x|| is the norm of the angular momentum and N the cross-track component of the perturbing accelerations. In the Hill model, these accelerations are the influence of Jupiter and of the degree 2 gravity field of Europa, characterised by their potentials U<sub>J</sub> and U<sub>P</sub> (see Eqs. (7) and (8)) so that

$$ N = \nabla_x(U_J + U_P) \cdot e_N. \qquad (11) $$

For near polar and quasi circular orbits the time derivative of the mean longitude of the ascending node can be expressed as

$$ \dot{\Omega} = -\frac{3}{2n} \left( J_2 n^2 \left( \frac{R_E}{a} \right)^2 + \frac{n_J^2}{2} \right) \cos(i), \qquad (12) $$

with n the mean motion of the probe around Europa, a the semi-major axis of the probe's orbit and R<sub>E</sub> the radius of Europa (see details in Appendix B). The orbits investigated in this paper result in $\dot{\Omega} \in [-0.74, 0.74]^{\circ}$/day. More specifically, we obtain $\dot{\Omega} = -0.1^{\circ}$/day $\approx -\dot{\Omega}_E$ for an 88.6°-inclined orbit. This implies $\dot{\beta}_{Earth} = 0^{\circ}$/day, i.e., a quasi-fixed geometry of the orbital probe plane with respect to Earth. As mentioned in Sec. 2.2, the inclination i<sub>0</sub> can be fixed at t<sub>0</sub>, which will be different from the mean inclination i<sub>m</sub>. The relation between i<sub>0</sub> and i<sub>m</sub> can be derived from the Gaussian perturbation equation for the inclination i

$$ \frac{di}{dt} = \frac{r \cos(u)}{h} N. \qquad (13) $$

Similarly to $\dot{\Omega}$ (see Appendix B), the time derivative of the inclination for near polar and quasi circular orbits can be expressed as

$$ \frac{di}{dt} = -4K_i \cos^2(u) \sin(2u_J), \qquad (14) $$

with

$$ K_i = \frac{3}{4} \frac{\sin(i)}{h} \left( \frac{1}{2} n_J^2 a^2 + \frac{3}{5} J_2 n^2 R_E^2 \right), \qquad (15) $$

and u<sub>J</sub> = n<sub>J</sub>t - Ω the argument of latitude of Jupiter around Europa.

Considering ˙u<sub>J</sub> = n<sub>J</sub> - $\dot{\Omega}$ ≈ n<sub>J</sub> ≪ n, Eq. (14) can be integrated via the first order perturbation calculation method (i.e., by setting K<sub>i</sub> = K<sub>im</sub>), which yields

$$ i(t) = i_m + K_{im} \frac{\cos(2u_J)}{n_J}, \qquad (16) $$

with i<sub>m</sub> the mean inclination. Because the initial conditions of the probe's orbit are fixed at (x<sub>0</sub>, 0, 0, 0, ˙y<sub>0</sub>, ˙z<sub>0</sub>) in the rotating frame (see Fig. 2), we also have u<sub>J</sub>(t<sub>0</sub>) = 180°, providing a simple relation between the mean inclination i<sub>m</sub> and the inclination i<sub>0</sub> at the initial time t<sub>0</sub>, such that

$$ i_0 = i_m + \frac{K_{im}}{n_J}. \qquad (17) $$

Eq. (17) allows us to choose the mean inclination during the refinement of the orbit detailed in Sec. 2.2. For example, in order to get a 117 km altitude orbit with a mean inclination of 89°, the orbit design should target an orbit with an initial inclination of 89.9°.

## 2.4. Scenario definition

In order to consistently compare the scientific value of each orbit, we considered several “scenarios” from which the orbit can be propagated in the simulation environment of the BSW using an detailed force model.

A scenario is defined as the combination of a given reference orbit (see Sec. 2.2) and an initial location on the cycle grid at a reference starting epoch. For practical reasons, we restrict the possible initial locations to the equatorial plane of Europa (see Fig. 2). Since an m:R RGTO completes R orbit revolutions before overlapping, there are R distinct ascending crossing points on the equatorial plane, as depicted in Fig. 4.

Therefore, after the orbit refinement (see Sec. 2.2) we have a list of state vectors X(t<sub>k</sub>) = (x<sub>k</sub>, y<sub>k</sub>, 0, ˙x<sub>k</sub>, ˙y<sub>k</sub>, ˙z<sub>k</sub>) in the rotating frame, with k = 0, ..., R-1. Defining a scenario amounts to selecting one of these state vectors with k = s and assigning a reference epoch to t<sub>s</sub>. This then fixes the initial longitude of the ascending node Ω, as well as β<sub>Earth</sub> (see Fig. 3).

Once the scenario is defined, the state vector X(t<sub>s</sub>) is converted from the rotating frame, where the force field of the Hill model is defined, to the inertial frame. The orbit is propagated from this vector to a full ephemeris and force model using the planetary extension of the BSW (see Sec. 3.1). Because the orbits are designed in a simpler force model, the ground track repetition condition is no longer precisely met. Moreover, for longer integration times, the considerable perturbations of Jupiter would cause a short lifetime of the mission i.e., the eccentricity will increase until the orbiter crashes on the surface of Europa.

To prevent this, the spacecraft will have to undergo more or less frequent orbit maintenance maneuvers. Our way of taking into account

<figure><img src="image_7.png" /><figcaption>Fig. 4. Distribution of the ascending orbit crossing points with respect to the equatorial plane for a 1:40 RGTO (top) and a 3:118 RGTO (bottom). The larger red dots denote the locations of the probe c<sub>k</sub>(k = 0, 1, 2) at the start of the first three propagation arcs in the case of c<sub>0</sub> = 1.</figcaption></figure>
such maneuvers is as follows. Starting from the initial conditions at $c_0$ we propagate the orbit in arcs with a length of approximately one Europa nodal day $D_n$ (≈ 3.55 Earth days). For an *m*:R RGTO the probe orbits *R/m* times during $D_n$.

This number is an integer only for *m* = 1, i.e., after one Europa day and *R* revolutions the probe will be at the same position on the equatorial plane again, see Fig. 4 (top). For *m* ≠ 1 the length of the propagated arc is either shortened or increased with respect to one Europa day, implying that the probe ends up on the equatorial plane at a (different) crossing point $c_1$. In Fig. 4 (bottom), the example of a 3:118 RGTO is shown, where the first arc is propagated for 41 revolutions. Then, for the second arc, the initial conditions at $c_1$ from the reference orbit are used for further propagating the orbit for an integer number of times close to *R/m*. In the example of the 3:118 RGTO the second arc was chosen to have 39 revolutions to end up at the equatorial plane at the crossing point $c_2$ (Fig. 4, bottom). We choose a new arc *k* to end after $[k \cdot R/m]$ revolutions counting from $c_0$.

This is repeated in total *m* times until the probe has performed a total of *R* revolutions and ends at the same crossing point $c_0$ again. This in turn is repeated to cover the entire assumed mission duration of 3 months. Using new initial conditions from the reference orbit for each arc guarantees that the orbit remains stable and an approximate RGTO for the entire mission duration and thus mimics an orbit maintenance maneuvers approximately every Europa day. Our approach will introduce orbit discontinuities but we found them to be < 3 km, which we do not consider critical for gravity field recovery results at our target resolution. The accumulated velocity change over 3 months amounts to 50 m/s, which is a realistic Δ$v$ budget (Blanc et al., 2020).

### 3. Simulation in the Bernese GNSS software

Each selected orbit scenario is then analyzed in the planetary extension of the Bernese GNSS Software (Dach et al., 2015; Arnold et al., 2015; Bertone et al., 2021). First, the orbit is propagated every Europa day starting from 1-May-2031, based on a reference force model (see Sec. 3.1), which will constitute our ground truth for later comparisons. Then, realistic 2-way, X-band Doppler tracking measurements are generated along the orbit (see Sec. 3.2).

These measurements are then used to reconstruct the orbit and to estimate the gravity field in a standard multi-arc least-squares (LS) process following the Celestial Mechanic Approach (Beutler et al., 2010), as described in Sec. 3.3. The quality of the orbit solution is discussed in Sec. 3.4.

#### 3.1. Force model

In a quasi inertial frame, the spacecraft equations of motion can be written as

$$ \ddot{r} = -\mu_E \frac{r}{r^3} + f(t, r, \dot{r}, q_1, \dots, q_d), \quad (18) $$

where *r* is the position vector of the spacecraft center of mass with respect to the center of mass of Europa, *r* its absolute value, a dot denotes time derivative and *f* collects all perturbation accelerations beyond the central

term (acceleration due to central point mass only). The $q_i$, $i = 1, \dots, d$, denote the parameters used to represent different kinds of orbit perturbations (e.g., empirical orbit parameters, gravity field parameters ...). Table 2 summarizes the accelerations *f* that we considered for orbit propagation. The integration of the equations of motion was performed via the collocation method described in Beutler (2005).

As part of this force model, we considered a synthetic gravity field model for Europa, whose potential *V* can be decomposed in a spherical harmonic expansion (Kaula and Street, 1967) as

$$ V(r, \lambda, \varphi) = \frac{\mu_E}{r} \sum_{n=2}^{n_{\max}} \sum_{m=0}^{n} \left( \frac{R_E}{r} \right)^n \bar{P}_{nm}(\sin \varphi) (\bar{C}_{nm} \cos m\lambda + \bar{S}_{nm} \sin m\lambda), \quad (19) $$

where $\lambda$ and $\varphi$ are longitude and latitude in a Europa-fixed reference system, $\bar{P}_{nm}$ are the fully normalized Legendre functions of degree *n* and order *m*, and $\bar{C}_{nm}$ and $\bar{S}_{nm}$ are the spherical harmonic coefficients.

We used the estimates of Anderson (1998) for coefficients up to degree and order (d/o) 2, while higher d/o were derived from the Moon's gravity field GRGM900C (Lemoine et al., 2014), appropriately scaled by the squared ratio of the surface gravity of Europa and of the Moon ($g_i = GM_i/R_i^2$, where *i* stands for Europa and the Moon), up to degree and order $n_{\max} = 90$. We label this synthetic reference gravity field as EUR-GLMo in the following.

Tidal deformations induced on Europa by Jupiter change the degree-2 gravity field coefficients according to Petit and Luzum (2010) as

$$ \Delta \bar{C}_{2m} - i \Delta \bar{S}_{2m} = \frac{k_{2m}}{5} \frac{\mu_J}{\mu_E} \left( \frac{R_E}{r_J} \right)^3 \bar{P}_{2m}(\sin \varphi_J) e^{-im\lambda_J}, \quad (20) $$

where $k_{2m}$ is the Love number for degree 2 and order *m*, $r_J$, $\lambda_J$ and $\varphi_J$ denote the spherical coordinates of Jupiter in a Europa-fixed reference system. $\mu_J$ denote the gravity constant of Jupiter. In this study, we assumed $k_{2m} = k_2 = 0.257$ (Wahr et al., 2006; Mazarico et al., 2015).

For the third-body attractions we considered the Sun, the planets and the other Galilean moons as point masses. As the influence of Jupiter is considerable, we also included Jupiter's zonal gravity field coefficients up to degree 6. Non-gravitational forces were not considered in this study, as their impact is negligible for our goals.

#### 3.2. Synthetic Doppler observables

Commonly used 2-way X-band Doppler observables are simulated along with the propagated orbit as measurements from Deep Space Network (DSN, Asmar and Renzetti, 1993) stations. The observables are computed following the formulation by Moyer (2003), under the following conditions: minimum spacecraft elevation of 10° above the horizon and eclipses from the Sun, Jupiter and Europa. Tropospheric and ionospheric delays are not considered. In this study, we consider the tracking of the 3 DSN stations (Canberra, Goldstone and Madrid) as baseline to reduce the number of free parameters and we tested that reducing the coverage to only one station doesn't significantly affect our results for the scope of this work. Light times are computed considering the Shapiro gravitational delay due to gravitational perturbation of the Sun, Earth and Jupiter on light propagation.

<table><caption>**Table 2**
Summary of accelerations included in our force model, see Eq. (18).</caption><thead><tr><th>Acceleration</th><th>Magnitude (m/s<sup>2</sup>)</th><th>Source</th></tr></thead><tbody><tr><td>Europa central term</td><td>1</td><td>Anderson (1998)</td></tr><tr><td>Higher degrees of Europa gravity field</td><td>1.10<sup>-3</sup></td><td>Anderson (1998) &amp; synthetic field</td></tr><tr><td>Jupiter-induced Tides</td><td>2.10<sup>-4</sup></td><td>Mazarico et al. (2015)</td></tr><tr><td>Jupiter central term</td><td>9.10<sup>-4</sup></td><td>JUP310 (Jacobson, 2013)</td></tr><tr><td>Higher degrees of Jupiter gravity field</td><td>7.10<sup>-5</sup></td><td>JUP310 (Jacobson, 2013)</td></tr><tr><td>Other Galilean moons (central terms)</td><td>1.10<sup>-8</sup>–1.10<sup>-6</sup></td><td>JUP310 (Jacobson, 2013)</td></tr><tr><td>Sun</td><td>7.10<sup>-10</sup></td><td>DE430 (Folkner et al., 2014)</td></tr><tr><td>Other planets</td><td>2.10<sup>-14</sup></td><td>DE430 (Folkner et al., 2014)</td></tr></tbody></table>
such maneuvers is as follows. Starting from the initial conditions at $c_0$ we propagate the orbit in arcs with a length of approximately one Europa nodal day $D_n$ (≈ 3.55 Earth days). For an *m*:R RGTO the probe orbits *R/m* times during $D_n$.

This number is an integer only for *m* = 1, i.e., after one Europa day and *R* revolutions the probe will be at the same position on the equatorial plane again, see Fig. 4 (top). For *m* ≠ 1 the length of the propagated arc is either shortened or increased with respect to one Europa day, implying that the probe ends up on the equatorial plane at a (different) crossing point $c_1$. In Fig. 4 (bottom), the example of a 3:118 RGTO is shown, where the first arc is propagated for 41 revolutions. Then, for the second arc, the initial conditions at $c_1$ from the reference orbit are used for further propagating the orbit for an integer number of times close to *R/m*. In the example of the 3:118 RGTO the second arc was chosen to have 39 revolutions to end up at the equatorial plane at the crossing point $c_2$ (Fig. 4, bottom). We choose a new arc *k* to end after $[k \cdot R/m]$ revolutions counting from $c_0$.

This is repeated in total *m* times until the probe has performed a total of *R* revolutions and ends at the same crossing point $c_0$ again. This in turn is repeated to cover the entire assumed mission duration of 3 months. Using new initial conditions from the reference orbit for each arc guarantees that the orbit remains stable and an approximate RGTO for the entire mission duration and thus mimics an orbit maintenance maneuvers approximately every Europa day. Our approach will introduce orbit discontinuities but we found them to be < 3 km, which we do not consider critical for gravity field recovery results at our target resolution. The accumulated velocity change over 3 months amounts to 50 m/s, which is a realistic Δ$v$ budget (Blanc et al., 2020).

### 3. Simulation in the Bernese GNSS software

Each selected orbit scenario is then analyzed in the planetary extension of the Bernese GNSS Software (Dach et al., 2015; Arnold et al., 2015; Bertone et al., 2021). First, the orbit is propagated every Europa day starting from 1-May-2031, based on a reference force model (see Sec. 3.1), which will constitute our ground truth for later comparisons. Then, realistic 2-way, X-band Doppler tracking measurements are generated along the orbit (see Sec. 3.2).

These measurements are then used to reconstruct the orbit and to estimate the gravity field in a standard multi-arc least-squares (LS) process following the Celestial Mechanic Approach (Beutler et al., 2010), as described in Sec. 3.3. The quality of the orbit solution is discussed in Sec. 3.4.

#### 3.1. Force model

In a quasi inertial frame, the spacecraft equations of motion can be written as

$$ \ddot{r} = -\mu_E \frac{r}{r^3} + f(t, r, \dot{r}, q_1, \dots, q_d), \quad (18) $$

where *r* is the position vector of the spacecraft center of mass with respect to the center of mass of Europa, *r* its absolute value, a dot denotes time derivative and *f* collects all perturbation accelerations beyond the central

term (acceleration due to central point mass only). The $q_i$, $i = 1, \dots, d$, denote the parameters used to represent different kinds of orbit perturbations (e.g., empirical orbit parameters, gravity field parameters ...). Table 2 summarizes the accelerations *f* that we considered for orbit propagation. The integration of the equations of motion was performed via the collocation method described in Beutler (2005).

As part of this force model, we considered a synthetic gravity field model for Europa, whose potential *V* can be decomposed in a spherical harmonic expansion (Kaula and Street, 1967) as

$$ V(r, \lambda, \varphi) = \frac{\mu_E}{r} \sum_{n=2}^{n_{\max}} \sum_{m=0}^{n} \left( \frac{R_E}{r} \right)^n \bar{P}_{nm}(\sin \varphi) (\bar{C}_{nm} \cos m\lambda + \bar{S}_{nm} \sin m\lambda), \quad (19) $$

where $\lambda$ and $\varphi$ are longitude and latitude in a Europa-fixed reference system, $\bar{P}_{nm}$ are the fully normalized Legendre functions of degree *n* and order *m*, and $\bar{C}_{nm}$ and $\bar{S}_{nm}$ are the spherical harmonic coefficients.

We used the estimates of Anderson (1998) for coefficients up to degree and order (d/o) 2, while higher d/o were derived from the Moon's gravity field GRGM900C (Lemoine et al., 2014), appropriately scaled by the squared ratio of the surface gravity of Europa and of the Moon ($g_i = GM_i/R_i^2$, where *i* stands for Europa and the Moon), up to degree and order $n_{\max} = 90$. We label this synthetic reference gravity field as EUR-GLMo in the following.

Tidal deformations induced on Europa by Jupiter change the degree-2 gravity field coefficients according to Petit and Luzum (2010) as

$$ \Delta \bar{C}_{2m} - i \Delta \bar{S}_{2m} = \frac{k_{2m}}{5} \frac{\mu_J}{\mu_E} \left( \frac{R_E}{r_J} \right)^3 \bar{P}_{2m}(\sin \varphi_J) e^{-im\lambda_J}, \quad (20) $$

where $k_{2m}$ is the Love number for degree 2 and order *m*, $r_J$, $\lambda_J$ and $\varphi_J$ denote the spherical coordinates of Jupiter in a Europa-fixed reference system. $\mu_J$ denote the gravity constant of Jupiter. In this study, we assumed $k_{2m} = k_2 = 0.257$ (Wahr et al., 2006; Mazarico et al., 2015).

For the third-body attractions we considered the Sun, the planets and the other Galilean moons as point masses. As the influence of Jupiter is considerable, we also included Jupiter's zonal gravity field coefficients up to degree 6. Non-gravitational forces were not considered in this study, as their impact is negligible for our goals.

#### 3.2. Synthetic Doppler observables

Commonly used 2-way X-band Doppler observables are simulated along with the propagated orbit as measurements from Deep Space Network (DSN, Asmar and Renzetti, 1993) stations. The observables are computed following the formulation by Moyer (2003), under the following conditions: minimum spacecraft elevation of 10° above the horizon and eclipses from the Sun, Jupiter and Europa. Tropospheric and ionospheric delays are not considered. In this study, we consider the tracking of the 3 DSN stations (Canberra, Goldstone and Madrid) as baseline to reduce the number of free parameters and we tested that reducing the coverage to only one station doesn't significantly affect our results for the scope of this work. Light times are computed considering the Shapiro gravitational delay due to gravitational perturbation of the Sun, Earth and Jupiter on light propagation.

<table><caption>**Table 2**
Summary of accelerations included in our force model, see Eq. (18).</caption><thead><tr><th>Acceleration</th><th>Magnitude (m/s<sup>2</sup>)</th><th>Source</th></tr></thead><tbody><tr><td>Europa central term</td><td>1</td><td>Anderson (1998)</td></tr><tr><td>Higher degrees of Europa gravity field</td><td>1.10<sup>-3</sup></td><td>Anderson (1998) &amp; synthetic field</td></tr><tr><td>Jupiter-induced Tides</td><td>2.10<sup>-4</sup></td><td>Mazarico et al. (2015)</td></tr><tr><td>Jupiter central term</td><td>9.10<sup>-4</sup></td><td>JUP310 (Jacobson, 2013)</td></tr><tr><td>Higher degrees of Jupiter gravity field</td><td>7.10<sup>-5</sup></td><td>JUP310 (Jacobson, 2013)</td></tr><tr><td>Other Galilean moons (central terms)</td><td>1.10<sup>-8</sup>–1.10<sup>-6</sup></td><td>JUP310 (Jacobson, 2013)</td></tr><tr><td>Sun</td><td>7.10<sup>-10</sup></td><td>DE430 (Folkner et al., 2014)</td></tr><tr><td>Other planets</td><td>2.10<sup>-14</sup></td><td>DE430 (Folkner et al., 2014)</td></tr></tbody></table>
The computation of Doppler observables is very sensitive to numerical noise. Already for the GRAIL data analysis in the BSW presented by Bertone et al. (2021), the handling of observation epochs was treated by separating the fractional part of the day from its integer part to increase numerical precision with double precision (64-bit) floats. However, when dealing with the tracking of planetary probes in the outer Solar System, the light time between the probe and a ground antenna is increased from a few seconds at the Moon to 30–50 min in the case of the Jupiter system. Consequently, we computed the light times in quadruple (128-bit) precision to avoid the occurrence of numerical errors. Results are then treated in double precision in the rest of the code, so that the performance is not heavily affected. In our controlled simulation environment, we reduced the orbit fit residual Root Mean Square (RMS) error by one order of magnitude (from $2.0 \times 10^{-3}$ mm/s to $1.3 \times 10^{-4}$ mm/s) by using localized quadruple precision to compute Doppler observables when no external source of noise are introduced. Following a numerical error model by Zannoni and Tortora (2013), we can assess an upper bound of $1.6 \times 10^{-3}$ mm/s for the numerical noise introduced by our orbit determination software for the tracking of probes in the Jovian system.

We then add Gaussian white noise consistent with 2-way X-band tracking data, simulated for a total of 3 months, with $\sigma_{obs} = 0.10$ mm/s (1-way) at 60s integration time, to encompass the major relevant noise sources (interplanetary plasma, troposphere, etc., Asmar et al., 2005).

### 3.3. Estimation process

Orbit and gravity field parameters are estimated following the Celestial Mechanics Approach (Beutler et al., 2010). As described in Arnold et al. (2015) and Bertone et al. (2021), this procedure is divided into an initial arc-wise orbit determination, followed by a generalized orbit improvement, including both orbit and geodetic parameters.

While the simulated orbits are generated in arcs of length close to 1 Europa day (see Sec. 2.4), in the orbit reconstruction we estimate 3 sets of initial conditions per such arc, for a total of 78 initial conditions for 26 Europa days ($\approx$ 3 months). This leads to (variable) arc lengths between 27 h and 30 h. For each orbit scenario, the first two arcs have the same length (28 h or 29 h), while the last ends at the beginning of the subsequent simulation arc (right after the maneuver), which makes its duration different but close to the other two arcs.

We add uncertainties with a standard deviation of $\sigma_{pos} = 50$ m, $\sigma_{vel} = 1$ mm/s to the initial position and velocity of the a priori orbits for each estimation arc. The equations of motion (18) are then numerically integrated based on the fixed background force model detailed in Sec. 3.1. The partial derivatives of the orbit with respect to each initial Keplerian osculating elements are obtained by numerically integrating the variational equations, and the ones with respect to the other estimated parameters are computed by solving definite integrals as outlined in Beutler (2005).

This allows for a least squares adjustment of arc-specific parameters (e.g., the six initial Keplerian osculating elements) to minimize the Doppler residuals (see Sec. 3.4). This process is repeated for at least 4 iterations, and until the relative change of RMS of Doppler residuals between one iteration and the next is below 0.5%.

Based on the updated orbits, the equations of motion and the variational equations are once more integrated to compute the partial derivatives with respect to both arc-specific and global parameters, thus allowing us to set up the normal equation systems (NEQs) from the Doppler observations for all parameters.

These NEQs are set up individually for each arc and then stacked for 26 Europa days. Arc-specific parameters are pre-eliminated (see, e.g., Dach et al., 2015) prior to NEQ stacking and the accumulated NEQ is inverted to retrieve corrections to the global parameters.

Table 3 summarizes our set of estimated parameters, with optional constraints. When using pseudo-stochastic pulses (instantaneous velocity changes, see Sect. 5.1), we set them every 60 min in all 3 directions

<caption>**Table 3**</caption>

Summary of the estimated parameters. Osculating orbital elements and $k_2$ were freely estimated ($k_{20} = k_{21} = k_{22}$). The Kaula constraint was only used for part of the solutions.

<table><thead><tr><th rowspan="2">Parameter</th><th colspan="2">Number of parameters</th><th rowspan="2">Constraint</th></tr><tr><th>Arc</th><th>Total</th></tr></thead><tbody><tr><td>Osculating elements</td><td>6</td><td>468</td><td>—</td></tr><tr><td>Pseudo-stochastic pulses</td><td>≈ 80</td><td>6384</td><td>1 mm/s</td></tr><tr><td>$k_2$ Love number</td><td>—</td><td>1</td><td>—</td></tr><tr><td>Gravity field coefficients</td><td>—</td><td>8278</td><td>Kaula ($K = 0.5$)</td></tr></tbody></table>

(radial, along-track and cross-track) and we constrained their amplitude to 1 mm/s. We also considered a Kaula law for some solutions to constrain gravity field coefficients to zero according to their degree $n$ with a variance of

$$ \sigma_K(n) = \frac{K}{n^2}. \qquad (21) $$

This results in tighter constraints on the high degree coefficients (which might diverge due to the limited by data coverage) than on the low-degree coefficients. K was empirically chosen to only prevent the high degree gravity field coefficients from diverging, and we were careful not to constrain the estimated parameters to their expected value.

### 3.4. Orbit determination quality

The arc-specific parameters are first estimated for every arc separately using the full degree-90 gravity field (EURGLMo) as a priori information. This first orbit recovery is repeated for at least 4 iterations until convergence is achieved. Using realistic Doppler observables described in Sec. 3.2 to improve the a priori orbits described in Sec. 3.3 allows to reduce 2-way Doppler residuals below 5 mHz when estimating only osculating elements, which corresponds to the applied Doppler noise (0.1 mm/s 1-way).

The simulation environment allows us to compare the reconstructed orbits with the true (simulated) orbits and to obtain a direct quality measure also on the orbit level.

The orbit differences are analyzed in a frame ($e_1, e_2, e_3$) related to the line of sight from the Earth, because of its relevance when using Doppler observations to reconstruct the orbit (Bonanno and Milani, 2002). $e_1$ is the unit vector from Europa to the Earth, $e_2$ is the normalized cross

<figure><img src="image_7.png" /><figcaption>**Fig. 5.** RMS of orbit differences for each arc in the ($e_1, e_2, e_3$) frame for a 5:197 RGTO ($i = 89^\circ, h = 134$ km, $\beta_{Earth} = -6^\circ$, using EURGLMo as gravity field model), with $\Delta = \| (e_1, e_2, e_3) \|$.</figcaption></figure>
The computation of Doppler observables is very sensitive to numerical noise. Already for the GRAIL data analysis in the BSW presented by Bertone et al. (2021), the handling of observation epochs was treated by separating the fractional part of the day from its integer part to increase numerical precision with double precision (64-bit) floats. However, when dealing with the tracking of planetary probes in the outer Solar System, the light time between the probe and a ground antenna is increased from a few seconds at the Moon to 30–50 min in the case of the Jupiter system. Consequently, we computed the light times in quadruple (128-bit) precision to avoid the occurrence of numerical errors. Results are then treated in double precision in the rest of the code, so that the performance is not heavily affected. In our controlled simulation environment, we reduced the orbit fit residual Root Mean Square (RMS) error by one order of magnitude (from $2.0 \times 10^{-3}$ mm/s to $1.3 \times 10^{-4}$ mm/s) by using localized quadruple precision to compute Doppler observables when no external source of noise are introduced. Following a numerical error model by Zannoni and Tortora (2013), we can assess an upper bound of $1.6 \times 10^{-3}$ mm/s for the numerical noise introduced by our orbit determination software for the tracking of probes in the Jovian system.

We then add Gaussian white noise consistent with 2-way X-band tracking data, simulated for a total of 3 months, with $\sigma_{obs} = 0.10$ mm/s (1-way) at 60s integration time, to encompass the major relevant noise sources (interplanetary plasma, troposphere, etc., Asmar et al., 2005).

### 3.3. Estimation process

Orbit and gravity field parameters are estimated following the Celestial Mechanics Approach (Beutler et al., 2010). As described in Arnold et al. (2015) and Bertone et al. (2021), this procedure is divided into an initial arc-wise orbit determination, followed by a generalized orbit improvement, including both orbit and geodetic parameters.

While the simulated orbits are generated in arcs of length close to 1 Europa day (see Sec. 2.4), in the orbit reconstruction we estimate 3 sets of initial conditions per such arc, for a total of 78 initial conditions for 26 Europa days ($\approx$ 3 months). This leads to (variable) arc lengths between 27 h and 30 h. For each orbit scenario, the first two arcs have the same length (28 h or 29 h), while the last ends at the beginning of the subsequent simulation arc (right after the maneuver), which makes its duration different but close to the other two arcs.

We add uncertainties with a standard deviation of $\sigma_{pos} = 50$ m, $\sigma_{vel} = 1$ mm/s to the initial position and velocity of the a priori orbits for each estimation arc. The equations of motion (18) are then numerically integrated based on the fixed background force model detailed in Sec. 3.1. The partial derivatives of the orbit with respect to each initial Keplerian osculating elements are obtained by numerically integrating the variational equations, and the ones with respect to the other estimated parameters are computed by solving definite integrals as outlined in Beutler (2005).

This allows for a least squares adjustment of arc-specific parameters (e.g., the six initial Keplerian osculating elements) to minimize the Doppler residuals (see Sec. 3.4). This process is repeated for at least 4 iterations, and until the relative change of RMS of Doppler residuals between one iteration and the next is below 0.5%.

Based on the updated orbits, the equations of motion and the variational equations are once more integrated to compute the partial derivatives with respect to both arc-specific and global parameters, thus allowing us to set up the normal equation systems (NEQs) from the Doppler observations for all parameters.

These NEQs are set up individually for each arc and then stacked for 26 Europa days. Arc-specific parameters are pre-eliminated (see, e.g., Dach et al., 2015) prior to NEQ stacking and the accumulated NEQ is inverted to retrieve corrections to the global parameters.

Table 3 summarizes our set of estimated parameters, with optional constraints. When using pseudo-stochastic pulses (instantaneous velocity changes, see Sect. 5.1), we set them every 60 min in all 3 directions

<caption>**Table 3**</caption>

Summary of the estimated parameters. Osculating orbital elements and $k_2$ were freely estimated ($k_{20} = k_{21} = k_{22}$). The Kaula constraint was only used for part of the solutions.

<table><thead><tr><th rowspan="2">Parameter</th><th colspan="2">Number of parameters</th><th rowspan="2">Constraint</th></tr><tr><th>Arc</th><th>Total</th></tr></thead><tbody><tr><td>Osculating elements</td><td>6</td><td>468</td><td>—</td></tr><tr><td>Pseudo-stochastic pulses</td><td>≈ 80</td><td>6384</td><td>1 mm/s</td></tr><tr><td>$k_2$ Love number</td><td>—</td><td>1</td><td>—</td></tr><tr><td>Gravity field coefficients</td><td>—</td><td>8278</td><td>Kaula (K = 0.5)</td></tr></tbody></table>

(radial, along-track and cross-track) and we constrained their amplitude to 1 mm/s. We also considered a Kaula law for some solutions to constrain gravity field coefficients to zero according to their degree $n$ with a variance of

$$ \sigma_K(n) = \frac{K}{n^2}. \qquad (21) $$

This results in tighter constraints on the high degree coefficients (which might diverge due to the limited by data coverage) than on the low-degree coefficients. K was empirically chosen to only prevent the high degree gravity field coefficients from diverging, and we were careful not to constrain the estimated parameters to their expected value.

### 3.4. Orbit determination quality

The arc-specific parameters are first estimated for every arc separately using the full degree-90 gravity field (EURGLMO) as a priori information. This first orbit recovery is repeated for at least 4 iterations until convergence is achieved. Using realistic Doppler observables described in Sec. 3.2 to improve the a priori orbits described in Sec. 3.3 allows to reduce 2-way Doppler residuals below 5 mHz when estimating only osculating elements, which corresponds to the applied Doppler noise (0.1 mm/s 1-way).

The simulation environment allows us to compare the reconstructed orbits with the true (simulated) orbits and to obtain a direct quality measure also on the orbit level.

The orbit differences are analyzed in a frame ($e_1, e_2, e_3$) related to the line of sight from the Earth, because of its relevance when using Doppler observations to reconstruct the orbit (Bonanno and Milani, 2002). $e_1$ is the unit vector from Europa to the Earth, $e_2$ is the normalized cross

<figure><img src="image_8.png" /><figcaption>**Fig. 5.** RMS of orbit differences for each arc in the ($e_1, e_2, e_3$) frame for a 5:197 RGTO ($i = 89^\circ, h = 134$ km, $\beta_{Earth} = -6^\circ$, using EURGLMO as gravity field model), with $\Delta = \| (e_1, e_2, e_3) \|$.</figcaption></figure>
product of $e_1$ with the position vector of the probe with respect to Europa, and $e_3$ is the cross product of $e_1$ and $e_2$.

In this frame, the $e_2$ direction, perpendicular to the line of sight and to the position vector, ends up being the worst determined direction (up to 100 times worse than along $e_1$ and $e_3$, see Fig. 5), regardless of the different parameters considered in Sec. 4. However, the quality of the recovery (< 1.3 m) is still at a level acceptable for most mission goals.

Using a perfect a priori force model and estimating only the oscillating orbital elements, the orbit can be recovered with a precision below 2 m. However, when co-estimating gravity field parameters, the quality of the orbit is degraded (see Fig. 6, bottom). The main orbit differences are still in the $e_2$ direction, but the RMS differences increase to a maximum of 6.3m. On the other hand, the Doppler residuals are reduced, which is a consequence of enlarging the parameter space by estimating gravity field parameters (see Fig. 6, top).

## 4. Impact on gravity field recovery

In this section, we investigate the quality of gravity field recovery for

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 6. RMS of Doppler residuals (top) and orbit differences (bottom) for every arc in the case of a 5:197 RGTO ($i = 89°$, $h = 134$ km, $\beta_{Earth} = -6°$). The red curve is the result of the fit considering a perfect a priori force model, and estimating only orbital parameters. The green curve corresponds to the estimation of the orbital elements and of the gravity field parameters up to d/o 90.</figcaption></figure>

multiple orbit scenarios defined in Sec. 2.4. We carefully analyze the influence of a set of parameters on the gravity field recovery, by isolating each of them as much as possible from the influence of the others. For computational time reasons, we use a perfect force model (i.e., the unperturbed simulated d/o 90 field) as a priori information for all our tests, to avoid iterating on the gravity field solution. We later show (in Sec. 5) for a chosen scenario that the gravity field can be recovered to a comparable level when starting from a degree-2 a priori gravity field (corresponding to our current knowledge), albeit for the price of more iterations.

We use the weighted RMS of the geoid height differences with respect to the reference gravity field as a global quality assessment. It is defined as

$$ (\Delta g)_{WRMS} = \sqrt{\frac{\sum_{\phi, \lambda} \cos(\phi) \Delta g_{\phi, \lambda}^2}{\sum_{\phi, \lambda} \cos(\phi)}}, \qquad (22) $$

where $\Delta g_{\phi, \lambda}$ is the difference of the geoid height at latitude $\phi$ and longitude $\lambda$ between a given gravity field solution and the reference gravity field. Gravity field solutions can also be compared in terms of difference degree amplitudes

$$ \Delta_n = \sqrt{\frac{1}{2n+1} \sum_{m=0}^{n} (\Delta \bar{C}_{nm}^2 + \Delta \bar{S}_{nm}^2)}, \qquad (23) $$

where $\Delta \bar{C}_{nm}$ and $\Delta \bar{S}_{nm}$ are the differences of the respective spherical harmonic coefficients. We also evaluate the error degree amplitudes, which are obtained from Eq. (23) by replacing the coefficient differences with the formal errors of the estimated coefficients.

In all cases presented here, the Love number $k_2$ is co-estimated, resulting in uncertainties < 1.3 × 10⁻⁵ (< 0.005% relative to the expected signal), which would allow to distinguish between an ocean-bearing and ocean-free Europa (Wu et al., 2001; Wahr et al., 2006). However, we found that the different scenarios investigated in this section only have a marginal impact on the recovery of $k_2$, so that we will not further discuss this parameter.

### 4.1. Ground tracks repetition

The repetition cycle of an RGTO directly impacts the ground surface coverage, as mentioned in Sec. 2.1. The shorter the repetition cycle, the larger the cycle intertrack (see Table 1), resulting in a lower the spatial resolution of the ground tracks, which can limit the resolution of the estimated gravity field solution.

We derive from Eq. (19) the full wavelength of a degree-n gravity field at the equator as $\lambda_E = \frac{2\pi R_E}{n}$. According to the Nyquist criterion, a sampling of $\frac{\lambda_E}{2}$ allows to recover down to the shortest wavelength of $\frac{2\pi R_E}{n}$. In other words, to avoid aliasing (Sneeuw, 2000) when estimating a degree-n gravity field, the spatial resolution (or cycle intertrack $\delta_{eq}$) of the ground tracks of an m:R RGTO (see Eq. (3)), must be larger than this minimum sampling, i.e., $\delta_{eq} \ge \frac{\lambda_E}{2}$. Hence, we can expect the maximum resolvable degree $n_{max}$ of a recovered gravity field from an m:R RGTO to be given by

$$ n_{max} = \frac{R}{2}. \qquad (24) $$

It can be deduced from this equation that for $R \ge 180$, we can estimate the gravity field coefficients up to d/o 90 without any aliasing. For an orbiter with an altitude between 100 km and 200 km, this corresponds to $m \ge 5$. As an example, for an orbit with an altitude $h$ close to 134 km and $\beta_{Earth} \approx 67°$, the gravity field can be estimated unconstrained up to d/o 62, for every $m \ge 5$ (see $m = 5$ and $m = 26$ in Fig. 7). Indeed, decreasing the repetition rate so that the ground tracks coverage becomes denser does not improve the quality of the gravity field solution.
product of $e_1$ with the position vector of the probe with respect to Europa, and $e_3$ is the cross product of $e_1$ and $e_2$.

In this frame, the $e_2$ direction, perpendicular to the line of sight and to the position vector, ends up being the worst determined direction (up to 100 times worse than along $e_1$ and $e_3$, see Fig. 5), regardless of the different parameters considered in Sec. 4. However, the quality of the recovery (< 1.3 m) is still at a level acceptable for most mission goals.

Using a perfect a priori force model and estimating only the oscillating orbital elements, the orbit can be recovered with a precision below 2 m. However, when co-estimating gravity field parameters, the quality of the orbit is degraded (see Fig. 6, bottom). The main orbit differences are still in the $e_2$ direction, but the RMS differences increase to a maximum of 6.3m. On the other hand, the Doppler residuals are reduced, which is a consequence of enlarging the parameter space by estimating gravity field parameters (see Fig. 6, top).

## 4. Impact on gravity field recovery

In this section, we investigate the quality of gravity field recovery for

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 6. RMS of Doppler residuals (top) and orbit differences (bottom) for every arc in the case of a 5:197 RGTO ($i = 89^\circ$, $h = 134$ km, $\beta_{Earth} = -6^\circ$). The red curve is the result of the fit considering a perfect a priori force model, and estimating only orbital parameters. The green curve corresponds to the estimation of the orbital elements and of the gravity field parameters up to d/o 90.</figcaption></figure>

multiple orbit scenarios defined in Sec. 2.4. We carefully analyze the influence of a set of parameters on the gravity field recovery, by isolating each of them as much as possible from the influence of the others. For computational time reasons, we use a perfect force model (i.e., the unperturbed simulated d/o 90 field) as a priori information for all our tests, to avoid iterating on the gravity field solution. We later show (in Sec. 5) for a chosen scenario that the gravity field can be recovered to a comparable level when starting from a degree-2 a priori gravity field (corresponding to our current knowledge), albeit for the price of more iterations.

We use the weighted RMS of the geoid height differences with respect to the reference gravity field as a global quality assessment. It is defined as

$$ (\Delta g)_{WRMS} = \sqrt{\frac{\sum_{\varphi, \lambda} \cos(\varphi) \Delta g_{\varphi, \lambda}^2}{\sum_{\varphi, \lambda} \cos(\varphi)}}, \qquad (22) $$

where $\Delta g_{\varphi, \lambda}$ is the difference of the geoid height at latitude $\varphi$ and longitude $\lambda$ between a given gravity field solution and the reference gravity field. Gravity field solutions can also be compared in terms of difference degree amplitudes

$$ \Delta_n = \sqrt{\frac{1}{2n+1} \sum_{m=0}^{n} (\Delta \bar{C}_{nm}^2 + \Delta \bar{S}_{nm}^2)}, \qquad (23) $$

where $\Delta \bar{C}_{nm}$ and $\Delta \bar{S}_{nm}$ are the differences of the respective spherical harmonic coefficients. We also evaluate the error degree amplitudes, which are obtained from Eq. (23) by replacing the coefficient differences with the formal errors of the estimated coefficients.

In all cases presented here, the Love number $k_2$ is co-estimated, resulting in uncertainties < 1.3 × 10⁻⁵ (< 0.005% relative to the expected signal), which would allow to distinguish between an ocean-bearing and ocean-free Europa (Wu et al., 2001; Wahr et al., 2006). However, we found that the different scenarios investigated in this section only have a marginal impact on the recovery of $k_2$, so that we will not further discuss this parameter.

### 4.1. Ground tracks repetition

The repetition cycle of an RGTO directly impacts the ground surface coverage, as mentioned in Sec. 2.1. The shorter the repetition cycle, the larger the cycle intertrack (see Table 1), resulting in a lower the spatial resolution of the ground tracks, which can limit the resolution of the estimated gravity field solution.

We derive from Eq. (19) the full wavelength of a degree-n gravity field at the equator as $\lambda_E = \frac{2\pi R_E}{n}$. According to the Nyquist criterion, a sampling of $\frac{\lambda_E}{2}$ allows to recover down to the shortest wavelength of $\frac{2\pi R_E}{n}$. In other words, to avoid aliasing (Sneeuw, 2000) when estimating a degree-n gravity field, the spatial resolution (or cycle intertrack $\delta_{eq}$) of the ground tracks of an m:R RGTO (see Eq. (3)), must be larger than this minimum sampling, i.e., $\delta_{eq} \ge \frac{\lambda_E}{2}$. Hence, we can expect the maximum resolvable degree $n_{max}$ of a recovered gravity field from an m:R RGTO to be given by

$$ n_{max} = \frac{R}{2}. \qquad (24) $$

It can be deduced from this equation that for $R \ge 180$, we can estimate the gravity field coefficients up to d/o 90 without any aliasing. For an orbiter with an altitude between 100 km and 200 km, this corresponds to $m \ge 5$. As an example, for an orbit with an altitude $h$ close to 134 km and $\beta_{Earth} \approx 67^\circ$, the gravity field can be estimated unconstrained up to d/o 62, for every $m \ge 5$ (see $m = 5$ and $m = 26$ in Fig. 7). Indeed, decreasing the repetition rate so that the ground tracks coverage becomes denser does not improve the quality of the gravity field solution.
<figure><img src="image_10.png" /><figcaption>**Fig. 7.** Dependency of recovered gravity field solutions on RGTO parameters. The synthetic gravity field (EURGLMo) is used as a reference. Difference (solid) and error (dashed) degree amplitudes of gravity fields recovered from *m*:R RGTO. For *m* = 1,2, a severe aliasing is visible for the entire spectrum.</figcaption></figure>

On the other hand, when estimating the gravity field up to d/o 90 from an *m*:R RGTO with *R* < 180, one expects an aliasing of the gravity field solution, increasing with decreasing values of *R* (see *m* = 1, 2, 3 in Fig. 7), meaning that some frequencies of the gravity field signal become indistinguishable. For the orbits considered in this study, aliasing is indeed visible as soon as one tries to estimate coefficients of degree larger than 20 × *m*.

A first approach to avoid aliasing would be to solve only for a gravity field solution up to d/o 20 × *m*. However, this necessarily induces an omission error, meaning that the signal from the higher degree coefficients (which are not estimated) leaks into the estimated gravity field coefficients. To prevent this, one can still estimate the gravity field up a higher degree (d/o 90 in our simulation scenario), while also constraining the gravity field using the Kaula law from Eq. (21). In order to

<figure><img src="image_11.png" /><figcaption>**Fig. 8.** Synthetic gravity field (EURGLMo) used as a reference. Difference (solid) and error (dashed) degree amplitudes of several gravity field solutions from a 3 : 118 RGTO. With a perfect a priori knowledge of the gravity field up to d/o 90, the coefficients are estimated up to d/o 90 (red) or 20 × *m* (= 60, green). Using only an a priori gravity field up to d/o 60, the coefficients are estimated up to d/o 60 (blue), or to 90 using a Kaula constraint (magenta).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 9.** Synthetic gravity field (EURGLMo) used as a reference. Difference (solid) and error (dashed) degree amplitudes of gravity field solutions recovered from *m*:R RGTO with *m* ≤ 5, using a Kaula constraint.</figcaption></figure>

validate this approach, we also compute an “artificial” solution, for which the omission error is removed by estimating the gravity field up to d/o 20 × *m*, but assuming a perfect knowledge of the higher degree coefficients (green curve in Fig. 8). Fig. 8 then confirms that the Kaula-constrained gravity field solution tends towards the “artificial” solution.

When using a Kaula constraint, the gravity field can be estimated up to d/o 62 with *m* = 3, 4 to the same precision than with *m* ≥ 5 (Fig. 9). A Kaula constraint also improves the solution for *m* ≤ 2, but the solution is still degraded in the low degrees. A tighter constraint could in principle further improve the results, but it would also bias our solution.

## 4.2. Earth β−angle

The angle $\beta_{Earth}$ between the orbital plane of the probe and the Earth direction plays a significant role in gravity field recovery. When this angle is below a certain threshold $\beta_{Earth,c}$, the probe gets behind Europa as seen from Earth. Then a part of the orbit is not covered by observations collected from stations on Earth (Fig. 10).

<figure><img src="image_13.png" /><figcaption>**Fig. 10.** A 26:1023 RGTO with $\beta_{Earth} = 58^{\circ}$ will be only partially visible from Earth, as shown in this longitudinal zoom of its ground track ($\varphi_c = 45^{\circ}$, $-5^{\circ} \le \lambda \le 5^{\circ}$). The few observations between $+\varphi_c$ and $-\varphi_c$ are still recovered thanks to the rotation of Europa around its rotational axis.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Fig. 7.** Dependency of recovered gravity field solutions on RGTO parameters. The synthetic gravity field (EURGLMo) is used as a reference. Difference (solid) and error (dashed) degree amplitudes of gravity fields recovered from *m*:R RGTO. For *m* = 1,2, a severe aliasing is visible for the entire spectrum.</figcaption></figure>

On the other hand, when estimating the gravity field up to d/o 90 from an *m*:R RGTO with *R* < 180, one expects an aliasing of the gravity field solution, increasing with decreasing values of *R* (see *m* = 1, 2, 3 in Fig. 7), meaning that some frequencies of the gravity field signal become indistinguishable. For the orbits considered in this study, aliasing is indeed visible as soon as one tries to estimate coefficients of degree larger than 20 × *m*.

A first approach to avoid aliasing would be to solve only for a gravity field solution up to d/o 20 × *m*. However, this necessarily induces an omission error, meaning that the signal from the higher degree coefficients (which are not estimated) leaks into the estimated gravity field coefficients. To prevent this, one can still estimate the gravity field up a higher degree (d/o 90 in our simulation scenario), while also constraining the gravity field using the Kaula law from Eq. (21). In order to

<figure><img src="image_12.png" /><figcaption>**Fig. 8.** Synthetic gravity field (EURGLMo) used as a reference. Difference (solid) and error (dashed) degree amplitudes of several gravity field solutions from a 3 : 118 RGTO. With a perfect a priori knowledge of the gravity field up to d/o 90, the coefficients are estimated up to d/o 90 (red) or 20 × *m* (= 60, green). Using only an a priori gravity field up to d/o 60, the coefficients are estimated up to d/o 60 (blue), or to 90 using a Kaula constraint (magenta).</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. 9.** Synthetic gravity field (EURGLMo) used as a reference. Difference (solid) and error (dashed) degree amplitudes of gravity field solutions recovered from *m*:R RGTO with *m* ≤ 5, using a Kaula constraint.</figcaption></figure>

validate this approach, we also compute an “artificial” solution, for which the omission error is removed by estimating the gravity field up to d/o 20 × *m*, but assuming a perfect knowledge of the higher degree coefficients (green curve in Fig. 8). Fig. 8 then confirms that the Kaula-constrained gravity field solution tends towards the “artificial” solution.

When using a Kaula constraint, the gravity field can be estimated up to d/o 62 with *m* = 3, 4 to the same precision than with *m* ≥ 5 (Fig. 9). A Kaula constraint also improves the solution for *m* ≤ 2, but the solution is still degraded in the low degrees. A tighter constraint could in principle further improve the results, but it would also bias our solution.

## 4.2. Earth β−angle

The angle $\beta_{Earth}$ between the orbital plane of the probe and the Earth direction plays a significant role in gravity field recovery. When this angle is below a certain threshold $\beta_{Earth,c}$, the probe gets behind Europa as seen from Earth. Then a part of the orbit is not covered by observations collected from stations on Earth (Fig. 10).

<figure><img src="image_14.png" /><figcaption>**Fig. 10.** A 26:1023 RGTO with $\beta_{Earth} = 58^{\circ}$ will be only partially visible from Earth, as shown in this longitudinal zoom of its ground track ($\varphi_c = 45^{\circ}$, $-5^{\circ} \le \lambda \le 5^{\circ}$). The few observations between $+\varphi_c$ and $-\varphi_c$ are still recovered thanks to the rotation of Europa around its rotational axis.</figcaption></figure>
For a given $\beta_{Earth} < \beta_{Earth,c}$, there exists a latitude band centered around the equator for which the orbit is not covered on the far side of Europa with respect to Earth. Considering a near-circular and quasi-polar orbit, the line of sight from Earth is blocked by Europa if the probe's latitude $\varphi$ fulfills $|\varphi| < \varphi_c$, with

$$ \cos(\varphi_c) \cos(\beta_{Earth}) \approx \sqrt{1 - \left(\frac{R_E}{R_E + h}\right)^2}. \qquad (25) $$

$\beta_{Earth,c}$ can be deduced from this equation with $\varphi_c = 0^{\circ}$, and is a function of the probe's altitude $h$, and the radius of Europa $R_E$ (see Fig. 11).

At worst, the probe cannot be observed from Earth for ~40% of the orbit with a completely edge-on orbit ($\beta_{Earth} = 0^{\circ}$) with an altitude between 100 km and 200 km. This highly affects the visible ground coverage. On the other hand, a completely face-on orbit ($\beta_{Earth} = 90^{\circ}$) will always be visible from Earth, but it is also much less suited to probe the gravity signal, as it will induce only relatively small velocity variations along the line of sight direction sensed by Doppler measurements. For a 89° inclined orbit, $\beta_{Earth}$ is quasi constant during the 3 months mission (see Sec. 2.3). A 26:1023 RGTO does not repeat in 3 months, and with an altitude $h = 133$ km, $\beta_{Earth,c} = 66.8^{\circ}$ (Eq. (25)). We considered several initial positions of the orbital plane, i.e. several different $\beta_{Earth}$ values, to investigate the influence of this parameter.

For $|\beta_{Earth}| > |\beta_{Earth,c}|$, a high $\beta_{Earth}$ degrades the estimation of the low

<figure><img src="image_14.png" /><figcaption>Fig. 11. Critical latitude $\varphi_c$ for which the line of sight is blocked on the far side of Europa with respect to Earth as a function of $\beta_{Earth}$. Larger $\beta_{Earth}$ limit this effect to a narrower set of latitudes around the equator.</figcaption></figure>

order gravity field coefficients (zonal and near zonal) as shown in Fig. 12. On the other hand, for $|\beta_{Earth}| < |\beta_{Earth,c}|$, the higher degrees coefficients are not well determined (see Fig. 13), because of the reduced ground surface coverage.

One could estimate the gravity field to a lower d/o, but this introduces an omission error, as shown in Sec. 4.1. A relatively loose Kaula constraint can be considered instead, to reduce the impact on the lower-degree coefficients. Fig. 14 clearly shows that for all degrees, the gravity field coefficients are better recovered from an edge-on orbit ($\beta_{Earth} = 0^{\circ}$) than from a more face-on orbit ($\beta_{Earth} = 90^{\circ}$).

In conclusion, despite the reduced number of observations, a lower $\beta_{Earth}$ improves the quality of the estimated gravity field using a Kaula constraint, as shown in Fig. 15. On the other hand, the overall quality of the recovered gravity field from an orbit with $\beta_{Earth}$ slightly larger than $\beta_{Earth,c}$ (e.g., $\beta_{Earth} \in [66^{\circ}, 75^{\circ}]$ for $h = 135$ km) is only slightly deteriorated in comparison to the case $\beta_{Earth} = 0^{\circ}$.

In a more realistic case where $\beta_{Earth}$ is not fixed, a combination of low $\beta_{Earth}$ and $\beta_{Earth,c}$ close to $\beta_{Earth,c}$ would be optimal to increase the precision of the recovered gravity field without relying on constraints to mitigate aliasing.

### 4.3. Altitude

We analyzed the impact of the probe's altitude on gravity field recovery for 5:R RGTOs. For $m \ge 5$, the impact of ground track repetition of the studied m:R RGTO is negligible (see Sec. 4.1), i.e., these results are

<figure><img src="image_15.png" /><figcaption>Fig. 13. Difference degree amplitude of unconstrained recovered gravity field solutions for different $\beta_{Earth}$ angles.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 12. Formal errors of the freely estimated gravity field solution recovered from a 26 : 1023 RGTO for several $\beta_{Earth}$ values.</figcaption></figure>
For a given $\beta_{Earth} < \beta_{Earth,c}$, there exists a latitude band centered around the equator for which the orbit is not covered on the far side of Europa with respect to Earth. Considering a near-circular and quasi-polar orbit, the line of sight from Earth is blocked by Europa if the probe's latitude $\varphi$ fulfills $|\varphi| < \varphi_c$, with

$$ \cos(\varphi_c) \cos(\beta_{Earth}) \approx \sqrt{1 - \left(\frac{R_E}{R_E + h}\right)^2}. \quad (25) $$

$\beta_{Earth,c}$ can be deduced from this equation with $\varphi_c = 0^{\circ}$, and is a function of the probe's altitude $h$, and the radius of Europa $R_E$ (see Fig. 11).

At worst, the probe cannot be observed from Earth for ~40% of the orbit with a completely edge-on orbit ($\beta_{Earth} = 0^{\circ}$) with an altitude between 100 km and 200 km. This highly affects the visible ground coverage. On the other hand, a completely face-on orbit ($\beta_{Earth} = 90^{\circ}$) will always be visible from Earth, but it is also much less suited to probe the gravity signal, as it will induce only relatively small velocity variations along the line of sight direction sensed by Doppler measurements. For a 89° inclined orbit, $\beta_{Earth}$ is quasi constant during the 3 months mission (see Sec. 2.3). A 26:1023 RGTO does not repeat in 3 months, and with an altitude $h = 133$ km, $\beta_{Earth,c} = 66.8^{\circ}$ (Eq. (25)). We considered several initial positions of the orbital plane, i.e. several different $\beta_{Earth}$ values, to investigate the influence of this parameter.

For $|\beta_{Earth}| > |\beta_{Earth,c}|$, a high $\beta_{Earth}$ degrades the estimation of the low

<figure><img src="image_15.png" /><figcaption>Fig. 11. Critical latitude $\varphi_c$ for which the line of sight is blocked on the far side of Europa with respect to Earth as a function of $\beta_{Earth}$. Larger $\beta_{Earth}$ limit this effect to a narrower set of latitudes around the equator.</figcaption></figure>

order gravity field coefficients (zonal and near zonal) as shown in Fig. 12. On the other hand, for $|\beta_{Earth}| < |\beta_{Earth,c}|$, the higher degrees coefficients are not well determined (see Fig. 13), because of the reduced ground surface coverage.

One could estimate the gravity field to a lower d/o, but this introduces an omission error, as shown in Sec. 4.1. A relatively loose Kaula constraint can be considered instead, to reduce the impact on the lower-degree coefficients. Fig. 14 clearly shows that for all degrees, the gravity field coefficients are better recovered from an edge-on orbit ($\beta_{Earth} = 0^{\circ}$) than from a more face-on orbit ($\beta_{Earth} = 90^{\circ}$).

In conclusion, despite the reduced number of observations, a lower $\beta_{Earth}$ improves the quality of the estimated gravity field using a Kaula constraint, as shown in Fig. 15. On the other hand, the overall quality of the recovered gravity field from an orbit with $\beta_{Earth}$ slightly larger than $\beta_{Earth,c}$ (e.g., $\beta_{Earth} \in [66^{\circ}, 75^{\circ}]$ for $h = 135$ km) is only slightly deteriorated in comparison to the case $\beta_{Earth} = 0^{\circ}$.

In a more realistic case where $\beta_{Earth}$ is not fixed, a combination of low $\beta_{Earth}$ and $\beta_{Earth,c}$ close to $\beta_{Earth,c}$ would be optimal to increase the precision of the recovered gravity field without relying on constraints to mitigate aliasing.

### 4.3. Altitude

We analyzed the impact of the probe's altitude on gravity field recovery for 5:R RGTOs. For $m \ge 5$, the impact of ground track repetition of the studied m:R RGTO is negligible (see Sec. 4.1), i.e., these results are

<figure><img src="image_16.png" /><figcaption>Fig. 13. Difference degree amplitude of unconstrained recovered gravity field solutions for different $\beta_{Earth}$ angles.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 12. Formal errors of the freely estimated gravity field solution recovered from a 26 : 1023 RGTO for several $\beta_{Earth}$ values.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Fig. 14.** Difference degree amplitude of recovered gravity field solutions with a Kaula constraint ($K = 0.5$) for different $\beta_{Earth}$ angles.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Fig. 15.** Weighted RMS of geoid height differences (in red) and lower latitude covered by the probe on the far side of Europa (in green) as a function of $\beta_{Earth}$. The geoid heights are computed from gravity field solutions up to d/o 60, estimated using a Kaula constraint.</figcaption></figure>

valid for any $m \ge 5$. With an inclination $i = 89^\circ$ we fixed $\beta_{Earth}$ to 69^\circ$, which is larger than $\beta_{Earth,c}$ for $h = 100$ km (see Fig. 11), to avoid any aliasing due to ground coverage issues (see Sec. 4.2) for any orbit altitude. We thus show only unconstrained solutions in this section.

Under these assumptions, the best gravity field solution can be obtained up to d/o 70 from a 5:202 RGTO, with an altitude of $h = 105$ km (Fig. 16). For higher altitudes (lower values of $R$), high degree gravity field coefficients cannot be properly estimated to the same extent. This behaviour is expected and a consequence of the spherical harmonics decomposition given by Eq. (19): the coefficients of degree $n$ are attenuated by the factor $(R_E/R)^n$. As a result, the higher the probe is, the less sensitive it is to the high degrees of the gravity field. For the same reason the impact of a higher altitude orbit is not the same for all degrees, e.g., for $R = 187$, at an altitude of $h = 194$ km, gravity field coefficients can only be estimated up to d/o 50.

Fig. 17 illustrates the global dependency of the gravity field quality with respect to the altitude of the orbit. Because $\beta_{Earth}$ is not precisely 69° for all these tests, the weighted RMS is not rigorously monotonous.

<figure><img src="image_19.png" /><figcaption>**Fig. 16.** Gravity field solutions for orbits with different altitudes with $m = 5$ and $\beta_{Earth}$ close to 69° during the 3 months of the simulated mission period. Synthetic gravity field (EURGLMo) as a reference. Difference (solid) and error (dashed) degree amplitudes.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Fig. 17.** Weighted RMS of geoid height differences between the reference gravity field (EURGLMo) and gravity field solutions for orbits with different altitude with $m = 5$ and $\beta_{Earth}$ close to 69°.</figcaption></figure>

Indeed, the gravity field solution is very sensitive to $\beta_{Earth}$ when being close to $\beta_{Earth,c}$ (see Fig. 15). Nevertheless, the impact of $\beta_{Earth}$ is minor at this scale, in comparison with the impact of the altitude.

One should also note that stronger orbital perturbations from Europa's gravity field would be experienced from very low altitude orbits (e.g. 100 km altitude). The number and magnitude of maneuvers needed to maintain the orbit would be reduced with higher altitude orbits. As a consequence, from a mission point of view, a compromise has to be reached between the resolution of the gravity field desired and the number of maneuvers (which impacts the number of clean observations, and the total $\Delta v$) to ensure a sufficient lifetime of the mission.

### 4.4. Inclination

As mentioned in Sec. 2.4, the choice of inclination impacts the variation of the $\beta_{Earth}$ angle. For a mission in the middle of 2031, a 89°-
<figure><img src="image_18.png" /><figcaption>**Fig. 14.** Difference degree amplitude of recovered gravity field solutions with a Kaula constraint ($K = 0.5$) for different $\beta_{Earth}$ angles.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Fig. 15.** Weighted RMS of geoid height differences (in red) and lower latitude covered by the probe on the far side of Europa (in green) as a function of $\beta_{Earth}$. The geoid heights are computed from gravity field solutions up to d/o 60, estimated using a Kaula constraint.</figcaption></figure>

valid for any $m \ge 5$. With an inclination $i = 89^\circ$ we fixed $\beta_{Earth}$ to 69^\circ$, which is larger than $\beta_{Earth,c}$ for $h = 100$ km (see Fig. 11), to avoid any aliasing due to ground coverage issues (see Sec. 4.2) for any orbit altitude. We thus show only unconstrained solutions in this section.

Under these assumptions, the best gravity field solution can be obtained up to d/o 70 from a 5:202 RGTO, with an altitude of $h = 105$ km (Fig. 16). For higher altitudes (lower values of $R$), high degree gravity field coefficients cannot be properly estimated to the same extent. This behaviour is expected and a consequence of the spherical harmonics decomposition given by Eq. (19): the coefficients of degree $n$ are attenuated by the factor $(R_E/R)^n$. As a result, the higher the probe is, the less sensitive it is to the high degrees of the gravity field. For the same reason the impact of a higher altitude orbit is not the same for all degrees, e.g., for $R = 187$, at an altitude of $h = 194$ km, gravity field coefficients can only be estimated up to d/o 50.

Fig. 17 illustrates the global dependency of the gravity field quality with respect to the altitude of the orbit. Because $\beta_{Earth}$ is not precisely 69° for all these tests, the weighted RMS is not rigorously monotonous.

<figure><img src="image_20.png" /><figcaption>**Fig. 16.** Gravity field solutions for orbits with different altitudes with $m = 5$ and $\beta_{Earth}$ close to 69° during the 3 months of the simulated mission period. Synthetic gravity field (EURGLMo) as a reference. Difference (solid) and error (dashed) degree amplitudes.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Fig. 17.** Weighted RMS of geoid height differences between the reference gravity field (EURGLMo) and gravity field solutions for orbits with different altitude with $m = 5$ and $\beta_{Earth}$ close to 69°.</figcaption></figure>

Indeed, the gravity field solution is very sensitive to $\beta_{Earth}$ when being close to $\beta_{Earth,c}$ (see Fig. 15). Nevertheless, the impact of $\beta_{Earth}$ is minor at this scale, in comparison with the impact of the altitude.

One should also note that stronger orbital perturbations from Europa's gravity field would be experienced from very low altitude orbits (e.g. 100 km altitude). The number and magnitude of maneuvers needed to maintain the orbit would be reduced with higher altitude orbits. As a consequence, from a mission point of view, a compromise has to be reached between the resolution of the gravity field desired and the number of maneuvers (which impacts the number of clean observations, and the total $\Delta v$) to ensure a sufficient lifetime of the mission.

### 4.4. Inclination

As mentioned in Sec. 2.4, the choice of inclination impacts the variation of the $\beta_{Earth}$ angle. For a mission in the middle of 2031, a 89°-
inclined orbit would enable a low variation of $\beta_{Earth}$ angle. This is very practical in the frame of this study, but it is not imperative from a mission point of view. On the contrary, the variation of the $\beta_{Earth}$ angle of a probe on a 80°-inclined orbit would reach up to 60° in 3 months. A combination of several orbital plane configurations with respect to Earth would actually be beneficial to take advantage of both high and low values of $\beta_{Earth}$. As mentioned in Sec. 4.2, high values of $\beta_{Earth}$ enable a denser coverage of Europa, and lower values of $\beta_{Earth}$ would improve the estimation of the lower degree of the gravity field.

Variations of the $\beta_{Earth}$ angle make a proper comparison of gravity field solutions recovered from between a polar and non-polar orbit delicate. We decided to fix the mean value of $\beta_{Earth}$ to 70° for both orbits. $\beta_{Earth}$ variation of the polar orbit is quite small, but in case of the 80° inclined orbit, we set $\beta_{Earth}$ to 47° at orbit injection (Fig. 18).

The consequence of a non-polar orbit is an unobserved gap in the polar regions of the celestial body. The probe in orbit will then be less sensitive to the zonal and near-zonal gravity field coefficients (van Gelderen and Koop, 1997). In terms of gravity field recovery, this means that the estimates of these coefficients will be degraded in comparison with the use of a polar orbit (Fig. 19). The impact on higher order gravity field coefficients is negligible. Some differences can also be seen in the near sectorials of high degrees, but this is a marginal effect due to the slightly lower altitude of the 5:196 non-polar RGTO with respect to the polar 5:197 RGTO.

<figure><img src="image_21.png" /><figcaption>**Fig. 18.** $\beta_{Earth}$ angle variation over 3 months for a polar orbit, starting at $\beta_{Earth,0} = 70^{\circ}$, and for an 80° inclined orbit, starting at $\beta_{Earth,0} = 47^{\circ}$.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**Fig. 19.** Formal errors of the gravity field solution recovered from an *m* = 5 polar RGTO (right) and from an *m* = 5 RGTO with a 10° polar gap.</figcaption></figure>

# 5. Gravity field recovery strategy

Once the optimal scenario has been chosen depending on mission goals, the gravity field can be improved from previous knowledge. In our two step procedure, the arc-wise parameters are first iteratively estimated for every arc, while fixing the global parameters (e.g., gravity field coefficients) to their a priori value. However, the current knowledge of Europa's gravity field is limited to the degree-2 coefficients only (Anderson, 1998). In order for the orbit fit to converge, additional parameters are needed in addition to the initial conditions to compensate for the significant lack of information on the a priori force model. Two solutions are considered in this section. The first one is to use pseudo-stochastic pulses in our estimation, and the other one is to co-estimate in a first iteration the low-degrees coefficients of the gravity field together with the orbital parameters. These solutions are compared with a reference solution, computed using a perfect a priori force model (up to d/o 90), labelled EURSOL01.

## 5.1. Pseudo-stochastic pulses

Pseudo-stochastic pulses are instantaneous velocity changes in pre-defined directions (radial, along-track and cross-track) (Beutler et al., 2010). We estimate a set of pulses every 60 min, while constraining them to $\sigma_p = 1$ mm/s. We use them as additional parameters to absorb modelling errors.

With the help of these pulses, the data fit for all the arcs converged within 8 iterations when starting with a degree-2 a priori gravity field. The level of convergence in terms of Doppler residuals is still very far from the precision one could expect with a better a priori force model (Sec. 3.4), but this is to be expected from the limited knowledge of the gravity field. The orbits are also quite far (tens of kilometers) away from the true orbit. However, after this first step the nominal procedure can be pursued. The global parameters are successfully estimated after stacking the 3 months arc-wise NEQs, pre-eliminating all the arc-wise parameters (orbital elements and pseudo-stochastic pulses) and finally solving for the gravity field coefficients up to d/o 90.

The complete procedure is then iterated using the new gravity field as a priori information. The arc-wise parameters (initial osculating elements and pseudo-stochastic pulses) are not updated from one global solution to the other. They are instead re-estimated for each arc from their initial a priori value. Within a few iterations of the full process, the gravity field can be estimated to the same level of precision as when starting from the reference gravity field EURGLMo. In the first iteration, the RMS values of the pulses in radial, along-track and cross track direction are (35, 46, 1) mm/s. The magnitude of the pulses decreases at each iteration. They reach the level of (1, 2, 0.1) mm/s after a few iterations, until they are not needed anymore to help the first orbit fit to converge. At this point, pulses are not needed anymore and we either apply tight constrains to 0 or fully avoid estimating them in order to preserve the low-degrees of our gravity field solution.

As mentioned in the previous sections, considering an edge-on orbit impacts the visible ground coverage. When considering an orbit for which no observations are lost because of occultations by Europa (i.e. $\beta_{Earth} > \beta_{Earth,c}$), only 3 to 4 iterations are needed to achieve the level of precision one would have when using a perfect a priori gravity field.

However, when considering an edge-on orbit, more iterations can be needed to converge. This is due to the reduced visible ground coverage. As an example, 5 more iterations are needed for a nearly edge-on orbit (see Fig. 20a) than for an orbit with $\beta_{Earth} = \beta_{Earth,c}$. The effect of the reduced observed ground coverage can already be seen in the first iteration. The coefficients with a degree larger than 40 are not well estimated, which is in agreement with the results in Sec. 4.2.

Here again, one could use a loose Kaula constraint to improve the gravity field solution after the first iteration. However, the convergence is still quite slow, even if the gravity field solution improves with respect to a free solution.
inclined orbit would enable a low variation of $\beta_{Earth}$ angle. This is very practical in the frame of this study, but it is not imperative from a mission point of view. On the contrary, the variation of the $\beta_{Earth}$ angle of a probe on a 80°-inclined orbit would reach up to 60° in 3 months. A combination of several orbital plane configurations with respect to Earth would actually be beneficial to take advantage of both high and low values of $\beta_{Earth}$. As mentioned in Sec. 4.2, high values of $\beta_{Earth}$ enable a denser coverage of Europa, and lower values of $\beta_{Earth}$ would improve the estimation of the lower degree of the gravity field.

Variations of the $\beta_{Earth}$ angle make a proper comparison of gravity field solutions recovered from between a polar and non-polar orbit delicate. We decided to fix the mean value of $\beta_{Earth}$ to 70° for both orbits. $\beta_{Earth}$ variation of the polar orbit is quite small, but in case of the 80° inclined orbit, we set $\beta_{Earth}$ to 47° at orbit injection (Fig. 18).

The consequence of a non-polar orbit is an unobserved gap in the polar regions of the celestial body. The probe in orbit will then be less sensitive to the zonal and near-zonal gravity field coefficients (van Gelderen and Koop, 1997). In terms of gravity field recovery, this means that the estimates of these coefficients will be degraded in comparison with the use of a polar orbit (Fig. 19). The impact on higher order gravity field coefficients is negligible. Some differences can also be seen in the near sectorials of high degrees, but this is a marginal effect due to the slightly lower altitude of the 5:196 non-polar RGTO with respect to the polar 5:197 RGTO.

<figure><img src="image_22.png" /><figcaption>**Fig. 18.** $\beta_{Earth}$ angle variation over 3 months for a polar orbit, starting at $\beta_{Earth,0} = 70^{\circ}$, and for an 80° inclined orbit, starting at $\beta_{Earth,0} = 47^{\circ}$.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Fig. 19.** Formal errors of the gravity field solution recovered from an *m* = 5 polar RGTO (right) and from an *m* = 5 RGTO with a 10° polar gap.</figcaption></figure>

# 5. Gravity field recovery strategy

Once the optimal scenario has been chosen depending on mission goals, the gravity field can be improved from previous knowledge. In our two step procedure, the arc-wise parameters are first iteratively estimated for every arc, while fixing the global parameters (e.g., gravity field coefficients) to their a priori value. However, the current knowledge of Europa's gravity field is limited to the degree-2 coefficients only (Anderson, 1998). In order for the orbit fit to converge, additional parameters are needed in addition to the initial conditions to compensate for the significant lack of information on the a priori force model. Two solutions are considered in this section. The first one is to use pseudo-stochastic pulses in our estimation, and the other one is to co-estimate in a first iteration the low-degrees coefficients of the gravity field together with the orbital parameters. These solutions are compared with a reference solution, computed using a perfect a priori force model (up to d/o 90), labelled EURSOL01.

## 5.1. Pseudo-stochastic pulses

Pseudo-stochastic pulses are instantaneous velocity changes in pre-defined directions (radial, along-track and cross-track) (Beutler et al., 2010). We estimate a set of pulses every 60 min, while constraining them to $\sigma_p = 1$ mm/s. We use them as additional parameters to absorb modelling errors.

With the help of these pulses, the data fit for all the arcs converged within 8 iterations when starting with a degree-2 a priori gravity field. The level of convergence in terms of Doppler residuals is still very far from the precision one could expect with a better a priori force model (Sec. 3.4), but this is to be expected from the limited knowledge of the gravity field. The orbits are also quite far (tens of kilometers) away from the true orbit. However, after this first step the nominal procedure can be pursued. The global parameters are successfully estimated after stacking the 3 months arc-wise NEQs, pre-eliminating all the arc-wise parameters (orbital elements and pseudo-stochastic pulses) and finally solving for the gravity field coefficients up to d/o 90.

The complete procedure is then iterated using the new gravity field as a priori information. The arc-wise parameters (initial osculating elements and pseudo-stochastic pulses) are not updated from one global solution to the other. They are instead re-estimated for each arc from their initial a priori value. Within a few iterations of the full process, the gravity field can be estimated to the same level of precision as when starting from the reference gravity field EURGLMo. In the first iteration, the RMS values of the pulses in radial, along-track and cross track direction are (35, 46, 1) mm/s. The magnitude of the pulses decreases at each iteration. They reach the level of (1, 2, 0.1) mm/s after a few iterations, until they are not needed anymore to help the first orbit fit to converge. At this point, pulses are not needed anymore and we either apply tight constrains to 0 or fully avoid estimating them in order to preserve the low-degrees of our gravity field solution.

As mentioned in the previous sections, considering an edge-on orbit impacts the visible ground coverage. When considering an orbit for which no observations are lost because of occultations by Europa (i.e. $\beta_{Earth} > \beta_{Earth,c}$), only 3 to 4 iterations are needed to achieve the level of precision one would have when using a perfect a priori gravity field.

However, when considering an edge-on orbit, more iterations can be needed to converge. This is due to the reduced visible ground coverage. As an example, 5 more iterations are needed for a nearly edge-on orbit (see Fig. 20a) than for an orbit with $\beta_{Earth} = \beta_{Earth,c}$. The effect of the reduced observed ground coverage can already be seen in the first iteration. The coefficients with a degree larger than 40 are not well estimated, which is in agreement with the results in Sec. 4.2.

Here again, one could use a loose Kaula constraint to improve the gravity field solution after the first iteration. However, the convergence is still quite slow, even if the gravity field solution improves with respect to a free solution.
<figure><img src="image_23.png" /><figcaption>**Fig. 20.** Difference degree amplitude of gravity field solutions estimated up to d/o 90 using pulses in the first iterations, and either reintroducing all gravity field parameters at every iteration (a) or reintroducing the gravity field parameters only up to d/o 40 for the first two iterations (b).</figcaption></figure>

An increase in convergence speed can be obtained by estimating gravity field coefficients up to d/o 90 but re-introducing the estimated gravity field solutions only up to d/o 40 as a priori for the following iteration (see Fig. 20b). Doing this for the first 2 iterations, within a total of 5 iterations, the gravity field parameters converge to the same level of precision than when using a perfect a priori gravity field.

One drawback of this approach is that estimating pulses every 60 min during 3 months significantly increases the total number of parameters to be estimated to 15131 parameters (1.73 times more, see Table 3). Here, pulses are only considered as a useful tool to temporarily compensate for force model deficiencies.

## 5.2. Co-estimation of the low-degrees gravity coefficients

Another solution is to estimate only the low-degree gravity field

coefficients along with the orbital elements in one common adjustment. Contrary to our nominal procedure, the orbital elements are thus never estimated alone. Normal equations systems (NEQs) including both orbit and gravity field parameters are set up for each arc but not solved, i.e., we do not compute an intermediate arc-wise orbit-only fit. The arc-wise NEQs are then stacked, the orbit parameters are again pre-eliminated, and the 3-months NEQ is inverted to compute the gravity field solution.

Estimating gravity field coefficients up to d/o 20 results in a total of 8746 parameters, which is approximately half the 17024 parameters estimated when employing pseudo-stochastic pulses, as detailed in Sec. 5.1.

With this method, one can estimate a reasonable medium-degree gravity field solution, in a single iteration, without using pulses (see Fig. 21). After this first iteration, the nominal procedure can be resumed until convergence: a first estimation of local parameters in an arc-wise fit, and then stacking all the NEQs to solve for a global orbit and gravity field solution.

## 5.3. Comparison

The two methods presented in sections 5.1 and 5.2 converge to the same level towards the gravity field reference solution EURSOL01, which is computed using a perfect a priori force model, as shown in the previous sections. Because of the larger number of parameters, the use of pulses increases the total processing time. The differences are presented in Table 4.

In the very first iteration, one can see that pulses are more efficient in reducing the RMS of Doppler residuals (Fig. 22 left). This is to be expected due of the higher number of parameters to absorb model deficiencies. On the other hand, the differences with respect to the true orbits are larger when estimating pulses than when co-estimating gravity field coefficients (Fig. 22 right). Re-injecting only the gravity field up to d/o 40 when using pulses markedly increases the convergence speed.

One can also note that when not estimating pulses anymore (6th iteration for (A) and 4th iteration for (B)) the Doppler residuals can temporarily increase, but then decrease even faster in the next iterations. This is not the case for the orbit differences: the RMS does not increase and the convergence is also faster after the pulses are not estimated anymore.

It is important to emphasize that both methods converge to the same

<figure><img src="image_24.png" /><figcaption>**Fig. 21.** Synthetic gravity field (EURGLMo) as a reference. Difference (solid) and error (dashed) degree amplitudes of gravity field solution recovered from m:R RGTO with $m \ge 5$, without pulses and with an increasing number of estimated gravity field coefficients.</figcaption></figure>
<figure><img src="image_24.png" /><figcaption>**Fig. 20.** Difference degree amplitude of gravity field solutions estimated up to d/o 90 using pulses in the first iterations, and either reintroducing all gravity field parameters at every iteration (a) or reintroducing the gravity field parameters only up to d/o 40 for the first two iterations (b).</figcaption></figure>

An increase in convergence speed can be obtained by estimating gravity field coefficients up to d/o 90 but re-introducing the estimated gravity field solutions only up to d/o 40 as a priori for the following iteration (see Fig. 20b). Doing this for the first 2 iterations, within a total of 5 iterations, the gravity field parameters converge to the same level of precision than when using a perfect a priori gravity field.

One drawback of this approach is that estimating pulses every 60 min during 3 months significantly increases the total number of parameters to be estimated to 15131 parameters (1.73 times more, see Table 3). Here, pulses are only considered as a useful tool to temporarily compensate for force model deficiencies.

## 5.2. Co-estimation of the low-degrees gravity coefficients

Another solution is to estimate only the low-degree gravity field

coefficients along with the orbital elements in one common adjustment. Contrary to our nominal procedure, the orbital elements are thus never estimated alone. Normal equations systems (NEQs) including both orbit and gravity field parameters are set up for each arc but not solved, i.e., we do not compute an intermediate arc-wise orbit-only fit. The arc-wise NEQs are then stacked, the orbit parameters are again pre-eliminated, and the 3-months NEQ is inverted to compute the gravity field solution.

Estimating gravity field coefficients up to d/o 20 results in a total of 8746 parameters, which is approximately half the 17024 parameters estimated when employing pseudo-stochastic pulses, as detailed in Sec. 5.1.

With this method, one can estimate a reasonable medium-degree gravity field solution, in a single iteration, without using pulses (see Fig. 21). After this first iteration, the nominal procedure can be resumed until convergence: a first estimation of local parameters in an arc-wise fit, and then stacking all the NEQs to solve for a global orbit and gravity field solution.

## 5.3. Comparison

The two methods presented in sections 5.1 and 5.2 converge to the same level towards the gravity field reference solution EURSOL01, which is computed using a perfect a priori force model, as shown in the previous sections. Because of the larger number of parameters, the use of pulses increases the total processing time. The differences are presented in Table 4.

In the very first iteration, one can see that pulses are more efficient in reducing the RMS of Doppler residuals (Fig. 22 left). This is to be expected due of the higher number of parameters to absorb model deficiencies. On the other hand, the differences with respect to the true orbits are larger when estimating pulses than when co-estimating gravity field coefficients (Fig. 22 right). Re-injecting only the gravity field up to d/o 40 when using pulses markedly increases the convergence speed.

One can also note that when not estimating pulses anymore (6th iteration for (A) and 4th iteration for (B)) the Doppler residuals can temporarily increase, but then decrease even faster in the next iterations. This is not the case for the orbit differences: the RMS does not increase and the convergence is also faster after the pulses are not estimated anymore.

It is important to emphasize that both methods converge to the same

<figure><img src="image_25.png" /><figcaption>**Fig. 21.** Synthetic gravity field (EURGLMo) as a reference. Difference (solid) and error (dashed) degree amplitudes of gravity field solution recovered from m:R RGTO with $m \ge 5$, without pulses and with an increasing number of estimated gravity field coefficients.</figcaption></figure>
<caption>**Table 4**</caption>

Number of iterations and computation time for all 3 methods for a 3-months d/o 90 gravity field recovery, starting with a d/o 2 a priori gravity field, considering that arc fits can be processed in parallel or not.

<table><thead><tr><th rowspan="2">Method</th><th rowspan="2">Number of iterations</th><th colspan="2">Time</th></tr><tr><th>seq.</th><th>par.</th></tr></thead><tbody><tr><td>Pulses (A)</td><td>11</td><td>66h12</td><td>5h10</td></tr><tr><td>Pulses (B) (remove d/o ≥ 40)</td><td>8</td><td>44h50</td><td>3h33</td></tr><tr><td>Low-degree co-estimation</td><td>8</td><td>39h17</td><td>2h26</td></tr></tbody></table>

<figure><img src="image_25.png" /><figcaption>**Fig. 22.** RMS of Doppler residuals (top) and orbit differences (bottom) with respect to iteration number for different methods. Pulses (A) refers to the use of pulses (up to 5th iteration) and to consistently re-inject the 90-degree gravity field solution for each iteration. Pulses (B) refers to re-injecting only coefficients up to d/o 40 and using pulses until the 3rd iteration. The dashed line represents the RMS reached when using a perfect a priori gravity field.</figcaption></figure>

level of precision toward the reference solution EURSOL01. Convergence seems nevertheless faster when co-estimating the low-degree coefficients in the first iteration.

## 6. Conclusions and discussion

Our study offers an overview of the relative influence of a set of orbit parameters on the quality of gravity field recovery from a low altitude near-polar circular orbit based on Doppler tracking from Earth. In this respect, we developed a toolbox to design Repetitive Ground Tracks Orbits with various orbital characteristics, and we used it to investigate the effect of several orbit configurations and scenarios on the recovery of orbit and geodetic parameters. These orbits allow for a proper characterization of the ground track coverage, and can also be adapted to non-repetitive orbits by choosing a repetition period longer than the mission duration. In this paper, we considered the case of the Galilean moon Europa, but the method could be applied to other celestial bodies in the Solar System. After detailing our simulation framework for the probe's orbit and Doppler observables, we described our approach to the determination of orbit and geodetic parameters, and we discussed the quality of the retrieved orbits. When considering a white noise with an amplitude of 0.10 mm/s at 60 s integration time for the simulated Doppler data, the orbit of a probe orbiting Europa on a polar orbit can be recovered with a precision < 2 m when using a perfect a priori gravity field model. However, when co-estimating gravity field coefficients the orbit quality is slightly degraded, as expected.

We investigated the impact of a set of orbital parameters (ground tracks coverage, $\beta_{Earth}$-angle, altitude, and inclination) on the gravity field recovery by isolating each of them to the extent possible. The results concerning the altitude and inclination confirm our expectations and validate our procedure. A low altitude orbit is beneficial to gravity field recovery, especially for high-degree coefficients. This has to be balanced with the operational preference for higher altitude orbits requiring less orbit maintenance. Also, a non-polar orbit will result in a polar gap with a negative impact on the zonal coefficients, as it is to be expected.

The number of Europa days until which the orbit ground tracks repeat also plays an important role on coverage and resolution. A high number of ground track overlaps can induce large gaps between the ground tracks, which result in a severe aliasing of the gravity field signal when estimating it to a high degree. However, when estimating a d/o 90 European gravity field, any orbit whose ground tracks repeat after 5 Europa days or more, will be equivalent to a non-repetitive orbit.

Furthermore, we analyzed the impact of the angle $\beta_{Earth}$ between the probe's orbital plane and the Earth direction. While an edge-on orbit is beneficial to estimate low degree gravity field coefficients (Doppler tracking is mostly sensitive along the Earth-probe direction), it also limits the fraction of the probe's orbit visible from Earth. This results in a reduced ground coverage and in aliasing at high degrees of the gravity field, which can be partially regularised by, e.g., a weak Kaula rule constraining. A combination between orbits with a low $\beta_{Earth}$ angle and more face-on orbits (allowing for a full coverage, but making Doppler less sensitive to orbital changes) is thus recommended to compute an accurate gravity field solution while minimizing regularization.

A common issue in the outer Solar System is our very limited current knowledge of planetary gravity fields, often limited to the bulk mass or to a few coefficients. This leads to very large Doppler residuals and to potential non-linearities and convergence issues. We compared two different strategies to estimate a full solution for Europa's gravity field when starting from the current a priori knowledge, which is limited to degree 2. One solution is to introduce pseudo-stochastic pulses to compensate for model deficiencies, and to check that their amplitude decreases after each iteration. We also tested co-estimating additional low-degree gravity field coefficients, to gradually reduce Doppler residuals before opening up the full parameter space. Both strategies lead to the same satisfactory accuracy in the recovery of the gravity field used in our simulations. One notable advantage of co-estimating low-degree coefficients is the reduced number of estimated parameters and reaching convergence with less iterations (which results in a reduced computation time).
<caption>**Table 4**</caption>

Number of iterations and computation time for all 3 methods for a 3-months d/o 90 gravity field recovery, starting with a d/o 2 a priori gravity field, considering that arc fits can be processed in parallel or not.

<table><thead><tr><th rowspan="2">Method</th><th rowspan="2">Number of iterations</th><th colspan="2">Time</th></tr><tr><th>seq.</th><th>par.</th></tr></thead><tbody><tr><td>Pulses (A)</td><td>11</td><td>66h12</td><td>5h10</td></tr><tr><td>Pulses (B) (remove d/o ≥ 40)</td><td>8</td><td>44h50</td><td>3h33</td></tr><tr><td>Low-degree co-estimation</td><td>8</td><td>39h17</td><td>2h26</td></tr></tbody></table>

<figure><img src="image_26.png" /><figcaption>**Fig. 22.** RMS of Doppler residuals (top) and orbit differences (bottom) with respect to iteration number for different methods. Pulses (A) refers to the use of pulses (up to 5th iteration) and to consistently re-inject the 90-degree gravity field solution for each iteration. Pulses (B) refers to re-injecting only coefficients up to d/o 40 and using pulses until the 3rd iteration. The dashed line represents the RMS reached when using a perfect a priori gravity field.</figcaption></figure>

level of precision toward the reference solution EURSOL01. Convergence seems nevertheless faster when co-estimating the low-degree coefficients in the first iteration.

## 6. Conclusions and discussion

Our study offers an overview of the relative influence of a set of orbit parameters on the quality of gravity field recovery from a low altitude near-polar circular orbit based on Doppler tracking from Earth. In this respect, we developed a toolbox to design Repetitive Ground Tracks Orbits with various orbital characteristics, and we used it to investigate the effect of several orbit configurations and scenarios on the recovery of orbit and geodetic parameters. These orbits allow for a proper characterization of the ground track coverage, and can also be adapted to non-repetitive orbits by choosing a repetition period longer than the mission duration. In this paper, we considered the case of the Galilean moon Europa, but the method could be applied to other celestial bodies in the Solar System. After detailing our simulation framework for the probe's orbit and Doppler observables, we described our approach to the determination of orbit and geodetic parameters, and we discussed the quality of the retrieved orbits. When considering a white noise with an amplitude of 0.10 mm/s at 60 s integration time for the simulated Doppler data, the orbit of a probe orbiting Europa on a polar orbit can be recovered with a precision < 2 m when using a perfect a priori gravity field model. However, when co-estimating gravity field coefficients the orbit quality is slightly degraded, as expected.

We investigated the impact of a set of orbital parameters (ground tracks coverage, $\beta_{Earth}$-angle, altitude, and inclination) on the gravity field recovery by isolating each of them to the extent possible. The results concerning the altitude and inclination confirm our expectations and validate our procedure. A low altitude orbit is beneficial to gravity field recovery, especially for high-degree coefficients. This has to be balanced with the operational preference for higher altitude orbits requiring less orbit maintenance. Also, a non-polar orbit will result in a polar gap with a negative impact on the zonal coefficients, as it is to be expected.

The number of Europa days until which the orbit ground tracks repeat also plays an important role on coverage and resolution. A high number of ground track overlaps can induce large gaps between the ground tracks, which result in a severe aliasing of the gravity field signal when estimating it to a high degree. However, when estimating a d/o 90 European gravity field, any orbit whose ground tracks repeat after 5 Europa days or more, will be equivalent to a non-repetitive orbit.

Furthermore, we analyzed the impact of the angle $\beta_{Earth}$ between the probe's orbital plane and the Earth direction. While an edge-on orbit is beneficial to estimate low degree gravity field coefficients (Doppler tracking is mostly sensitive along the Earth-probe direction), it also limits the fraction of the probe's orbit visible from Earth. This results in a reduced ground coverage and in aliasing at high degrees of the gravity field, which can be partially regularised by, e.g., a weak Kaula rule constraining. A combination between orbits with a low $\beta_{Earth}$ angle and more face-on orbits (allowing for a full coverage, but making Doppler less sensitive to orbital changes) is thus recommended to compute an accurate gravity field solution while minimizing regularization.

A common issue in the outer Solar System is our very limited current knowledge of planetary gravity fields, often limited to the bulk mass or to a few coefficients. This leads to very large Doppler residuals and to potential non-linearities and convergence issues. We compared two different strategies to estimate a full solution for Europa's gravity field when starting from the current a priori knowledge, which is limited to degree 2. One solution is to introduce pseudo-stochastic pulses to compensate for model deficiencies, and to check that their amplitude decreases after each iteration. We also tested co-estimating additional low-degree gravity field coefficients, to gradually reduce Doppler residuals before opening up the full parameter space. Both strategies lead to the same satisfactory accuracy in the recovery of the gravity field used in our simulations. One notable advantage of co-estimating low-degree coefficients is the reduced number of estimated parameters and reaching convergence with less iterations (which results in a reduced computation time).
Our results illustrate how decisive the orbit choice can be for the recovery of the gravity field of a celestial body. We offer operative indications to inform preliminary mission design and for the subsequent mission specific analysis, usually a compromise between multiple science objectives and mission constraints. Some aspects were neglected in this study, such as non-gravitational accelerations, or a detailed observation schedule of the ground stations, which could be investigated by future works. Additional measurements could further improve orbit determination such as those provided by an on-board accelerometer or by a laser altimeter.

## Author contributions

**William Desprats:** Conceptualization, Methodology, Software, Writing - Original Draft, Investigation.

**Stefano Bertone:** Conceptualization, Software, Writing - Original Draft, Writing - Review & Editing.

**Daniel Arnold:** Conceptualization, Writing - Review & Editing, Supervision, Funding acquisition.

**Adrian Jäggı:** Conceptualization, Writing - Review & Editing, Supervision.

**Michel Blanc:** Writing - Review & Editing, Funding acquisition.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Data availability

Data will be made available on request.

## Acknowledgments

Calculations were performed on UBELIX, the HPC cluster at the University of Bern (http://www.id.unibe.ch/hpc).

This study has been funded with the support of the Swiss National Science Foundation (SNSF) grant #200021_185056 "Callisto geodesy: A simulation study to support further space missions to the Jovian system".

Stefano Bertone acknowledges support by NASA under award number 80GSFC21M0002.

## Appendix A. Differential correction

Our differential correction procedure is adapted from Russell (2006) and Pavlak (2013). The equations of motion of the Hill model provided by Eq. (5) are invariant under the transformation $t \to -t$, $y \to -y$, $z \to -z$. Indeed when applying this transformation, we obtain

$$ \mathbf{x}_1 = (x_1, y_1, z_1) = (x_0, -y_0, -z_0) $$

$$ \dot{\mathbf{x}}_I = (\dot{x}_1, \dot{y}_1, \dot{z}_1) = (-\dot{x}_0, \dot{y}_0, \dot{z}_0) $$

$$ \ddot{\mathbf{x}}_I = (\ddot{x}_1, \ddot{y}_1, \ddot{z}_1) = (\ddot{x}_0, -\ddot{y}_0, -\ddot{z}_0) $$

It means that the trajectory of a probe starting at $x(t_0) = (x_0, y_0, z_0)$ is a 180° rotation around the x-axis of a backward integrated trajectory of a probe starting at $(x_0, -y_0, -z_0)$. This property is called axi-symmetry (around the x-axis).

Because of this symmetry, if an orbit starts on the x-axis, and if the orbit is perpendicular to the x-axis, i.e.,

$$ \mathbf{x}(t_0) = (x_0, 0, 0) \qquad (A.1) $$

$$ \dot{\mathbf{x}}(t_0) = (0, \dot{y}_0, \dot{z}_0), \qquad (A.2) $$

the resulting propagated orbit and its image (180° rotation around the x-axis) will form one continuous trajectory. Enforcing the condition

$$ y_T = z_T = \dot{x}_T = 0 \qquad (A.2) $$

ensures that the trajectory crosses perpendicularly the x-axis at $t = T$. Then, the trajectory and its image will result in an orbit which closes after 2T.

We consequently decided to look for such axi-symmetric m:R Repetitive Ground Tracks Orbits, with the initial conditions (A.1). In order to find these orbits, we perform a differential correction, which we then derive by means of the Newton-Raphson method. Let us define the free variable Y as

$$ \mathbf{Y} = \begin{pmatrix} x_0 \\ \dot{y}_0 \\ \dot{z}_0 \\ T \end{pmatrix}, \qquad (A.3) $$

and the function F(Y) as

$$ \mathbf{F}(\mathbf{Y}) = \begin{pmatrix} y_T \\ z_T \\ \dot{x}_T \end{pmatrix}, \qquad (A.4) $$

where $y_T$, $z_T$, $\dot{x}_T$ are obtained by propagating (A.1) to $t = T$. The constraint (A.2) can then be enforced by $\mathbf{F}(\mathbf{Y}) = 0$, which can be solved iteratively, so that for every iteration k we get

$$ \mathbf{F}(\mathbf{Y}^k) = -\mathbf{D}\mathbf{F}(\mathbf{Y}^k) \cdot \delta\mathbf{Y}^k, \qquad (A.5) $$
Our results illustrate how decisive the orbit choice can be for the recovery of the gravity field of a celestial body. We offer operative indications to inform preliminary mission design and for the subsequent mission specific analysis, usually a compromise between multiple science objectives and mission constraints. Some aspects were neglected in this study, such as non-gravitational accelerations, or a detailed observation schedule of the ground stations, which could be investigated by future works. Additional measurements could further improve orbit determination such as those provided by an on-board accelerometer or by a laser altimeter.

## Author contributions

**William Desprats:** Conceptualization, Methodology, Software, Writing - Original Draft, Investigation.

**Stefano Bertone:** Conceptualization, Software, Writing - Original Draft, Writing - Review & Editing.

**Daniel Arnold:** Conceptualization, Writing - Review & Editing, Supervision, Funding acquisition.

**Adrian Jäggı:** Conceptualization, Writing - Review & Editing, Supervision.

**Michel Blanc:** Writing - Review & Editing, Funding acquisition.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Data availability

Data will be made available on request.

## Acknowledgments

Calculations were performed on UBELIX, the HPC cluster at the University of Bern (http://www.id.unibe.ch/hpc).

This study has been funded with the support of the Swiss National Science Foundation (SNSF) grant #200021_185056 "Callisto geodesy: A simulation study to support further space missions to the Jovian system".

Stefano Bertone acknowledges support by NASA under award number 80GSFC21M0002.

## Appendix A. Differential correction

Our differential correction procedure is adapted from Russell (2006) and Pavlak (2013). The equations of motion of the Hill model provided by Eq. (5) are invariant under the transformation $t \to -t$, $y \to -y$, $z \to -z$. Indeed when applying this transformation, we obtain

$$ \mathbf{x}_1 = (x_1, y_1, z_1) = (x_0, -y_0, -z_0) $$

$$ \dot{\mathbf{x}}_I = (\dot{x}_1, \dot{y}_1, \dot{z}_1) = (-\dot{x}_0, \dot{y}_0, \dot{z}_0) $$

$$ \ddot{\mathbf{x}}_I = (\ddot{x}_1, \ddot{y}_1, \ddot{z}_1) = (\ddot{x}_0, -\ddot{y}_0, -\ddot{z}_0) $$

It means that the trajectory of a probe starting at $x(t_0) = (x_0, y_0, z_0)$ is a 180° rotation around the x-axis of a backward integrated trajectory of a probe starting at $(x_0, -y_0, -z_0)$. This property is called axi-symmetry (around the x-axis).

Because of this symmetry, if an orbit starts on the x-axis, and if the orbit is perpendicular to the x-axis, i.e.,

$$ \begin{aligned} \mathbf{x}(t_0) &= (x_0, 0, 0) \\ \dot{\mathbf{x}}(t_0) &= (0, \dot{y}_0, \dot{z}_0), \end{aligned} \qquad (A.1) $$

the resulting propagated orbit and its image (180° rotation around the x-axis) will form one continuous trajectory. Enforcing the condition

$$ y_T = z_T = \dot{x}_T = 0 \qquad (A.2) $$

ensures that the trajectory crosses perpendicularly the x-axis at $t = T$. Then, the trajectory and its image will result in an orbit which closes after 2T.

We consequently decided to look for such axi-symmetric m:R Repetitive Ground Tracks Orbits, with the initial conditions (A.1). In order to find these orbits, we perform a differential correction, which we then derive by means of the Newton-Raphson method. Let us define the free variable Y as

$$ \mathbf{Y} = \begin{pmatrix} x_0 \\ \dot{y}_0 \\ \dot{z}_0 \\ T \end{pmatrix}, \qquad (A.3) $$

and the function F(Y) as

$$ \mathbf{F}(\mathbf{Y}) = \begin{pmatrix} y_T \\ z_T \\ \dot{x}_T \end{pmatrix}, \qquad (A.4) $$

where $y_T$, $z_T$, $\dot{x}_T$ are obtained by propagating (A.1) to $t = T$. The constraint (A.2) can then be enforced by $\mathbf{F}(\mathbf{Y}) = 0$, which can be solved iteratively, so that for every iteration k we get

$$ \mathbf{F}(\mathbf{Y}^k) = -\mathbf{D}\mathbf{F}(\mathbf{Y}^k) \cdot \delta\mathbf{Y}^k, \qquad (A.5) $$
with

$$ \delta \mathbf{Y} = \begin{pmatrix} \delta x_0 \\ \delta \dot{y}_0 \\ \delta \dot{z}_0 \\ \delta T \end{pmatrix} \qquad (A.6) $$

and

$$ \begin{aligned} \mathbf{DF}(\mathbf{Y}) &= \frac{\partial \mathbf{F}}{\partial \mathbf{Y}} \\ &= \begin{pmatrix} \frac{\partial y}{\partial x_0} & \frac{\partial y}{\partial \dot{y}_0} & \frac{\partial y}{\partial \dot{z}_0} & \frac{\partial y}{\partial T} \\ \frac{\partial z}{\partial x_0} & \frac{\partial z}{\partial \dot{y}_0} & \frac{\partial z}{\partial \dot{z}_0} & \frac{\partial z}{\partial T} \\ \frac{\partial \dot{x}}{\partial x_0} & \frac{\partial \dot{x}}{\partial \dot{y}_0} & \frac{\partial \dot{x}}{\partial \dot{z}_0} & \frac{\partial \dot{x}}{\partial T} \end{pmatrix} \\ &= \begin{pmatrix} \mathbf{\Phi}_{21} & \mathbf{\Phi}_{25} & \mathbf{\Phi}_{26} & \dot{y}_T \\ \mathbf{\Phi}_{31} & \mathbf{\Phi}_{35} & \mathbf{\Phi}_{36} & \dot{z}_T \\ \mathbf{\Phi}_{41} & \mathbf{\Phi}_{45} & \mathbf{\Phi}_{46} & \ddot{x}_T \end{pmatrix}, \end{aligned} \qquad (A.7) $$

with $\mathbf{\Phi}(t, t_0)$ the state transition matrix (STM), which maps an initial perturbation of $(x(t_0), \dot{x}(t_0))$ at $t_0$ to the resulting perturbation at $t$, such that

$$ \mathbf{\Phi}(t, t_0) = \begin{pmatrix} \frac{\partial x}{\partial x_0} & \frac{\partial x}{\partial \dot{y}_0} & \frac{\partial x}{\partial \dot{z}_0} & \frac{\partial x}{\partial \dot{x}_0} & \frac{\partial x}{\partial \dot{y}_0} & \frac{\partial x}{\partial \dot{z}_0} \\ \frac{\partial y}{\partial x_0} & \frac{\partial y}{\partial \dot{y}_0} & \frac{\partial y}{\partial \dot{z}_0} & \frac{\partial y}{\partial \dot{x}_0} & \frac{\partial y}{\partial \dot{y}_0} & \frac{\partial y}{\partial \dot{z}_0} \\ \frac{\partial z}{\partial x_0} & \frac{\partial z}{\partial \dot{y}_0} & \frac{\partial z}{\partial \dot{z}_0} & \frac{\partial z}{\partial \dot{x}_0} & \frac{\partial z}{\partial \dot{y}_0} & \frac{\partial z}{\partial \dot{z}_0} \\ \frac{\partial \dot{x}}{\partial x_0} & \frac{\partial \dot{x}}{\partial \dot{y}_0} & \frac{\partial \dot{x}}{\partial \dot{z}_0} & \frac{\partial \dot{x}}{\partial \dot{x}_0} & \frac{\partial \dot{x}}{\partial \dot{y}_0} & \frac{\partial \dot{x}}{\partial \dot{z}_0} \\ \frac{\partial \dot{y}}{\partial x_0} & \frac{\partial \dot{y}}{\partial \dot{y}_0} & \frac{\partial \dot{y}}{\partial \dot{z}_0} & \frac{\partial \dot{y}}{\partial \dot{x}_0} & \frac{\partial \dot{y}}{\partial \dot{y}_0} & \frac{\partial \dot{y}}{\partial \dot{z}_0} \\ \frac{\partial \dot{z}}{\partial x_0} & \frac{\partial \dot{z}}{\partial \dot{y}_0} & \frac{\partial \dot{z}}{\partial \dot{z}_0} & \frac{\partial \dot{z}}{\partial \dot{x}_0} & \frac{\partial \dot{z}}{\partial \dot{y}_0} & \frac{\partial \dot{z}}{\partial \dot{z}_0} \end{pmatrix} $$

The latter is integrated along the orbit with the initial condition $\mathbf{\Phi}(t_0, t_0) = \mathbf{I}_6$. Equation (A.5) is under-determined. We decided to add a constraint by fixing the inclination $i_0$ at $t_0$. Using Eq. (A.1) the inclination can be written as

$$ \tan(i_0) = \frac{\dot{z}_0}{\dot{y}_0 + n_J x_0}, \qquad (A.8) $$

where $n_J$ is the mean motion of Europa around Jupiter. We can already apply this constraint to Eq. (A.5), but we rather first reduce the problem. We fix the half-period $T$ by numerically integrating the trajectory until the orbit crosses $R$ times the equatorial plane at $t = T$, enforcing then the constraint $z_T = 0$. Introducing $z_T = 0$ into Eq. (A.5) and considering (A.6) and (A.7) then yields

$$ z_T = \mathbf{\Phi}_{31} \delta x + \mathbf{\Phi}_{35} \delta \dot{y}_0 + \mathbf{\Phi}_{36} \delta \dot{z}_0 + \dot{z}_T \delta T \stackrel{!}{=} 0, $$

so that

$$ \delta T = -\frac{1}{\dot{z}_T} (\mathbf{\Phi}_{31} \delta x + \mathbf{\Phi}_{35} \delta \dot{y}_0 + \mathbf{\Phi}_{36} \delta \dot{z}_0). $$

We can rewrite Eq. (A.5), omitting the k indices, with a reduced function

$$ \tilde{\mathbf{F}}(\mathbf{Y}) = \begin{pmatrix} y_T \\ \dot{x}_T \end{pmatrix} = - \left( \begin{pmatrix} \mathbf{\Phi}_{21} & \mathbf{\Phi}_{25} & \mathbf{\Phi}_{26} \\ \mathbf{\Phi}_{41} & \mathbf{\Phi}_{45} & \mathbf{\Phi}_{46} \end{pmatrix} - \frac{1}{\dot{z}_T} \begin{pmatrix} \dot{y}_T \\ \ddot{x}_T \end{pmatrix} \begin{pmatrix} \mathbf{\Phi}_{31} & \mathbf{\Phi}_{35} & \mathbf{\Phi}_{36} \end{pmatrix} \right) \begin{pmatrix} \delta x_0 \\ \delta \dot{y}_0 \\ \delta \dot{z}_0 \end{pmatrix}. \qquad (A.9) $$
with

$$ \delta \mathbf{Y} = \begin{pmatrix} \delta x_0 \\ \delta \dot{y}_0 \\ \delta \dot{z}_0 \\ \delta T \end{pmatrix} \qquad (A.6) $$

and

$$ \begin{aligned} \mathbf{DF}(\mathbf{Y}) &= \frac{\partial \mathbf{F}}{\partial \mathbf{Y}} \\ &= \begin{pmatrix} \frac{\partial y}{\partial x_0} & \frac{\partial y}{\partial \dot{y}_0} & \frac{\partial y}{\partial \dot{z}_0} & \frac{\partial y}{\partial T} \\ \frac{\partial z}{\partial x_0} & \frac{\partial z}{\partial \dot{y}_0} & \frac{\partial z}{\partial \dot{z}_0} & \frac{\partial z}{\partial T} \\ \frac{\partial \dot{x}}{\partial x_0} & \frac{\partial \dot{x}}{\partial \dot{y}_0} & \frac{\partial \dot{x}}{\partial \dot{z}_0} & \frac{\partial \dot{x}}{\partial T} \end{pmatrix} \\ &= \begin{pmatrix} \Phi_{21} & \Phi_{25} & \Phi_{26} & \dot{y}_T \\ \Phi_{31} & \Phi_{35} & \Phi_{36} & \dot{z}_T \\ \Phi_{41} & \Phi_{45} & \Phi_{46} & \ddot{x}_T \end{pmatrix}, \end{aligned} \qquad (A.7) $$

with $\Phi(t, t_0)$ the state transition matrix (STM), which maps an initial perturbation of $(x(t_0), \dot{x}(t_0))$ at $t_0$ to the resulting perturbation at $t$, such that

$$ \mathbf{\Phi}(t, t_0) = \begin{pmatrix} \frac{\partial x}{\partial x_0} & \frac{\partial x}{\partial \dot{y}_0} & \frac{\partial x}{\partial \dot{z}_0} & \frac{\partial x}{\partial \dot{x}_0} & \frac{\partial x}{\partial \dot{y}_0} & \frac{\partial x}{\partial \dot{z}_0} \\ \frac{\partial y}{\partial x_0} & \frac{\partial y}{\partial \dot{y}_0} & \frac{\partial y}{\partial \dot{z}_0} & \frac{\partial y}{\partial \dot{x}_0} & \frac{\partial y}{\partial \dot{y}_0} & \frac{\partial y}{\partial \dot{z}_0} \\ \frac{\partial z}{\partial x_0} & \frac{\partial z}{\partial \dot{y}_0} & \frac{\partial z}{\partial \dot{z}_0} & \frac{\partial z}{\partial \dot{x}_0} & \frac{\partial z}{\partial \dot{y}_0} & \frac{\partial z}{\partial \dot{z}_0} \\ \frac{\partial \dot{x}}{\partial x_0} & \frac{\partial \dot{x}}{\partial \dot{y}_0} & \frac{\partial \dot{x}}{\partial \dot{z}_0} & \frac{\partial \dot{x}}{\partial \dot{x}_0} & \frac{\partial \dot{x}}{\partial \dot{y}_0} & \frac{\partial \dot{x}}{\partial \dot{z}_0} \\ \frac{\partial \dot{y}}{\partial x_0} & \frac{\partial \dot{y}}{\partial \dot{y}_0} & \frac{\partial \dot{y}}{\partial \dot{z}_0} & \frac{\partial \dot{y}}{\partial \dot{x}_0} & \frac{\partial \dot{y}}{\partial \dot{y}_0} & \frac{\partial \dot{y}}{\partial \dot{z}_0} \\ \frac{\partial \dot{z}}{\partial x_0} & \frac{\partial \dot{z}}{\partial \dot{y}_0} & \frac{\partial \dot{z}}{\partial \dot{z}_0} & \frac{\partial \dot{z}}{\partial \dot{x}_0} & \frac{\partial \dot{z}}{\partial \dot{y}_0} & \frac{\partial \dot{z}}{\partial \dot{z}_0} \end{pmatrix} $$

The latter is integrated along the orbit with the initial condition $\Phi(t_0, t_0) = \mathbf{I}_6$. Equation (A.5) is under-determined. We decided to add a constraint by fixing the inclination $i_0$ at $t_0$. Using Eq. (A.1) the inclination can be written as

$$ \tan(i_0) = \frac{\dot{z}_0}{\dot{y}_0 + n_J x_0}, \qquad (A.8) $$

where $n_J$ is the mean motion of Europa around Jupiter. We can already apply this constraint to Eq. (A.5), but we rather first reduce the problem. We fix the half-period $T$ by numerically integrating the trajectory until the orbit crosses $R$ times the equatorial plane at $t = T$, enforcing then the constraint $z_T = 0$. Introducing $z_T = 0$ into Eq. (A.5) and considering (A.6) and (A.7) then yields

$$ z_T = \Phi_{31} \delta x + \Phi_{35} \delta \dot{y}_0 + \Phi_{36} \delta \dot{z}_0 + \dot{z}_T \delta T \stackrel{!}{=} 0, $$

so that

$$ \delta T = -\frac{1}{\dot{z}_T} (\Phi_{31} \delta x + \Phi_{35} \delta \dot{y}_0 + \Phi_{36} \delta \dot{z}_0). $$

We can rewrite Eq. (A.5), omitting the k indices, with a reduced function

$$ \tilde{\mathbf{F}}(\mathbf{Y}) = \begin{pmatrix} y_T \\ \dot{x}_T \end{pmatrix} = - \left( \begin{pmatrix} \Phi_{21} & \Phi_{25} & \Phi_{26} \\ \Phi_{41} & \Phi_{45} & \Phi_{46} \end{pmatrix} - \frac{1}{\dot{z}_T} \begin{pmatrix} \dot{y}_T \\ \ddot{x}_T \end{pmatrix} \begin{pmatrix} \Phi_{31} & \Phi_{35} & \Phi_{36} \end{pmatrix} \right) \begin{pmatrix} \delta x_0 \\ \delta \dot{y}_0 \\ \delta \dot{z}_0 \end{pmatrix}. \qquad (A.9) $$
We then reduce the number of free variables by using Eq. (A.8) to compute $\dot{y}_0$ such that

$$ \dot{y}_0 = \dot{z}_0 \cot(i_0) - n_J x_0. \qquad (A.10) $$

By introducing Eq. (A.10) into Eq. (A.9), we finally get

$$ \tilde{\mathbf{F}}(\mathbf{Y}) = -\mathbf{A} \cdot \delta \tilde{\mathbf{Y}}, \qquad (A.11) $$

with

$$ \mathbf{A} = \begin{pmatrix} \Phi_{21} & \Phi_{26} \\ \Phi_{41} & \Phi_{46} \end{pmatrix} + \begin{pmatrix} \Phi_{25} \\ \Phi_{45} \end{pmatrix} (-n_J \cot(i_0)) - \frac{1}{\dot{z}_T} \begin{pmatrix} \dot{y}_T \\ \ddot{x}_T \end{pmatrix} \left( \begin{pmatrix} \Phi_{31} & \Phi_{36} \\ \Phi_{35} & \Phi_{36} \end{pmatrix} + \Phi_{35} (-n_J \cot(i_0)) \right), $$

and the correction to the reduced free variable vector $\tilde{\mathbf{Y}}$

$$ \delta \tilde{\mathbf{Y}} = \begin{pmatrix} \delta x_0 \\ \delta \dot{z}_0 \end{pmatrix}. $$

To summarize, at every iteration the orbit and the STM are propagated until the orbit crosses the equatorial plane R times, defining the half-period T. The resulting $\Phi(T, t_0)$, $\dot{x}(T)$ and $\ddot{x}(T)$ are used to compute corrections on $x_0$ and $\dot{z}_0$ based on Eq. (A.11), while $\dot{y}_0$ is recomputed using Eq. (A.10). This procedure is repeated until convergence to the desired level.

## Appendix B. Evolution of the orbital elements

In this appendix, we derive the evolution over time of the inclination $i$ and of the longitude of the ascending node $\Omega$. We use the Gaussian perturbation equations (Beutler, 2005) to express their time derivatives as a function of the forces perturbing the Keplerian motion. In our case, the perturbing forces are the third-body attraction of Jupiter, and the force induced by the degree-2 gravity field of Europa. The derivatives of both the inclination and the longitude of the ascending node only depend on the cross-track direction of these perturbing forces, such that

$$ \begin{aligned} \frac{di}{dt} &= \frac{r \cos i}{h} (N_J + N_P) \\ \frac{d\Omega}{dt} &= \frac{r \sin u}{h \sin i} (N_J + N_P), \end{aligned} \qquad (B.1) $$

with

$$ N_J = \nabla_x (U_J)_{\text{loc}} \cdot \mathbf{e}_N $$

$$ N_P = \nabla_x (U_P)_{\text{loc}} \cdot \mathbf{e}_N, $$

where loc indicates the “local orbital frame” (i.e., Radial, Tangential, Normal). The potential $U_J$ and $U_P$ in the rotating frame rot are then given by Eq. (7) and by Eq. (8) rewritten as

$$ \begin{aligned} U_J &= \frac{3x^2 - r^2}{2} n_J^2 \\ U_P &= \frac{\mu_E R_E^2 J_2}{r^2} \frac{g(\mathbf{x}_{\text{rot}})}{r^2}, \end{aligned} \qquad (B.2) $$

with

$$ \mathbf{x}_{\text{rot}} = (x, y, z) $$

$$ g(\mathbf{x}_{\text{rot}}) = 7x^2 - 2y^2 - 5z^2. $$

## Appendix B.1. Relation between the rotating frame, the inertial frame and the local orbital frame

Let us neglect the obliquity of Europa (< 0.1°) and the inclination of Europa's orbital plane with respect to the ecliptic (< 2°). The rotating frame is then obtained by a rotation around the z-axis from the inertial frame. Assuming that at epoch t = 0, Europa and Jupiter lie on the x-axis of the inertial frame, we can write

$$ \mathbf{x}_{\text{rot}} = \mathbf{R}_3(n_J t) \cdot \mathbf{x}_{\text{in}}, $$

where $\mathbf{x}_{\text{rot}}$ and $\mathbf{x}_{\text{in}}$ are the components of the position vectors in the rotating frame and in the inertial frame, and $\mathbf{R}_i(\theta)$ is the 3 × 3 matrix representing a
We then reduce the number of free variables by using Eq. (A.8) to compute $\dot{y}_0$ such that

$$ \dot{y}_0 = \dot{z}_0 \cot(i_0) - n_J x_0. \qquad (A.10) $$

By introducing Eq. (A.10) into Eq. (A.9), we finally get

$$ \tilde{\mathbf{F}}(\mathbf{Y}) = -\mathbf{A} \cdot \delta \tilde{\mathbf{Y}}, \qquad (A.11) $$

with

$$ \mathbf{A} = \begin{pmatrix} \Phi_{21} & \Phi_{26} \\ \Phi_{41} & \Phi_{46} \end{pmatrix} + \begin{pmatrix} \Phi_{25} \\ \Phi_{45} \end{pmatrix} (-n_J \cot(i_0)) - \frac{1}{\dot{z}_T} \begin{pmatrix} \dot{y}_T \\ \ddot{x}_T \end{pmatrix} \left( \begin{pmatrix} \Phi_{31} & \Phi_{36} \\ \Phi_{35} & \Phi_{36} \end{pmatrix} + \Phi_{35} (-n_J \cot(i_0)) \right), $$

and the correction to the reduced free variable vector $\tilde{\mathbf{Y}}$

$$ \delta \tilde{\mathbf{Y}} = \begin{pmatrix} \delta x_0 \\ \delta \dot{z}_0 \end{pmatrix}. $$

To summarize, at every iteration the orbit and the STM are propagated until the orbit crosses the equatorial plane R times, defining the half-period T. The resulting $\Phi(T, t_0)$, $\dot{x}(T)$ and $\ddot{x}(T)$ are used to compute corrections on $x_0$ and $\dot{z}_0$ based on Eq. (A.11), while $\dot{y}_0$ is recomputed using Eq. (A.10). This procedure is repeated until convergence to the desired level.

## Appendix B. Evolution of the orbital elements

In this appendix, we derive the evolution over time of the inclination $i$ and of the longitude of the ascending node $\Omega$. We use the Gaussian perturbation equations (Beutler, 2005) to express their time derivatives as a function of the forces perturbing the Keplerian motion. In our case, the perturbing forces are the third-body attraction of Jupiter, and the force induced by the degree-2 gravity field of Europa. The derivatives of both the inclination and the longitude of the ascending node only depend on the cross-track direction of these perturbing forces, such that

$$ \begin{aligned} \frac{di}{dt} &= \frac{r \cos i}{h} (N_J + N_P) \\ \frac{d\Omega}{dt} &= \frac{r \sin u}{h \sin i} (N_J + N_P), \end{aligned} \qquad (B.1) $$

with

$$ N_J = \nabla_x (U_J)_{\text{loc}} \cdot e_N $$

$$ N_P = \nabla_x (U_P)_{\text{loc}} \cdot e_N, $$

where loc indicates the “local orbital frame” (i.e., Radial, Tangential, Normal). The potential $U_J$ and $U_P$ in the rotating frame rot are then given by Eq. (7) and by Eq. (8) rewritten as

$$ \begin{aligned} U_J &= \frac{3x^2 - r^2}{2} n_J^2 \\ U_P &= \frac{\mu_E R_E^2 J_2}{r^2} \frac{g(\mathbf{x}_{\text{rot}})}{r^2}, \end{aligned} \qquad (B.2) $$

with

$$ \mathbf{x}_{\text{rot}} = (x, y, z) $$

$$ g(\mathbf{x}_{\text{rot}}) = 7x^2 - 2y^2 - 5z^2. $$

## Appendix B.1. Relation between the rotating frame, the inertial frame and the local orbital frame

Let us neglect the obliquity of Europa (< 0.1°) and the inclination of Europa's orbital plane with respect to the ecliptic (< 2°). The rotating frame is then obtained by a rotation around the z-axis from the inertial frame. Assuming that at epoch t = 0, Europa and Jupiter lie on the x-axis of the inertial frame, we can write

$$ \mathbf{x}_{\text{rot}} = \mathbf{R}_3(n_J t) \cdot \mathbf{x}_{\text{in}}, $$

where $\mathbf{x}_{\text{rot}}$ and $\mathbf{x}_{\text{in}}$ are the components of the position vectors in the rotating frame and in the inertial frame, and $\mathbf{R}_i(\theta)$ is the 3 × 3 matrix representing a
rotation about the i-th coordinate axis by an angle θ.

Additionally, the transformation between the inertial frame and the local orbital frame can be written as

$$ \mathbf{x}_{\text{loc}} = \mathbf{R}_3(u) \cdot \mathbf{R}_1(i) \cdot \mathbf{R}_3(\Omega) \cdot \mathbf{x}_{\text{in}}, $$

which yields

$$ \mathbf{x}_{\text{loc}} = \mathbf{R}_{\text{tot}}(u, i, u_J) \cdot \mathbf{x}_{\text{rot}} \qquad (B.3) $$

with $\mathbf{R}_{\text{tot}}(u, i, u_J) = \mathbf{R}_3(u) \cdot \mathbf{R}_1(i) \cdot \mathbf{R}_3(-u_J)$ and $u_J = n_J t - \Omega$.

Once the gradients of the potentials are computed in the rotating frame, they need to be converted to the local orbital frame, so that

$$ \nabla_x(U_J + U_P)_{\text{loc}} = \mathbf{R}_{\text{tot}}(u, i, u_J) \cdot \nabla_x(U_J + U_P)_{\text{rot}}. $$

Additionally, it follows from Eq. (B.3) that $\mathbf{x}_{\text{rot}}$ is

$$ x = r(\cos u \cos u_J + \sin u \cos i \sin u_J) $$

$$ y = r(-\cos u \sin u_J + \sin u \cos i \cos u_J) $$

$$ z = r(\sin u \sin i). $$

(B.4)

## Appendix B.2. The perturbing accelerations

In the rotating frame, the accelerations due to the third body and to the degree-2 gravity field are respectively

$$ \nabla_x U_J = \frac{n_J^2}{r^2} \begin{pmatrix} 2x \\ -y \\ -z \end{pmatrix} $$

$$ \nabla_x U_P = -\frac{\mu_E R_E^2 J_2}{r^2} \begin{pmatrix} \frac{x}{r} \left( \frac{5g(x, y, z)}{r^2} - 14 \right) \\ \frac{y}{r} \left( \frac{5g(x, y, z)}{r^2} + 4 \right) \\ \frac{z}{r} \left( \frac{5g(x, y, z)}{r^2} + 10 \right) \end{pmatrix}. \qquad (B.5) $$

In the local orbital frame of the probe around Europa, the cross-track components of these accelerations are computed using Eqs. (B.3) and (B.4) in Eq. (B.5), so that

$$ N_J = -N_{J,0} \sin i (3 \cos u \sin 2u_J + 3 \sin u \cos i (1 - \cos 2u_J)) $$

$$ N_P = -N_{P,0} \sin i (3 \cos u \sin 2u_J + \sin u \cos i (5 - 3 \cos 2u_J)), $$

(B.6)

with

$$ N_{J,0} = \frac{1}{2} n_J^2 r $$

$$ N_{P,0} = \frac{3}{5} \frac{\mu_E}{r^2} \frac{R_E^2}{r^2} J_2. $$

In the case of near-circular orbits, $r \approx a$ and $n^2 a^3 = \mu_E$, which gives

$$ N_{J,0} = \frac{1}{2} n_J^2 a $$

$$ N_{P,0} = \frac{3}{5} n^2 \frac{R_E^2}{a} J_2. $$

## Appendix B.3. Expression of the derivatives

The cross track component of all the perturbing forces considered here is obtained by summing $N_J$ and $N_P$ from Eq. (B.6), so that
rotation about the i-th coordinate axis by an angle θ.

Additionally, the transformation between the inertial frame and the local orbital frame can be written as

$$ \mathbf{x}_{\text{loc}} = \mathbf{R}_3(u) \cdot \mathbf{R}_1(i) \cdot \mathbf{R}_3(\Omega) \cdot \mathbf{x}_{\text{in}}, $$

which yields

$$ \mathbf{x}_{\text{loc}} = \mathbf{R}_{\text{tot}}(u, i, u_J) \cdot \mathbf{x}_{\text{rot}} \qquad (B.3) $$

with $\mathbf{R}_{\text{tot}}(u, i, u_J) = \mathbf{R}_3(u) \cdot \mathbf{R}_1(i) \cdot \mathbf{R}_3(-u_J)$ and $u_J = n_J t - \Omega$.

Once the gradients of the potentials are computed in the rotating frame, they need to be converted to the local orbital frame, so that

$$ \nabla_x(U_J + U_P)_{\text{loc}} = \mathbf{R}_{\text{tot}}(u, i, u_J) \cdot \nabla_x(U_J + U_P)_{\text{rot}}. $$

Additionally, it follows from Eq. (B.3) that $\mathbf{x}_{\text{rot}}$ is

$$ x = r(\cos u \cos u_J + \sin u \cos i \sin u_J) $$

$$ y = r(-\cos u \sin u_J + \sin u \cos i \cos u_J) $$

$$ z = r(\sin u \sin i). $$

(B.4)

## Appendix B.2. The perturbing accelerations

In the rotating frame, the accelerations due to the third body and to the degree-2 gravity field are respectively

$$ \nabla_x U_J = \frac{n_J^2}{r^2} \begin{pmatrix} 2x \\ -y \\ -z \end{pmatrix} $$

$$ \nabla_x U_P = -\frac{\mu_E R_E^2 J_2}{r^2} \begin{pmatrix} \frac{x}{r} \left( \frac{5g(x, y, z)}{r^2} - 14 \right) \\ \frac{y}{r} \left( \frac{5g(x, y, z)}{r^2} + 4 \right) \\ \frac{z}{r} \left( \frac{5g(x, y, z)}{r^2} + 10 \right) \end{pmatrix}. \qquad (B.5) $$

In the local orbital frame of the probe around Europa, the cross-track components of these accelerations are computed using Eqs. (B.3) and (B.4) in Eq. (B.5), so that

$$ N_J = -N_{J,0} \sin i (3 \cos u \sin 2u_J + 3 \sin u \cos i (1 - \cos 2u_J)) $$

$$ N_P = -N_{P,0} \sin i (3 \cos u \sin 2u_J + \sin u \cos i (5 - 3 \cos 2u_J)), $$

(B.6)

with

$$ N_{J,0} = \frac{1}{2} n_J^2 r $$

$$ N_{P,0} = \frac{3}{5} \frac{\mu_E}{r^2} \frac{R_E^2}{r^2} J_2. $$

In the case of near-circular orbits, $r \approx a$ and $n^2 a^3 = \mu_E$, which gives

$$ N_{J,0} = \frac{1}{2} n_J^2 a $$

$$ N_{P,0} = \frac{3}{5} n^2 \frac{R_E^2}{a} J_2. $$

## Appendix B.3. Expression of the derivatives

The cross track component of all the perturbing forces considered here is obtained by summing $N_J$ and $N_P$ from Eq. (B.6), so that
$$N_J + N_P = -\sin i (3(N_{J,0} + N_{P,0})\cos u \sin 2u_J + \sin u \cos i ((3N_{J,0} + 5N_{P,0}) - 3(N_{J,0} + N_{P,0})\cos 2u_J)),$$

and which we can finally substitute into Eq. (B.1). For near-polar orbits ($\cos i \approx 0$), the time-derivative of the inclination (B.1) is then

$$\frac{di}{dt} = -4K_i \cos^2 u \sin 2u_J, \qquad (B.7)$$

with

$$K_i = \frac{3}{4} \frac{r \sin i}{h} (N_{J,0} + N_{P,0}) = \frac{3}{4} \frac{\sin i}{h} \left( \frac{1}{2} n_J^2 a^2 + \frac{3}{5} n^2 R_E^2 J_2 \right).$$

Similarly, the time-derivative of the longitude of the ascending node can be written as

$$\frac{d\Omega}{dt} = -\frac{r}{h} (3(N_{J,0} + N_{P,0})\sin u \cos u \sin 2u_J + \sin^2 u \cos i ((3N_{J,0} + 5N_{P,0}) - 3(N_{J,0} + N_{P,0})\cos 2u_J))$$

The argument of latitude of the probe u and the argument of latitude of Jupiter around Europa uJ are respectively Tn-periodic and Dn-periodic. If we consider a m:R RGTO, we can use the periodicity condition Eq. (1), to isolate the secular term relevant for the long term evolution as

$$\begin{aligned} \left(\frac{d\Omega}{dt}\right)_{\text{secular}} &= \frac{1}{RT_n} \int_0^{RT_n} \frac{d\Omega}{dt} dt \\ &= -\frac{a}{2h} (3N_{J,0} + 5N_{P,0})\cos i \\ &= -\frac{3}{2n} \left(\frac{1}{2} n_J^2 + n^2 \left(\frac{R_E}{a}\right)^2 J_2\right) \cos i. \end{aligned} \qquad (B.8)$$

## References

Anderson, J.D., 1998. Europa's differentiated internal structure: inferences from four Galileo encounters. Science 281, 2019–2022. https://doi.org/10.1126/science.281.5385.2019.

Archinal, B.A., Acton, C.H., A'Hearn, M.F., Conrad, A., Consolmagno, G.J., Duxbury, T., Hestroffer, D., Hilton, J.L., Kirk, R.L., Klioner, S.A., McCarthy, D., Meech, K., Oberst, J., Ping, J., Seidelmann, P.K., Tholen, D.J., Thomas, P.C., Williams, I.P., 2018. Report of the IAU working group on cartographic coordinates and rotational elements: 2015. Celestial Mech. Dyn. Astron. 130. https://doi.org/10.1007/s10569-017-9805-5.

Arnold, D., Bertone, S., Jäggi, A., Beutler, G., Mervart, L., 2015. GRAIL gravity field determination using the Celestial Mechanics Approach. Icarus 261, 182–192. https://doi.org/10.1016/j.icarus.2015.08.015.

Asmar, S., Renzetti, N., 1993. The Deep Space Network as an Instrument for Radio Science Research (Technical Report).

Asmar, S.W., Armstrong, J.W., Iess, L., Tortora, P., 2005. Spacecraft Doppler tracking: noise budget and accuracy achievable in precision radio science observations. Radio Sci. 40. https://doi.org/10.1029/2004rs003101.

Bertone, S., Arnold, D., Girardin, V., Lasser, M., Meyer, U., Jäggi, A., 2021. Assessing reduced-dynamic parametrizations for GRAIL orbit determination and the recovery of independent lunar gravity field solutions. Earth Space Sci. 8. https://doi.org/10.1029/2020ea001454.

Beutler, G., 2005. Methods of Celestial Mechanics, vol. I. Springer Berlin, Heidelberg Physical, Mathematical, and Numerical Principles. https://doi.org/10.1007/b138225.

Beutler, G., Jäggi, A., Mervart, L., Meyer, U., 2010. The celestial mechanics approach: theoretical foundations. J. Geodes. 84, 605–624. https://doi.org/10.1007/s00190-010-0401-7.

Bignami, G., Cargill, P., Schutz, B., Turon, C., 2005. Cosmic Vision: Space Science for Europe 2015-2025. ESA Publ. Division.

Blanc, M., Prieto-Ballesteros, O., André, N., Gomez-Elvira, J., Jones, G., Sterken, V., Desprats, W., Gurvits, L.I., Khurana, K., Balmino, G., Blöcker, A., Broquet, R., Bunce, E., Cavel, C., Choblet, G., Colins, G., Coradini, M., Cooper, J., Dirkx, D., Fontaine, D., Garnier, P., Gaudin, D., Hartogh, P., Hussmann, H., Genova, A., Iess, L., Jäggi, A., Kempf, S., Krupp, N., Lara, L., Lasue, J., Lainey, V., Leblanc, F., Lebreton, J.P., Longobardo, A., Lorenz, R., Martins, P., Martins, Z., Marty, J.C., Masters, A., Mimoun, D., Palumba, E., Parro, V., Regnier, P., Saur, J., Schutte, A., Sittler, E.C., Spohn, T., Srama, R., Stephan, K., Szegő, K., Tosi, F., Vance, S., Wagner, R., Hoolst, T.V., Volwerk, M., Wahlund, J.E., Westall, F., Wurz, P., 2020. Joint Europa Mission (JEM): a multi-scale study of Europa to characterize its

habitability and search for extant life. Planet. Space Sci. 193, 104960. https://doi.org/10.1016/j.pss.2020.104960.

Bonanno, C., Milani, A., 2002. Symmetries and rank deficiency in the orbit determination around another planet. Celestial Mech. Dyn. Astron. 83, 17–33. https://doi.org/10.1023/a:1020110226569.

Böttcher, T., Huber, L., Corre, L.L., Leitner, J., McCarthy, D., Nilsson, R., Teixeira, C., Araujo, S.V., Wilson, R.C., Adjali, F., Altenburg, M., Briani, G., Buchas, P., Postollec, A.L., Meier, T., 2009. The HADES mission concept – astrobiological survey of Jupiter’s icy moon Europa. Int. J. Astrobiol. 8, 321–329. https://doi.org/10.1017/s147355040999022x.

Boutonnet, A., Varga, G., 2020. JUICE: Definition of the GCO with Cycles. ESA. Technical Report. ftp://spiftp.esac.esa.int/data/ANCDR/JUICE/man/JUICE_CReMA_Issue_4_Revision_2/Memo_20200429_JUICE_GCO_Cycles.pdf.

Cappuccio, P., Hickey, A., Durante, D., Benedetto, M.D., Iess, L., Marchi, F.D., Plainaki, C., Milillo, A., Mura, A., 2020. Ganymede’s gravity, tides and rotational state from JUICE’s 3GM experiment simulation. Planet. Space Sci. 187, 104902. https://doi.org/10.1016/j.pss.2020.104902.

Cinelli, M., Ciri, C., Ortore, E., 2015. Polynomial equations for science orbits around Europa. Celestial Mech. Dyn. Astron. 122, 199–212. https://doi.org/10.1007/s10569-015-9616-5.

Clark, K., Tan-Wang, G., Boldt, J., Greeley, R., Jun, I., Lock, R., Ludwinski, J., Pappalardo, R., Houten, T.V., Yan, T., 2009. Return to Europa: overview of the jupiter Europa orbiter mission. Adv. Space Res. 48, 629–650. https://doi.org/10.1109/aero.2009.4839315.

Dach, R., Lutz, S., Walser, P., Fridez, P., 2015. Bernese GNSS Software Version 5.2. University of Bern, Bern Open Publishing. https://doi.org/10.7892/boris.72297.

Folkner, W.M., Williams, J.G., Boggs, D.H., Park, R.S., Kuchynka, P., 2014. The planetary and lunar ephemerides DE430 and DE431. Interplan. Netw. Prog. Rep. 196, 42–196.

van Gelderen, M., Koop, R., 1997. The use of degree variances in satellite gradiometry. J. Geodes. 71, 337–343. https://doi.org/10.1007/s001900050101.

Grasset, O., Dougherty, M.K., Coustenis, A., Bunce, E.J., Erd, C., Titov, D., Blanc, M., Coates, A., Drossart, P., Fletcher, L.N., Hussmann, H., Jaumann, R., Krupp, N., Lebreton, J.P., Prieto-Ballesteros, O., Tortora, P., Tosi, F., Hoolst, T.V., 2013. JUpiter ICy moons Explorer (JUICE): an ESA mission to orbit Ganymede and to characterise the Jupiter system. Planet. Space Sci. 78, 1–21. https://doi.org/10.1016/j.pss.2012.12.002.

Jacobson, R.A., 2013. JUP310. ftp://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/satellites/a_old_versions/jup310.bsp.

Kaula, W.M., Street, R.E., 1967. Theory of satellite geodesy: applications of satellites to geodesy. Blaisdell 20, 101–101. https://doi.org/10.1063/1.3033941.
$$N_J + N_P = -\sin i (3(N_{J,0} + N_{P,0})\cos u \sin 2u_J + \sin u \cos i ((3N_{J,0} + 5N_{P,0}) - 3(N_{J,0} + N_{P,0})\cos 2u_J)),$$

and which we can finally substitute into Eq. (B.1). For near-polar orbits ($\cos i \approx 0$), the time-derivative of the inclination (B.1) is then

$$\frac{di}{dt} = -4K_i \cos^2 u \sin 2u_J, \qquad (B.7)$$

with

$$K_i = \frac{3}{4} \frac{r \sin i}{h} (N_{J,0} + N_{P,0}) = \frac{3}{4} \frac{\sin i}{h} \left( \frac{1}{2} n_J^2 a^2 + \frac{3}{5} n^2 R_E^2 J_2 \right).$$

Similarly, the time-derivative of the longitude of the ascending node can be written as

$$\frac{d\Omega}{dt} = -\frac{r}{h} (3(N_{J,0} + N_{P,0})\sin u \cos u \sin 2u_J + \sin^2 u \cos i ((3N_{J,0} + 5N_{P,0}) - 3(N_{J,0} + N_{P,0})\cos 2u_J))$$

The argument of latitude of the probe u and the argument of latitude of Jupiter around Europa uJ are respectively Tn-periodic and Dn-periodic. If we consider a m:R RGTO, we can use the periodicity condition Eq. (1), to isolate the secular term relevant for the long term evolution as

$$\begin{aligned} \left(\frac{d\Omega}{dt}\right)_{\text{secular}} &= \frac{1}{RT_n} \int_0^{RT_n} \frac{d\Omega}{dt} dt \\ &= -\frac{a}{2h} (3N_{J,0} + 5N_{P,0})\cos i \\ &= -\frac{3}{2n} \left(\frac{1}{2} n_J^2 + n^2 \left(\frac{R_E}{a}\right)^2 J_2\right) \cos i. \end{aligned} \qquad (B.8)$$

## References

Anderson, J.D., 1998. Europa's differentiated internal structure: inferences from four Galileo encounters. Science 281, 2019–2022. https://doi.org/10.1126/science.281.5385.2019.

Archinal, B.A., Acton, C.H., A'Hearn, M.F., Conrad, A., Consolmagno, G.J., Duxbury, T., Hestroffer, D., Hilton, J.L., Kirk, R.L., Klioner, S.A., McCarthy, D., Meech, K., Oberst, J., Ping, J., Seidelmann, P.K., Tholen, D.J., Thomas, P.C., Williams, I.P., 2018. Report of the IAU working group on cartographic coordinates and rotational elements: 2015. Celestial Mech. Dyn. Astron. 130. https://doi.org/10.1007/s10569-017-9805-5.

Arnold, D., Bertone, S., Jäggi, A., Beutler, G., Mervart, L., 2015. GRAIL gravity field determination using the Celestial Mechanics Approach. Icarus 261, 182–192. https://doi.org/10.1016/j.icarus.2015.08.015.

Asmar, S., Renzetti, N., 1993. The Deep Space Network as an Instrument for Radio Science Research (Technical Report).

Asmar, S.W., Armstrong, J.W., Iess, L., Tortora, P., 2005. Spacecraft Doppler tracking: noise budget and accuracy achievable in precision radio science observations. Radio Sci. 40. https://doi.org/10.1029/2004rs003101.

Bertone, S., Arnold, D., Girardin, V., Lasser, M., Meyer, U., Jäggi, A., 2021. Assessing reduced-dynamic parametrizations for GRAIL orbit determination and the recovery of independent lunar gravity field solutions. Earth Space Sci. 8. https://doi.org/10.1029/2020ea001454.

Beutler, G., 2005. Methods of Celestial Mechanics, vol. I. Springer Berlin, Heidelberg Physical, Mathematical, and Numerical Principles. https://doi.org/10.1007/b138225.

Beutler, G., Jäggi, A., Mervart, L., Meyer, U., 2010. The celestial mechanics approach: theoretical foundations. J. Geodes. 84, 605–624. https://doi.org/10.1007/s00190-010-0401-7.

Bignami, G., Cargill, P., Schutz, B., Turon, C., 2005. Cosmic Vision: Space Science for Europe 2015-2025. ESA Publ. Division.

Blanc, M., Prieto-Ballesteros, O., André, N., Gomez-Elvira, J., Jones, G., Sterken, V., Desprats, W., Gurvits, L.I., Khurana, K., Balmino, G., Blöcker, A., Broquet, R., Bunce, E., Cavel, C., Choblet, G., Colins, G., Coradini, M., Cooper, J., Dirkx, D., Fontaine, D., Garnier, P., Gaudin, D., Hartogh, P., Hussmann, H., Genova, A., Iess, L., Jäggi, A., Kempf, S., Krupp, N., Lara, L., Lasue, J., Lainey, V., Leblanc, F., Lebreton, J.P., Longobardo, A., Lorenz, R., Martins, P., Martins, Z., Marty, J.C., Masters, A., Mimoun, D., Palumba, E., Parro, V., Regnier, P., Saur, J., Schutte, A., Sittler, E.C., Spohn, T., Srama, R., Stephan, K., Szegő, K., Tosi, F., Vance, S., Wagner, R., Hoolst, T.V., Volwerk, M., Wahlund, J.E., Westall, F., Wurz, P., 2020. Joint Europa Mission (JEM): a multi-scale study of Europa to characterize its

habitability and search for extant life. Planet. Space Sci. 193, 104960. https://doi.org/10.1016/j.pss.2020.104960.

Bonanno, C., Milani, A., 2002. Symmetries and rank deficiency in the orbit determination around another planet. Celestial Mech. Dyn. Astron. 83, 17–33. https://doi.org/10.1023/a:1020110226569.

Böttcher, T., Huber, L., Corre, L.L., Leitner, J., McCarthy, D., Nilsson, R., Teixeira, C., Araujo, S.V., Wilson, R.C., Adjali, F., Altenburg, M., Briani, G., Buchas, P., Postollec, A.L., Meier, T., 2009. The HADES mission concept – astrobiological survey of Jupiter’s icy moon Europa. Int. J. Astrobiol. 8, 321–329. https://doi.org/10.1017/s147355040999022x.

Boutonnet, A., Varga, G., 2020. JUICE: Definition of the GCO with Cycles. ESA. Technical Report. ftp://spiftp.esac.esa.int/data/ANCDR/JUICE/man/JUICE_CReMA_Issue_4_Revision_2/Memo_20200429_JUICE_GCO_Cycles.pdf.

Cappuccio, P., Hickey, A., Durante, D., Benedetto, M.D., Iess, L., Marchi, F.D., Plainaki, C., Milillo, A., Mura, A., 2020. Ganymede’s gravity, tides and rotational state from JUICE’s 3GM experiment simulation. Planet. Space Sci. 187, 104902. https://doi.org/10.1016/j.pss.2020.104902.

Cinelli, M., Ciri, C., Ortore, E., 2015. Polynomial equations for science orbits around Europa. Celestial Mech. Dyn. Astron. 122, 199–212. https://doi.org/10.1007/s10569-015-9616-5.

Clark, K., Tan-Wang, G., Boldt, J., Greeley, R., Jun, I., Lock, R., Ludwinski, J., Pappalardo, R., Houten, T.V., Yan, T., 2009. Return to Europa: overview of the jupiter Europa orbiter mission. Adv. Space Res. 48, 629–650. https://doi.org/10.1109/aero.2009.4839315.

Dach, R., Lutz, S., Walser, P., Fridez, P., 2015. Bernese GNSS Software Version 5.2. University of Bern, Bern Open Publishing. https://doi.org/10.7892/boris.72297.

Folkner, W.M., Williams, J.G., Boggs, D.H., Park, R.S., Kuchynka, P., 2014. The planetary and lunar ephemerides DE430 and DE431. Interplan. Netw. Prog. Rep. 196, 42–196.

van Gelderen, M., Koop, R., 1997. The use of degree variances in satellite gradiometry. J. Geodes. 71, 337–343. https://doi.org/10.1007/s001900050101.

Grasset, O., Dougherty, M.K., Coustenis, A., Bunce, E.J., Erd, C., Titov, D., Blanc, M., Coates, A., Drossart, P., Fletcher, L.N., Hussmann, H., Jaumann, R., Krupp, N., Lebreton, J.P., Prieto-Ballesteros, O., Tortora, P., Tosi, F., Hoolst, T.V., 2013. JUpiter ICy moons Explorer (JUICE): an ESA mission to orbit Ganymede and to characterise the Jupiter system. Planet. Space Sci. 78, 1–21. https://doi.org/10.1016/j.pss.2012.12.002.

Jacobson, R.A., 2013. JUP310. ftp://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/satellites/a_old_versions/jup310.bsp.

Kaula, W.M., Street, R.E., 1967. Theory of satellite geodesy: applications of satellites to geodesy. Blaisdell 20, 101–101. https://doi.org/10.1063/1.3033941.
Kivelson, M.G., Khurana, K.K., Russell, C.T., Volwerk, M., Walker, R.J., Zimmer, C., 2000. Galileo magnetometer measurements: a stronger case for a subsurface ocean at Europa. Science 289, 1340–1343. https://doi.org/10.1126/science.289.5483.1340.

Lara, M., Russell, R., 2007. Computation of a science orbit about Europa. J. Guid. Control Dynam. 30, 259–263. https://doi.org/10.2514/1.22493.

Lemoine, F.G., Goossens, S., Sabaka, T.J., Nicholas, J.B., Mazarico, E., Rowlands, D.D., Loomis, B.D., Chinn, D.S., Neumann, G.A., Smith, D.E., Zuber, M.T., 2014. GRGM900c: a degree 900 lunar gravity model from GRAIL primary and extended mission data. Geophys. Res. Lett. 41, 3382–3389. https://doi.org/10.1002/2014GL060027.

Mazarico, E., Genova, A., Neumann, G.A., Smith, D.E., Zuber, M.T., 2015. Simulated recovery of Europa's global shape and tidal Love numbers from altimetry and radio tracking during a dedicated flyby tour. Geophys. Res. Lett. 42, 3166–3173. https://doi.org/10.1002/2015gl063224.

Moyer, T.D., 2003. Formulation for Observed and Computed Values of Deep Space Network Data Types for Navigation. https://doi.org/10.1002/0471728470.

National Research Council, 2011. Vision and Voyages for Planetary Science in the Decade 2013–2022. National Academies Press.

Ortore, E., Circi, C., Cinelli, M., 2015. Optimal orbits around Ganymede for the JUICE mission. Aero. Sci. Technol. 46, 282–286. https://doi.org/10.1016/j.ast.2015.07.021.

Park, R.S., Bills, B., Buffington, B.B., Folkner, W.M., Konopliv, A.S., Martin-Mur, T.J., Mastrodemos, N., McElrath, T.P., Riedel, J.E., Watkins, M.M., 2015. Improved detection of tides at Europa with radiometric and optical tracking during flybys. Planet. Space Sci. 112, 10–14. https://doi.org/10.1016/j.pss.2015.04.005.

Pavlak, T.A., 2013. Trajectory Design and Orbit Maintenance Strategies in Multi-Body Dynamical Regimes. Ph.D. thesis. Purdue University.

Petit, G., Luzum, B., 2010. IERS Conventions (2010). IERS Technical Note 36,1.

Phillips, C.B., Pappalardo, R.T., 2014. Europa clipper mission concept: exploring Jupiter's ocean moon. EOS Transac. Am. Geophy. Union 95, 165–167. https://doi.org/10.1002/2014eo200002.

Roth, L., Saur, J., Retherford, K.D., Strobel, D.F., Feldman, P.D., McGrath, M.A., Nimmo, F., 2014. Transient water vapor at Europa's south Pole. Science 343, 171–174. https://doi.org/10.1126/science.1247051.

Russell, R.P., 2006. Global search for planar and three-dimensional periodic orbits near Europa. J. Astronaut. Sci. 54, 199–226. https://doi.org/10.1007/bf03256483.

Sneeuw, N., 2000. A Semi-analytical Approach to Gravity Field Analysis from Satellite Observations. Technische Universität München. Dissertation.

Verma, A.K., Margot, J.L., 2018. Expected precision of Europa Clipper gravity measurements. Icarus 314, 35–49. https://doi.org/10.1016/j.icarus.2018.05.018.

Wahr, J.M., Zuber, M.T., Smith, D.E., Lunine, J.I., 2006. Tides on Europa, and the thickness of Europa's icy shell. J. Geophys. Res.: Planets 111. https://doi.org/10.1029/2006je002729.

Wu, X., Bar-Sever, Y.E., Folkner, W.M., Williams, J.G., Zumberge, J.F., 2001. Probing Europa's hidden ocean from tidal effects on orbital dynamics. Geophys. Res. Lett. 28, 2245–2248. https://doi.org/10.1029/2000gl012814.

Zannoni, M., Tortora, P., 2013. Numerical error in interplanetary orbit determination software. J. Guid. Control Dynam. 36, 1008–1018. https://doi.org/10.2514/1.59294.
Kivelson, M.G., Khurana, K.K., Russell, C.T., Volwerk, M., Walker, R.J., Zimmer, C., 2000. Galileo magnetometer measurements: a stronger case for a subsurface ocean at Europa. Science 289, 1340–1343. https://doi.org/10.1126/science.289.5483.1340.

Lara, M., Russell, R., 2007. Computation of a science orbit about Europa. J. Guid. Control Dynam. 30, 259–263. https://doi.org/10.2514/1.22493.

Lemoine, F.G., Goossens, S., Sabaka, T.J., Nicholas, J.B., Mazarico, E., Rowlands, D.D., Loomis, B.D., Chinn, D.S., Neumann, G.A., Smith, D.E., Zuber, M.T., 2014. GRGM900c: a degree 900 lunar gravity model from GRAIL primary and extended mission data. Geophys. Res. Lett. 41, 3382–3389. https://doi.org/10.1002/2014GL060027.

Mazarico, E., Genova, A., Neumann, G.A., Smith, D.E., Zuber, M.T., 2015. Simulated recovery of Europa's global shape and tidal Love numbers from altimetry and radio tracking during a dedicated flyby tour. Geophys. Res. Lett. 42, 3166–3173. https://doi.org/10.1002/2015gl063224.

Moyer, T.D., 2003. Formulation for Observed and Computed Values of Deep Space Network Data Types for Navigation. https://doi.org/10.1002/0471728470.

National Research Council, 2011. Vision and Voyages for Planetary Science in the Decade 2013–2022. National Academies Press.

Ortore, E., Circi, C., Cinelli, M., 2015. Optimal orbits around Ganymede for the JUICE mission. Aero. Sci. Technol. 46, 282–286. https://doi.org/10.1016/j.ast.2015.07.021.

Park, R.S., Bills, B., Buffington, B.B., Folkner, W.M., Konopliv, A.S., Martin-Mur, T.J., Mastrodemos, N., McElrath, T.P., Riedel, J.E., Watkins, M.M., 2015. Improved detection of tides at Europa with radiometric and optical tracking during flybys. Planet. Space Sci. 112, 10–14. https://doi.org/10.1016/j.pss.2015.04.005.

Pavlak, T.A., 2013. Trajectory Design and Orbit Maintenance Strategies in Multi-Body Dynamical Regimes. Ph.D. thesis. Purdue University.

Petit, G., Luzum, B., 2010. IERS Conventions (2010). IERS Technical Note 36,1.

Phillips, C.B., Pappalardo, R.T., 2014. Europa clipper mission concept: exploring Jupiter's ocean moon. EOS Transac. Am. Geophy. Union 95, 165–167. https://doi.org/10.1002/2014eo200002.

Roth, L., Saur, J., Retherford, K.D., Strobel, D.F., Feldman, P.D., McGrath, M.A., Nimmo, F., 2014. Transient water vapor at Europa's south Pole. Science 343, 171–174. https://doi.org/10.1126/science.1247051.

Russell, R.P., 2006. Global search for planar and three-dimensional periodic orbits near Europa. J. Astronaut. Sci. 54, 199–226. https://doi.org/10.1007/bf03256483.

Sneeuw, N., 2000. A Semi-analytical Approach to Gravity Field Analysis from Satellite Observations. Technische Universität München. Dissertation.

Verma, A.K., Margot, J.L., 2018. Expected precision of Europa Clipper gravity measurements. Icarus 314, 35–49. https://doi.org/10.1016/j.icarus.2018.05.018.

Wahr, J.M., Zuber, M.T., Smith, D.E., Lunine, J.I., 2006. Tides on Europa, and the thickness of Europa's icy shell. J. Geophys. Res.: Planets 111. https://doi.org/10.1029/2006je002729.

Wu, X., Bar-Sever, Y.E., Folkner, W.M., Williams, J.G., Zumberge, J.F., 2001. Probing Europa's hidden ocean from tidal effects on orbital dynamics. Geophys. Res. Lett. 28, 2245–2248. https://doi.org/10.1029/2000gl012814.

Zannoni, M., Tortora, P., 2013. Numerical error in interplanetary orbit determination software. J. Guid. Control Dynam. 36, 1008–1018. https://doi.org/10.2514/1.59294.