Decoupled and coupled moons' ephemerides estimation strategies -- Application to the JUICE mission

When reconstructing natural satellites' ephemerides from space missions' tracking data, the dynamics of the spacecraft and natural bodies are often solved for separately, in a decoupled manner. Alternatively, the ephemeris generation and spacecraft orbit determination can be performed concurrently. This method directly maps the available data set to the estimated parameters' covariances while fully accounting for all dynamical couplings. It thus provides a statistically consistent solution to the estimation problem, whereas this is not directly ensured with the decoupled strategy. For the Galilean moons in particular, the JUICE mission provides a unique opportunity for ephemerides improvement. For such a dynamically coupled problem, choosing between the two strategies will be influential. This paper provides a detailed, explicit formulation for the coupled approach, before comparing the performances of the two state estimation methods for the JUICE test case. To this end, we used both decoupled and coupled models on simulated JUICE radiometric data. We compared the resulting covariances for the Galilean moons' states, and showed that the decoupled approach yields slightly lower formal errors for the moons' tangential positions. On the other hand, the coupled model can reduce the state uncertainties by more than one order of magnitude in the radial direction. It also proved more sensitive to the dynamical coupling between Io, Europa and Ganymede, allowing the solutions for the first two moons to fully benefit from JUICE orbital phase around Ganymede. However, many issues remain to be solved before a concurrent estimation strategy can be successfully applied to reconstruct the moons' dynamics over long timescales. Nonetheless, our analysis highlights promising ephemerides improvements and thus motivates future efforts to reach a coupled state solution for the Galilean moons.


Introduction
The upcoming JUICE mission 1 (JUpiter ICy moons Explorer) will focus on the three Galilean moons Europa, Ganymede and Callisto. The JUICE spacecraft is expected to arrive in the Jovian system in 2031, with a launch planned in 2023. It will first execute a series of flybys (2, 7 and 21 flybys at Europa, Ganymede and Callisto, respectively), from 2032 to 2034. JUICE will then initiate its orbital phase around Ganymede, with an eccentricity ranging from 0.6 to 0 (GEO and GCO500 phases: Ganymede Elliptic Orbit and Ganymede Circular Orbit with an altitude of 5000 km, respectively). In May 2035, after a second elliptical phase, the spacecraft will eventually enter its final circular orbit at 500 km altitude (denoted GCO500), for a nominal period of 4 months. This mission profile, displayed in Figure 1 and adopted in the rest of this paper, was obtained from the version 5.0 of the CReMA (Consolidated Report on the Mission Analysis) 2 .
The JUICE spacecraft will carry one dedicated radio science instrument (3GM: Gravity and Geophysics of Jupiter and the Galilean Moons, e.g. Di , which will provide highly accurate range and Doppler measurements (see Section 3.3). These 3GM observations will be complemented by the PRIDE experiment (Planetary Radio Interferometry and Doppler Experiment, e.g. Gurvits et al., 2013). The latter does not require any additional onboard instrument and uses tracking or 3GM radiometric signals to derive angular position measurements of the spacecraft with respect to the ICRF (International Celestial Reference Frame), as well as supplementary Doppler observables (e.g. Duev et al., 2012;Bocanegra-Bahamón et al., 2018;Molera Calvés et al., 2021). The radiometric data to be acquired by both 3GM and PRIDE are expected to contribute to a more accurate determination of the Galilean moons' states (Dirkx et al., 2016(Dirkx et al., , 2017Lari and Milani, 2019;Cappuccio et al., 2020). Improved ephemerides are crucial to better understand the long-term thermal-orbital evolution of these moons, which is strongly driven by tidal dissipation in both Jupiter and the satellites themselves (Peale, 1999;Hussmann and Spohn, 2004;Greenberg, 2010;Hay et al., 2020). The moons' ephemerides provide a natural way to extract the current rates of tidal dissipation, through the observed migration rates of the satellites (e.g. Lainey and Tobie, 2005;Lainey et al., 2009). Furthermore, a better characterisation of tidal dissipation mechanisms can provide tighter constraints on the moons' interiors, which is critical to investigate sub-surface ocean's properties (or confirm the existence of a putative ocean for Callisto, e.g. Lunine, 2017).
For natural satellites' ephemerides, the estimations of the spacecraft's and moons' dynamics are typically not performed in a coupled manner (e.g. Rosenblatt et al., 2008;Durante et al., 2019). Instead, when ephemerides are to be determined from flybys, the spacecraft trajectory with respect to the central body (i.e. body at which the flyby is performed) is determined from the available tracking data, along with the central body's state at the flyby epoch. The per-flyby state solutions for the natural body define the so-called normal points, which are then used in a second global estimation to reconstruct the longterm dynamics of this body (e.g. Durante et al., 2019). If needed, a unified model may also be used, in which the spacecraft dynamics are determined in a multi-arc fashion, and the natural bodies' dynamics in a single-arc fashion, during a single inversion. Such an approach, used for instance by Dirkx et al. (2018); Lari and Milani (2019), has the advantage of automatically incorporating all dynamical couplings, as well as the full sensitivity to physical parameters of both types of dynamics. The decoupled strategy, on the other hand, reconstructs the moons' dynamics from the normal points, which only capture the moons' kinematics (and not dynamics) at each flyby. However, while desirable, a coupled solution for the spacecraft's and moons' states is not always achievable in practice (e.g. Durante et al., 2019). It indeed requires the moons' dynamical models (with respect to the planet) to be consistent over both short and long timescales, to the accuracy level of the spacecraft's dynamics (with respect to the moon). Here, short and long timescales respectively refer to typical flyby duration (i.e. a few hours) and entire mission timeline (i.e. several years, so still short with respect to system evolution's timescale).
For the JUICE mission, we will be presented with a unique situation: the mission profile indeed calls for a combination of flybys around multiple satellites and an extended orbit phase around Ganymede, which was never before performed in a planetary mission. Additionally, three of the four Galilean moons are in resonance, making the estimation of the different moons' dynamics strongly coupled, with the added complication that JUICE will provide a strong imbalance in data for these three moons. As a result, the estimation of ephemerides from JUICE-only data is (close to be) an ill-posed mathematical problem (Dirkx et al., 2017). Due to the complexity and novelty of the mission profile, and to the strong dynamical couplings that are involved, the concurrent single-and multi-arc estimation strategy appears particularly well-suited for the JUICE test case. In this paper, we compare the simulated state estimation solutions obtained with both the decoupled and coupled approaches, to quantify the impact of the adopted estimation strategy.
We limit ourselves to a covariance analysis, complemented by a deterministic simulation performed as verification. As already highlighted, the practical applicability of the coupled method to the JUICE mission is however not guaranteed, as bringing the dynamical models fidelity down to the required accuracy level will be really challenging. By definition, these issues cannot be addressed by a covariance analysis, as the resulting formal uncertainties do not account for inaccuracies in the dynamical models used for the moons and the JUICE spacecraft, or in the models representing the observations' errors. Our paper thus assesses which uncertainty levels could be obtained with a coupled estimation, provided that our dynamical models are accurate enough for a viable solution to be achieved. The formal uncertainties we obtain therefore quantify the coupled strategy's requirements in terms of dynamical modelling accuracy. Besides comparison pur-poses, precisely evaluating the performance of the decoupled method is thus also crucial in case the obtention of a global coupled solution for the Galilean moons remains beyond (current) modelling capabilities. It must be noted that modelling issues, while not directly addressed by our covariance analysis, remain nonetheless deeply relevant for this study and will therefore be extensively discussed (see Sections 2.4 and 5.2). More generally, the limitations and scope of our formal analysis will be further detailed in Section 2.4.
The details of the coupled model, as well as the issues associated with its implementation and application, are not found in the literature, despite its application in a number of past studies (e.g. Dirkx et al., 2018;Lari and Milani, 2019;Magnanini, 2021). Therefore, we choose to provide a detailed exposition of our coupled method in Section 2, completed by a shorter description of the decoupled approach. The models used to either propagate the moons' and spacecraft's dynamics or simulate the JUICE radiometric observations are then described in Section 3. Section 4 presents the results of our comparative analysis of the coupled and decoupled estimation strategies, first for the flyby phase only, and then for the entire JUICE mission. Finally, Section 5 discusses in more detail the strengths and challenges of both estimation methods, before our conclusions are summarised in Section 6.

Estimation framework
This section describes the whole estimation process, for both the coupled and decoupled approaches introduced in Section 1. The complete formulation for the coupled single-and multi-arc state estimation model, still missing in the literature, is provided in Section 2.2. For the sake of completeness, our implementation of the decoupled strategy for the JUICE case is given in Section 2.3. This section thus directly highlights the main differences between the two estimation strategies.
All methods described in the following were implemented in our TU Delft Astrodynamics Toolbox (Tudat) software 3 , and are therefore freely available.

Covariance analysis
We first briefly review the propagation of the variational equations and describe how covariance matrices are generated and propagated in our simulations, as typically implemented in any estimation process, either singleand/or multi-arc.

Variational equations formulation
The variational equations describe how the dynamics of the system are influenced by the parameters to be estimated. In the following, we adopt the nomenclature of 3 Documentation: https://tudat-space.readthedocs.io Full source code: https://github.com/tudat-team/tudat-bundle Montenbruck and Gill (2000). The state vector is denoted as y and is propagated numerically from the initial time t 0 usingẏ where p is a vector of parameters influencing the system's dynamics or the observations, and f represents the dynamical model (described in Section 3). Unless otherwise indicated, all states are expressed in a reference frame with inertial orientation (e.g. J2000). We stress that, in a general formulation, the states y need not be translational states, but may be any type of dynamics, of any number of bodies (see Mazarico et al. (2017); Dirkx et al. (2019) for an example of coupled translational-rotational dynamics estimation of multiple bodies).
The state transition matrix Φ(t, t 0 ) and sensitivity matrix S(t) are defined as The differential equations used to solve for Φ and S are termed the variational equations, and are given by with the following initial conditions: where n and n p represent the sizes of the state vector y and parameter vector p, respectively. The single-arc and multi-arc formulations are essentially identical, with the sole difference that the multi-arc solution is obtained by subsequent, independent, integrations of Eqs. (1), (4) and (5). A variant of the multi-arc method, referred to as the constrained multi-arc approach, uses the fact that the arcwise state estimates obtained for a given body should be consistent to further constrain the estimation solution (Alessi et al., 2012;Serra et al., 2018). In our analysis, we however chose to limit ourselves to the unconstrained multi-arc estimation. During the first part of the mission, the flybys are indeed temporarily distant, such that propagating information from previous arcs would not efficiently constrain the JUICE spacecraft's state. When arcs are contiguous (i.e. orbital phase around Ganymede), the high quality of the estimation solution anyway undermines the use of multi-arc constraints.

Propagated covariance
Let h(T, q) denote the set of all modelled observations generated up to a time T . The design matrix H(T, q) associated with these observations is then formed by computing with q a vector containing the estimated parameters (e.g. Montenbruck and Gill, 2000;Milani and Gronchi, 2010). It usually includes initial states parameters, represented by the vector y 0 , and a subset of the parameters p influencing the dynamical or observational models. To simplify the notations, the vector of estimated parameters q will be divided as q = [y 0 ; p] in the following. It should however be noted that the exact definition of q depends on the estimation model, and that y 0 and p might not directly incorporate all initial states, dynamical and observational models. More details on how y 0 and p are precisely defined for both the coupled and decoupled models will be provided in Sections 2.2, 2.3 and 3.3. The covariance matrix of q obtained using data up to time T is denoted P qq (T ) and is given by where P qq,0 is the a priori covariance matrix of the parameters q (see Section 2.3.3 for a priori knowledge in our JUICE test case). The matrix W(T ) contains the weights associated with all observations up to time T . In most cases, it is set as a diagonal matrix with W ii = σ −2 h,i , implicitly assuming the measurement uncertainties to be uncorrelated. This is however not the case in every estimation step of the decoupled model, as will be discussed in Sections 2.3.1 and 2.3.3. σ h,i denotes the uncertainty of observation i. The covariance P qq (T ) can be used to compute the covariance of the state y at any later time t. We refer to this propagated covariance as P yy (t, T ) and define it as where Φ and S are the state transition and sensitivity matrices obtained through Eqs. 4, 5, 6 and 7. For the covariances in Eqs. (9) and (10), the formal errors are obtained from the square root of the diagonal elements of P qq and P yy , respectively.

Coupled single-and multi-arc estimation
This section describes the extension of the variational equations introduced in Section 2.1.1 to the concurrent estimation of single-and multi-arc states. The formulation specifics are detailed in Section 2.2.3 for the JUICE case.

General principle
The coupled strategy relies on the concurrent estimation of the spacecraft orbit and natural bodies' ephemerides, as well as of all parameters influencing the dynamics and/or observations (vector q in Section 2.1). The natural bodies' dynamics, on the other hand, are reconstructed over a single arc. More precisely, the spacecraft orbit is solved for in an arc-wise manner (along with any observations-or spacecraft-related parameters, e.g. biases, accelerometer calibrations factors, etc.). Such a coupled model allows us to directly and robustly link observation strategies, data quality, mission profile, etc. to the final formal uncertainties in natural bodies' ephemerides and dynamical parameters (tidal dissipation, gravity field coefficients). Figure 2a provides a schematic visualisation of this coupled approach, taking the JUICE flyby phase as an example.
A preliminary framework for a concurrent single-and multi-arc estimation was briefly described by Dirkx et al. (2018). A similar method was also used by Lari and Milani (2019) to study the onset of chaos in the dynamics of the JUICE spacecraft. The following sections provide the detailed and complete formulation for such a coupled estimation procedure for single-and multi-arc dynamics.

Coupled variational equations
We denote the single-arc state vector y S (t), of size n s and associated with the initial time t 0 . The multi-arc state, on the other hand, is designated by y M (t) and the arc-wise initial time is noted t i for arc i, with i = 1...N for N arcs. The size of each arc-wise state vector is n m .
The multi-arc state function y M (t) is defined as where y M,i (t) refers to the state at time t during arc i, and t i denotes the end time of arc i. It should be noted that the multi-arcs need not be contiguous, and gaps may exist in the arc-wise solutions to Eq. (1). Eq. (11) is therefore only defined if an arc i exists that satisfies Eq. (12) at the time t. The full state function is given as a combination of the single-and multi-arc states at time t, as follows: For our estimation, we require a linearised model for the change in y(t) induced by a variation in the parameters p and in the full vector of initial states y 0 , to compute S and Φ, respectively: Looking at Eqs. (13) and (14), the size of the full initial state vector to be estimated (y 0 ) is different from that of Flyby + 2 Flyby + 1

Flyby Moon
Step 1: JUICE's and moon's initial  the state function y(t), as the former combines all singlearc and multi-arc initial states, while y(t) only includes the single-arc states and the multi-arc state of the current arc. We note that p may affect the single-or multi-arc dynamics solutions, or both. The only limitation imposed on the dynamics (Eq. 1) is that the differential equation for y S (t) must be independent of y M (t). The opposite is not true, y M (t) being allowed to (and in our case does) depend on y S (t). These assumptions only hold if the masses of the multi-arc bodies are negligible with respect to the single-arc bodies'. This is typically the case as the spacecraft's dynamics are generally propagated in a multi-arc manner, while the bodies included in the single-arc solution are often natural bodies (see Section 2.2.3). We use Φ SS (t, t 0 ), Φ M M,i (t, t i ) to refer to the singleand multi-arc state transition matrices, respectively. Similarly, the single-and multi-arc sensitivity matrices are denoted as S S (t) and S M,i (t). We note that, for the multi-arc case, the parameter vector p can include local parameters that influence the dynamics of a single arc i only, as well as global parameters that affect all arcs.
The state transition matrix is defined as the derivative of the current state y(t) (Eq. 13) with respect to both the single-arc initial state y S (t 0 ) at time t 0 , and the arc-wise initial states y M,i , at the beginning t i of each arc. The full state transition matrix, noted Φ(t; t 0 , t i ), and sensitivity matrix S(t) can thus be written as , where we introduced the coupling term into the state transition matrix. The zero entries in the first row of Eq. (15) directly result from the dynamics of y S (t) being independent of y M (t). For clarification purposes, the dimensions of these zero blocks are specified as subscripts.
To obtain the numerical solution to the coupled variational equations, we first propagate the single-arc dynamics and variational equations, to obtain y S , Φ SS and S S . For the multi-arc formulation, the differential equations for Φ M M are unchanged compared to the classical (decoupled) approach. However, to compute the full state transition and sensitivity matrices, we need a formulation for the coupling term Φ M S,i which incorporates the influence of single-arc dynamics on the multi-arc dynamics, as follows: 5 Similarly, we require a formulation for S M S,i , given by This completes the required formulation for the differential equations governing the evolution of Eqs. (15) and (16). From the single-arc propagation, we retrieve y S , Φ SS and S S , appearing in Eqs (18) and (19), and use as a given when solving the multi-arc dynamics and variational equations.
The advantage of this approach, with two separate integrations to fully populate the coupled Φ and S matrices, is that different numerical settings may be used for the single-and multi-arc segment. In particular, for the case of coupled natural body and spacecraft dynamics estimation, one will typically require a much smaller time-step for propagating the spacecraft than for the natural bodies (as well as possibly a different integrator).

Formulation for the JUICE mission
In Section 2.2.2 we presented our general framework for propagating coupled single-and multi-arc variational equations. We now discuss specific details of the formulation for the JUICE mission. Our single-arc state vector is defined in a planetocentric reference frame as where the index 0 refers to Jupiter, and indices 1,2,3,4 correspond to Io, Europa, Ganymede and Callisto, respectively, following Dirkx et al. (2016).
Only the spacecraft's dynamics are solved for in arcwise manner, such that the multi-arc state vector for arc i can simply be written as We use 'sc' to denote properties relating to the JUICE spacecraft, while j i designates the index j of the central body during arc i. The reference frame origin is selected as the moon where a flyby is performed during the flyby phase (Europa, Ganymede or Callisto), and as Ganymede during the orbital phase. Solving the coupled variational equations provides solutions for the derivatives ∂x (j i ) sc,i ∂ * , which describe changes in the moon-centered state of the JUICE spacecraft. However, to evaluate our design matrix H (see Eq. 8), we need to account for the variations in the observed position of the spacecraft, often expressed in an inertial frame (e.g. Solar System Barycentre). As a result, the dynamics of the moons influence the observed position of the spacecraft in two distinct manners: • the dynamical contribution, through the bottom-left block of Eq. (15), • the kinematic or indirect contribution, through the variations in the moons' states with respect to the reference frame used for the observed spacecraft's motion.
This methodology automatically allows the incorporation of parameters that directly influence both the spacecraft's and moon's dynamics. Principally, this concerns the moons' spherical harmonic coefficients. Consistently propagating S(t) for the full system ensures that the covariance of the moons' initial states is robustly propagated to later epochs (see Section 2.1.2).

Decoupled single-and multi-arc estimation
To complement the description of the coupled estimation method in Section 2.2, the decoupled strategy is now discussed. As this approach does not differ from textbook formulations (e.g. Montenbruck and Gill, 2000;Milani and Gronchi, 2010), less details are provided and we directly address the JUICE case specifically. For our comparative analysis, it is however crucial to make both the decoupled and coupled formulations explicit, to highlight their main differences.

General principle
The decoupled estimation is performed in two separate steps, as shown in Figure 2b. The spacecraft's and natural bodies' dynamics are first solved for concurrently, as in the coupled case, but in a multi-arc manner. Only the dynamical coupling between the spacecraft and the central body is thus accounted for in this estimation step (while all dynamical couplings are included to propagate the moons' states, see Section 3.2). Since the natural bodies' states are independently estimated for each arc, the adopted dynamical model need not be consistent over long timescales.
This first estimation step therefore provides arc-wise estimated states for the central bodies. These so-called normal points are then used as observables in a second step, which aims at reconstructing the natural bodies' dynamics on a more global scale. More precisely, a normal point is defined as the central moon's cartesian state components (vector of size 6) and associated covariances (6by-6 matrices), determined with respect to Jupiter at the time of closest approach. The covariances P qq for the arc-wise initial states, resulting from the first estimation step, determine the weights W (see Eq. 9) assigned to each normal point in the second step. The matrix W is thus not exactly diagonal in this particular case (see Section 2.1). It instead shows non-zeros, diagonally-centered, 6-by-6 blocks containing the 6-by-6 normal points' covariances. The entire two-step decoupled estimation process is depicted in Figure 2b, using JUICE flybys as an example. While the coupled model estimates all parameters concurrently, different sets of estimated parameters are defined for the two steps of the decoupled method. An obvious example are the spacecraft's states, which are determined in a multi-arc manner in the first step but are absent from the second step, when reconstructing the global solution for the moons (see Figure 2b). As previously mentioned, it must be noted that the first step of the decoupled model only estimates the state of the central moon j when determining the normal point for a flyby around that moon. The state uncertainties of the other moons are not accounted for, and our decoupled estimation strategy might thus yield slightly too optimistic formal errors for the arc-wise state of moon j. However, as a verification, we ran an additional analysis for the JUICE test case including the other moons' states as consider parameters in the normal points determination process. As an indication, we provide the results obtained for two of the JUICE flybys in Appendix A. This verification showed that these uncertainties have a negligible impact on the normal points solution when using tracking arcs of 8 hours only around each flyby (see Section 3.3). This assumption should however be revisited if longer tracking arcs were to be considered for the JUICE flybys, as the influence of the state uncertainties for the non-central moons is then expected to increase. A more complete discussion on the estimated parameters for our JUICE analysis is provided in Section 3.3 (see Table 3).
2.3.2. Decoupled variational equations for the JUICE mission As described in Section 2.3.1, the spacecraft's and moon's dynamics are first reconstructed in a multi-arc fashion. For each arc i, the initial translational state to be estimated is thus defined as where j i again refers to the index of the central moon for arc i.
In practice, all arcs sharing the same central moon are combined, to allow some dynamical parameters to be estimated globally alongside the arc-wise states (e.g. gravity field coefficients of the central moon). For each moon j, the full initial state is thus built by concatenating the corresponding multi-arc states, as follows: ...
with N j the number of arcs with moon j as central body.
The arc-wise state transition matrix Φ i (t, t i ) can be derived from Eq. (22) as Eq. (25) shows some similarities with Eq. (15), but also clearly highlights major differences between the coupled and decoupled formulations. In particular, Φ sc,ji (t, t i ) also represents a coupling term, but expressed in a multi-arc fashion and with respect to the central moon j i only: as opposed to Eq. (17). The variational equations provided above apply to the first, arc-wise estimation step of the decoupled strategy (see Section 2.3.1). The second phase, in which a singlearc estimation is performed to reconstruct the long-term dynamics of the Galilean moons, follows the regular singlearc approach. The associated variational equations are therefore not detailed in this paper.

A priori knowledge strategy
As shown in Equation 9, prior knowledge is accounted for in the estimation by means of the a priori covariance matrix P qq,0 . Appropriate a priori values for all estimated parameters, referred to as default a priori covariances, are further discussed in Section 3.3 and are combined in a diagonal matrix P 0 (shortened notation for P qq,0 ).
For the decoupled case in particular, the moons' arcwise state solutions first determined at the beginning of each flyby strongly depend on these a priori constraints (see results in Section 4.1.3). Using the same default a priori values for all arc-wise moon states, derived from the existing ephemerides solutions, would be a rather conservative approach. It indeed neglects the iterative improvement achievable by progressively including more flybys in the estimation. Even if the observations processed by the estimation remain the same, some additional information is incorporated in the multi-arc model to improve the solution, namely that the arc-wise state solutions for a given moon j are not completely independent from one another. They indeed belong to a single body's trajectory and should thus be dynamically consistent. Such an update strategy for the a priori contraints on the moons' states can be compared to the multi-arc constrained approach for the spacecraft's orbit determination (e.g. Alessi et al., 2012), but applied to the moons' arc-wise states instead of the spacecraft's. It must be noted that using this a priori update strategy introduces some correlations between the arc-wise state components of moon j (i.e. between the different normal points determined for this moon). The off-diagonal blocks of the weight matrix W are therefore not filled with zeros anymore.
Focusing on the N j arcs with moon j as central body, more realistic a priori covariances can be derived for arc k by propagating the covariance obtained for arc k − 1 up to the beginning of arc k. This propagated covariance is denoted as P k→k+1 0 in the following. Some state components may nonetheless be poorly constrained by the previous arc's estimation, thus yielding unrealistically large a priori errors in certain directions. The a priori matrix P k 0 for arc k is thus built as a combination of the default and propagated a priori covariance matrices, as follows: is the state transition matrix for moon j, computed from the start of arc k − 1 to the beginning of the current arc k. This propagation scheme is initialised with the default a priori matrix (so P 0→1 0 = 0, i.e. matrix filled with zeros).
Iterating on the a priori knowledge for the moons' arcwise states requires to run the first step of the decoupled estimation multiple times, gradually increasing the number of arcs being processed. The final outcome of the multi-arc estimation (i.e. normal points and global parameters' estimates, see Section 2.3) is reached when all N j arcs associated with moon j are included. This process is schematically summarised in Figure 3.
It must be stressed that, in the strategy described above, we only propagate the covariances between moon j's own state components from arc k − 1 to the beginning of arc k. We thus neglect the influence that uncertainties in other moons' states could have on the propagated a priori covariance for moon j. Discarding the contribution of the other moons is consistent with the philosophy of the decoupled estimation, in which only the central moon's state is determined for each arc. It should however be noted that the a priori used for the normal points might therefore be slightly too optimistic.
The impact of the a priori information on the parameters solution, and especially the effect of the above updating strategy for the arc-wise states, will be further investigated and discussed in Section 4.1.3. For each estimated parameter, the contribution of the a priori information to the solution c q can be evaluated as follows (e.g. Floberghagen, 2001): where I is to the identify matrix, while P and P −1 0 refer to the final and a priori covariance matrices, respectively (P qq and P −1 qq,0 in Eq. 9). A c q equal to 1 indicates that the parameter's estimation relies entirely on the observations, while a value of 0 means that it is based on a priori information.

Flyby 3
Flyby 2 Flyby 1 Moon → → 2D projection of the a priori covariance for each arc 2D projection of the estimated covariance i for arc (colour-coded based on which a priori information is used) Covariance propagation Default a priori constraints (see Section 3.3) Updated a priori constraints ( i propagated to beginning of arc + 1) → + Figure 3: Schematic representation of the iterative strategy for the a priori covariances, applied to the first step of the decoupled method (i.e. normal points determination). The dashed ellipses represent the two-dimensional projection of the a priori covariance for the arc-wise position of the moon. The solid ellipses display the two-dimensional projection of the estimated covariance after inversion, for each arc. These solid covariance ellipses are colour-coded to indicate which a priori values are used. The default a priori information P 0 is represented in green. A dashed black arrow indicates a covariance propagation (i.e. update of the a priori covariance in this case). This representation only illustrates the update principle and is not to scale.

Scope of the comparative analysis
As mentioned in Section 1, we limit ourselves to a covariance analysis in this study. We compare the performances of the decoupled and coupled models by analysing the formal errors and correlations obtained in both cases. It is important to stress that the coupled model, by concurrently accounting for all dynamical couplings and sensitivities, directly maps the simulated observations set to the estimated parameters' covariances. Assuming that the fidelity of both our dynamical and measurement error models is sufficient, the resulting formal errors and correlations are therefore considered to provide a good statistical representation of the estimation solution, while this is not directly true for the decoupled method. Our study characterises how much the solution obtained by decoupling the spacecraft's and moons' state estimation departs from the covariances given by the coupled approach, regarded as statistically consistent.
For any estimation, true errors obtained from real data after completing the iterative least-squares estimation are however larger than the formal errors provided by a covariance analysis. Differences between true and formal errors originate from non-white measurement noise, as well as inaccuracies in the models used for the spacecraft's and planetary system's dynamics. For the JUICE mission data analysis, this situation may even be more severe than for any previous natural satellite's ephemeris determination, due to the much better data quality and subsequent higher requirements on dynamical modelling.
In practice, both the decoupled and coupled methods are limited by the dynamical model fidelity, so that the true errors would be larger than formal ones in the two cases. The goal of this study is to determine at which point and to what extent the coupled estimation would be beneficial for ephemerides determination, as well as to quantify the dynamical model requirements to achieve this (see Section 1). Dynamical mismodelling would nonetheless influence the decoupled and coupled solutions differently, and represent a major challenge for the applicability of the coupled model in particular. These modelling issues will therefore be discussed in more detail in Section 5.

Dynamical and Observation Models
In the following section, we discuss the settings and models used for our simulated covariance analysis for Galilean satellites' ephemerides from JUICE tracking data. The spacecraft's and moons' dynamical models are summarised in Sections 3.1 and 3.2, respectively, while the estimation and observations settings are described in Section 3.3.

Spacecraft Dynamics
When propagating the dynamics of the spacecraft during arc i, the following accelerations were taken into account: • spherical harmonic acceleration of the central moon j i , expanded up to degree l j and order m j (with l 1 /m 1 = 2/2, l 2 /m 2 = 4/4, l 3 /m 3 = 12/12, l 4 /m 4 = 6/6). Higher degrees might be accessible from JUICE data, especially for Ganymede, but were purposely not included in our analysis, which primarily focuses on ephemerides determination.
• point-mass acceleration of moon k, for each k ∈ {1, 2, 3, 4}, with k ̸ = j i • spherical harmonic acceleration of Jupiter, expanded up to degree l 0 = 8 and order 0 (zonal terms only), • point-mass accelerations exerted by Saturn and the Sun, • cannonball radiation pressure acceleration due to the Sun's radiation, • arc-wise constant (in RTN frame) empirical acceleration, representing errors in the accelerometer calibration.
We adopted the same models for the environment (gravity fields, ephemerides, etc.) as Dirkx et al. (2016), with a number of exceptions: we used the Jupiter gravity field from Iess et al. (2018), and CReMA 5.0 for the JUICE orbit 4 , released by ESA in the form of Spice kernels (Acton Jr, 1996).

Moon Dynamics
When propagating the dynamics of the Galilean moons, we used similar models as Lainey et al. (2004); Dirkx et al. (2016), taking into account • the mutual spherical harmonic acceleration between Jupiter and each moon j, with the gravity field of Jupiter expanded up to degree 8 and order 0, and that of the moons to degree and order 2, • the mutual spherical harmonic accelerations between all moons j and k, with the fields of the bodies expanded up to degree and order 2, • the point-mass accelerations due to Saturn and the Sun, • the acceleration exerted on each moon j due to tidal dissipation in Jupiter forced by moon j, • the acceleration on each moon j due to tidal dissipation in moon j forced by Jupiter. The influence of the tides raised by body k on body j on the system's dynamics was modelled by time-varying corrections applied to the spherical harmonics expansion of the body j' gravity field, as done in (e.g. Dirkx et al., 2016Dirkx et al., , 2017.

Estimation Settings
As the tracking configurations differ between the flyby and orbital phases of the JUICE mission, different estimation settings were used. An 8-hour tracking arc was defined for each flyby, centered at the time of closest approach. For the orbital phase, we simulated 8 hours of tracking per day. In practice, the JUICE spacecraft will be tracked from three stations of the European Space Tracking network (ESTRACK), the main one being Malargue, which is as of yet the only one enabling both X-and Ka-band tracking. However, we assumed that the other two will also be able to handle Ka-band tracking by the time the JUICE spacecraft arrives in the Jovian system. We thus considered 8 hours per day of almost continuous tracking, except during occultations or for elevations lower than 15 deg (as in e.g. Di Magnanini, 2021).
Tracking arcs of two days, separated by three days without tracking, were used as the nominal tracking configuration for the orbital phase. Nonetheless, the sensitivity of the estimation solution to these tracking settings was investigated by considering one day-and one weekarcs (results are presented in Section 4.2.3). The three days interval between two tracking arcs was merely used to reduce the computational load of our simulations. We verified that adding some buffer between tracking arcs did not affect the resulting formal uncertainties and, most importantly, the way the coupled solution compares to the decoupled one.
For each arc, we simulated both Doppler and range observables which are measurements, in the line of sight direction, of the spacecraft's position and velocity with respect to a ground station, respectively. Doppler observables were modelled with a noise level of 15 µm/s at an integration time of 60 s, while range observables have a noise level of 20 cm. This is quite precise but should actually be a conservative value, given the 1 cm range accuracy achieved by the BepiColumbo mission (e.g. Genova et al., 2021). For selected passes, as was done by Dirkx et al. (2017), we also simulated VLBI observables (lateral position of the target spacecraft) following methodology described by Pogrebenko et al. (2004) and Duev et al. (2016), with a noise level of 0.5 nrad. Doppler data were generated as unbiased, while we included arc-wise biases for both the range and VLBI observables. It should be noted that range and Doppler data are obtained in a topocentric frame, while VLBI observations are measured in the ICRF. For both the flyby and orbital phase, observations are subject to constraints on ground station visibility (occultation, Sun angle).
When presenting and discussing our results in the rest of this paper, the estimated states will generally be expressed in the RTN frame: the x-axis points from the central body towards the spacecraft or moon, the z-axis is aligned with the normal to the orbital plane and the yaxis completes the reference frame. In the following, they are referred to as the radial, tangential and normal directions, respectively.
In our simulations, we estimated the following set of parameters: • arc-wise JUICE initial states x • global or arc-wise Galilean moons' initial states x j (t 0 /t i ) (j = 1..4). The a priori uncertainty in position was set to 15 km in the three RTN directions. For the velocity components, we used the differences between the latest IMCCE and JPL ephemerides (NOE-5-2021 5 and JUP365 6 , respectively) as a con-5 https://ftp.imcce.fr/pub/ephem/satel/NOE/JUPITER/ 6 https://ssd.jpl.nasa.gov/sats/ephem/ servative a priori. These a apriori values are provided in Table 1.
• gravitational parameters of Galilean moons µ j (j = 1..4), using the a priori uncertainties provided in Schubert et al. (2004) and reported in Table 2.
• gravity field coefficients of Galilean moons C lm , up to degree and order 2, 4, 12 (6 when considering the flyby phase only) and 6, for Io, Europa, Ganymede and Callisto respectively. As a priori constraints, we used the formal uncertainties by Schubert et al. (2004) forC 20 andC 22 , which are given in Table 2. We applied Kaula's rule with K = 10 −5 for the remaining gravity field coefficients (σ = K/l 2 , Kaula, 1966).
• arc-wise accelerometer bias calibration factors c i , with the a priori constraint set to 10 −7 m·s −2 (10 times larger than in Cappuccio et al., 2020).
• arc-wise biases for range observables, with an a priori uncertainty fixed to 0.25 m.
• arc-wise biases for VLBI observables. We set the bias constraint at 0.5 nrad in both right ascension and declination (Charlot et al., 2020). Table 3 specifies whether a parameter is to be estimated globally or in an arc-wise manner. It highlights important differences between the two estimation methods, but also between the two steps of the decoupled approach. It must be stressed that the moons' gravity field coefficients are only included in the second step of the decoupled approach to account for the influence of uncertainties in the moons' gravity fields on the propagated state solutions (see Eq. 10). This is merely a way to avoid obtaining too optimistic formal errors because part of the uncertainties sources would be omitted. The a priori values for these coefficients are directly taken from the formal errors obtained after the first estimation step, and the gravity field solutions are actually not improved further by the second step, compared to these a prioris.
As shown by Dirkx et al. (2016Dirkx et al. ( , 2017, the influence of Jupiter's state and gravity field uncertainties on the estimation results was considered negligible in the post-Juno era , and these parameters were therefore not determined in our simulations. Tidal dissipation parameters were also excluded from the list of parameters to estimate in this preliminary study, keeping the focus of our analysis primarily on state estimation methods and on the resulting solutions for both the spacecraft and the moons.

Results
This section presents the results of our comparative covariance analyses, performed with both the coupled and decoupled estimation models (Sections 2.2 and 2.3). We   It should first be highlighted that the estimation problem is very close to being ill-posed, with extremely high condition number for the normal equations (Eq. 9). The exact values of the formal errors provided in the coming section should thus be treated cautiously. However, it must be noted that, as a verification, we also performed a deterministic least-squares estimation for the coupled case, to bring confidence in the formal uncertainties level (see Appendix B). It proves that our implementation of the lesser documented coupled model is correct, and that the obtained formal errors would be representative of the true errors under the assumptions of a covariance analysis (perfect dynamical and observational models, to be further discussed in Section 5). Our results therefore remain insightful, especially since we focus on comparing two estimation strategies (and not on absolute error values).
Given the near ill-posedness of the Galilean moons' state estimation problem, it is worth stressing that the condition number is higher when using the decoupled method, which is an important disadvantage of this approach. When reconstructing the moons' long-term dy-namics from the normal points, extremely high correlations between the position and velocity components in the radial and tangential directions even made the estimation problem non-invertible at first. Eventually, only the normal points' positions were therefore added as observables in the second step of the decoupled method (Section 2.3.1), to partially eliminate these correlations. In most other analyses, the normal points also include the central moon's position only (e.g. Durante et al., 2019;Di Ruscio et al., 2020;Di Ruscio, 2021)

Flybys phase only
This section presents the covariance analysis results obtained from observations simulated over the JUICE flyby phase only. We first discuss and compare the resulting formal errors in Galilean moons' states (Section 4.1.1), and in gravity field coefficients (Section 4.1.2). The sensitivity of the estimation solutions to the a priori covariances for the moons' initial states is then investigated in Section 4.1.3.

State estimation
Focusing on the first step of the decoupled estimation strategy, the uncertainties of the normal points generated for each flyby are displayed in Figure 4. For each moon, the reduction in the normal points' uncertainties as the number of flybys increases is clear. This is a direct consequence of updating the a priori knowledge for the moons' states, ensuring that each arc benefits from the previous ones (Section 2.3.3). Figure 4 shows that the position of the central moon is much better determined in the radial and tangential directions than it is in the normal direction (i.e. out-of-plane). Interestingly, the errors in radial and tangential positions are of similar orders of magnitude, especially for Callisto. This is due to high correlations between the tangential and radial state components. As an indication, Figure 5 shows the absolute correlations obtained when generating the normal points for the first flybys at Europa, Ganymede and Callisto. Callisto's first normal point is much less correlated than Europa's and Ganymede's. Looking at Figure 4, it appears that this does not indicate a good normal point determination, but, on the contrary, is due to the fact that the first, relatively high altitude flyby performed at Callisto (see Figure 1) does not allow the estimation to significantly improve the normal point determination compared to a priori values. The estimated state components for this normal point remain rather uncorrelated (since the default a priori covariances assume no correlation between parameters), but the associated formal uncertainties are quite large. The other flybys performed around Callisto, supported by the adopted a priori update strategy (see Section 2.3.3), progressively improve the quality of the normal points determination (see Figure 4), but also yield much higher normal points' correlations, which then become comparable to ones displayed for Europa and Ganymede in Figure 5. The global solutions for the Galilean moons' dynamics, either reconstructed from the normal points shown in Figure 4 in the decoupled case, or directly outputted by the coupled model, are displayed in Figure 6. The 1σ uncertainties in the moons' states, estimated at the beginning of the flyby phase from all flybys' data, were propagated through the 2030-2038 time period following the methodology presented in Section 2.1.2. The local uncertainty reductions in the propagated solutions clearly indicate when the flybys are performed.
Comparing both solutions in Figure 6, the coupled method leads to lower formal errors in the moons' radial positions (one or two orders of magnitude lower than in the decoupled case). On the other hand, the uncertainties in the moons' tangential positions are comparable between the two estimation strategies, and can even be locally slightly lower with the decoupled approach. These results follow from the similar formal error levels in the normal points' radial and tangential positions (see Figure 4), which translates into uncertainties of comparable orders of magnitude in both directions when reconstructing the global ephemerides solution. On the contrary, the coupled approach is able to more efficiently decorrelate the moon's radial and along-track motion.
Finally, differences between the coupled and decoupled solutions are not so significant in the normal direction and seem more arbitrary: the decoupled method performs slightly better for Io and Ganymede, while the converse is true for Europa and Callisto. The main difference between the two estimation approaches originates from the decoupled strategy only accounting for the dynamical coupling between the moons in the second step, when trying to reconstruct the moons' dynamics using the kinematic information contained in the normal points. The coupled model, on the other hand, includes all dynamical effects at once. As the strong dynamical coupling between the moons mostly manifests itself in the moons' orbital plane, the results obtained in the normal (i.e. out-of-plane) direction are less sensitive to the choice of estimation method.
While tidal dissipation parameters were intentionally excluded from our comparative state estimation analysis, preliminary insights can still be extrapolated from the expected ephemerides quality. Especially, an accurate determination of the moons' along-track positions is crucial to investigate tidal dissipation effects, through the secular change in mean motion that they induce. While the decoupled estimation led to slightly lower formal errors for the tangential positions of the moons (see Figure 6), we should keep in mind that the uncertainties obtained with the coupled model are thought to be more statistically representative, as mentioned in Section 2.4. It should therefore be highlighted that the determination of the moons' states in the tangential direction might be too optimistic in the decoupled case, possibly translating into lower formal errors for the tidal parameters. The coupled strategy, on the other hand, may prove beneficial to achieve realistic errors when trying to estimate Jupiter's dissipation at the forcing frequencies of the moons. This would be particularly important to further investigate whether Callisto is caught in a tidal resonance lock (Fuller et al., 2016;Lainey et al., 2020). It is however essential to stress that, with both models, achieving this tangential uncertainty level will be most challenging, and the results may well be limited by dynamical model error, as further discussed in Section 5.2.

Gravity field estimation
The formal errors for the moons' gravity field coefficients obtained with the decoupled and coupled approaches are provided in Figure 7, superimposed with their a priori values (Section 3.3). The limited number of flybys at Europa (2) and Ganymede (7) does not allow the estimation to significantly improve the gravity field solution with respect to the a priori knowledge. This is true for both the coupled and decoupled models, and neither of them seems to perform systematically and distinctly bet- ter for these two moons: for the rare coefficients for which an improvement is noticed compared to the a priori constraints, the formal errors are sometimes lower with one method, sometimes with the other. Results are different in Callisto's case: the 21 flybys, performed at lower altitudes and condensed over a short period of time (Figure 1), help to determine the gravity field well beyond a priori values. It is interesting to note that our formal uncertainties for Callisto's gravity field are comparable with those obtained in Di . As shown in Figure 7, the gravity field uncertainties given by the decoupled solution are slightly larger than those achieved with the coupled model. For some of Callisto's flybys, the decoupled method indeed performs poorly in decorrelating the moon's and spacecraft's arcwise states, which degrades the gravity field estimation. We must nonetheless stress that these differences remain small, the formal errors obtained with the two estimation strategies still being of comparable orders of magnitude.
The decoupled model was actually used in several gravity field estimation studies, mainly to avoid the challenges arising from the reconstruction of dynamically consistent ephemerides over long timescales, as the dynamical model fidelity could not be brought to the required level (e.g. Durante et al., 2019). Our results verified that the decoupled errors for the moons' gravity coefficients are also not too optimistic compared to the coupled solution. This analysis thus confirmed that opting for the normal points strategy for the JUICE flyby phase would not notably affect the gravity fields solution. This is particularly relevant for Callisto, while Ganymede's and Europa's gravity field determination will also significantly benefit from JUICE's orbital phase and the Europa Clipper mission, respectively.

Sensitivity to a priori knowledge
The decoupled solution was found to strongly depend on the a priori constraint applied to the moons' states before each flyby. As an experiment, we discarded the update strategy presented in Section 2.3.3 and applied the same default a priori covariances to all arcs. The normal points approach then led to rather different results, reported in Table 4. All position uncertainties get significantly larger when the state knowledge is not conveyed from one arc to the next. Results are the worst for the moons' radial and tangential positions, with errors increased by more than one order of magnitude. Updating the a priori information after each flyby is thus critical if realistic uncertainties are to be achieved with a decoupled approach. In particular, it progressively helps to decorrelate the central moon's and spacecraft's arc-wise dynamics.
When computing the observations' contribution to the solution using Equation 28, the average c q value for the moons' positions drops from c q >0.98 when using updated a priori covariances to ≈ 0.40 with the default ones (except for the first flyby of each moon for which no updated a priori is available and c q thus remains close to 1). This confirms that the a priori information then becomes predominant and significantly helps the solution.
On the contrary, the coupled solution is not noticeably affected by the adopted a priori values for the moons' states (c q ≈ 1), and thus appears significantly more robust. It also relies on a more straightforward, update-free strategy as it only uses the default a priori values (see Section 3.3).
It must be noted that the a priori knowledge for the moons' states, while driving the quality of the decoupled state estimation, has no notable impact on the gravity field solution, irrespective of the selected estimation method. This again shows that the main drawbacks of the normal points strategy do not significantly influence the estimated gravity fields. It confirms that the decoupled method is a good alternative when focusing on gravity field determination, in agreement with conclusions drawn in Section 4.1.2.

Flyby and orbital phases combined
We extended the tracking data set to include the orbital phase at Ganymede in addition to the flybys, again  . The gravity coefficients are plotted along the x-axis: for both the cosine C im and sine S im coefficients, they are grouped by degree and plotted by increasing order m. The label C/S im indicates where the coefficients of degree i start, and the order m of these coefficients progressively increases until the start of the next group of coefficients (degree i + 1).

State estimation
To apply the decoupled estimation to the orbital phase, we generated one normal point per tracking arc and determined the central moon's arc-wise state and the associated covariances at the centre of the arc. The propagated uncertainties in the Galilean moons' states are shown in Figure  8, for both the decoupled and coupled solutions. Callisto's state being almost only constrained by the flyby phase, the formal errors for this moon are similar to those discussed in Section 4.1, when excluding the orbital phase, and are therefore not discussed in the following (see Figures 6 and   8).
It should be stressed that Ganymede's formal errors fall below 1 m during JUICE orbit, in both the decoupled and coupled cases ( Figure 8). However, current dynamical models are likely far from being accurate enough to represent sub-meter level effects. Therefore, achieving the presented level of errors in reality will require these models to be rigorously adapted and validated. The implications of these modelling limitations, which differ for the coupled and decoupled solutions, will be further discussed in Section 5.
Comparing Figures 6 and 8 directly highlights the ephemerides improvement provided by the orbital phase. For both the coupled and decoupled solutions, the decrease in Ganymede's position uncertainties during JUICE orbit is clear. These results also clearly illustrate the strong dynamical coupling between the three innermost Galilean moons: the errors reduction for Io and Europa with respect to Figure 4.1.1 is indeed achieved by collecting more observations close to Ganymede. Figure 8 also confirms the flybys-based conclusions discussed in Section 4.1.1. In particular, the coupled estimation still provides a noticeable improvement in the radial direction, while errors in the moons' tangential positions are overall lower with the decoupled method. It is however interesting to note that these trends are accentuated for some moons, and attenuated for others. For Ganymede, the tangential position uncertainties are noticeably lower with the decoupled model during the orbital phase, and the error reduction in the radial direction achieved by the coupled solution remains limited. The opposite is observed for Io and Europa: the coupled solution's radial position uncertainties are on average one to two orders of magnitude lower than in the decoupled case, while the errors level remains comparable in the tangential direction. This is caused by differences in how each method captures the strong dynamical coupling between Io, Europa, and Ganymede. In the decoupled approach, the radiometric data collected during JUICE orbit are used to generate normal points solely at Ganymede. At first, these observations thus exclusively improve our knowledge of Ganymede's local states. The coupling between the three moons is only introduced in the second step of the decoupled estimation (Section 2.2.1), and Io's and Europa's solutions therefore benefit from the orbital phase in an indirect way, through very accurate normal points generated at Ganymede. On the contrary, the coupled model directly uses all data to estimate the four moons' states concurrently, and provides the most statistically accurate mapping of data uncertainty covariance to parameter covariance (see Section 2.4). In the coupled case, the solution improvement provided by the orbital phase is thus more evenly spread between the three innermost moons.

Gravity field estimation
The results and conclusions regarding the moons' estimated gravity fields are similar to those discussed in Section 4.1.2. For Europa's and Callisto's estimated gravity coefficients, there is actually no noticeable difference with respect to the flybys phase's results, which was expected since the orbital phase at Ganymede does not constrain other moons' gravity fields. The solution for Ganymede is however significantly improved by the orbital phase, as shown in Figure 9. It should be noted that these results rely on a simplified estimation setup, and that gravity field studies based on 3GM data from JUICE's orbital phase estimate Ganymede's gravity along with the moon's rotational parameters, Love numbers, etc. (e.g. Cappuccio et al., 2020). Nonetheless, the order of magnitude of the formal uncertainties reported in Figure 9 are in agreement with those obtained in dedicated 3GM studies (Cappuccio et al., 2020;De Marchi et al., 2021).
In our simulations, limited differences between the coupled and decoupled cases could still be detected from the flybys-based results, at least for Callisto whose gravity coefficients could be estimated beyond their a priori values (see Figure 7). However, such discrepancies between the two models are not observed anymore (for any moon) once the orbital phase is included, which is why Figure  9 only displays the coupled solution. Compared to the flybys, JUICE's orbital phase generates large amount of data, continuously collected over a longer period of time (as opposed to discrete arcs at each flyby). The contribution of the orbital phase thus completely dominates the gravity field solution (see Figure 9). The longer tracking arcs (i.e. 2 days, with 8 hours of tracking per day, instead of 8 hours only for the flybys) and the much larger numbers of observations allow both methods to properly decorrelate the spacecraft's and moon's dynamics, which explains why nearly identical uncertainties are obtained for Ganymede's gravity field in each case. This confirms conclusions drawn in Section 4.1.2, according to which the adopted state estimation strategy does not significantly influence the gravity field solutions.

Sensitivity to tracking settings
We re-ran our simulations with varying arc duration for the orbital phase, to investigate the sensitivity of each estimation method to the tracking configuration. As mentioned in Section 3.3, three test cases were considered with arcs of one day, two days (nominal) and one week, respectively. When shortening the tracking arcs from two days to one, the errors grow larger for both the decoupled and coupled solutions. The opposite is true when increasing the arc duration to a full week, as longer arcs allow a better decoupling of JUICE's and Ganymede's dynamics. As expected, this effect is the smallest for Callisto's uncertainties.
Importantly, the relative evolution of the decoupled and coupled solutions as the arc duration varies provides valuable insights into the fundamental differences between the two state estimation approaches. Figure 10 shows ratios in formal errors resulting from different arc durations. Except for Io's and Ganymede's normal positions, the improvement provided by longer arcs is systematically weaker in the decoupled case. This is caused by decoupling the spacecraft's and moons' state estimations. As already mentioned, lengthening the tracking arcs generally helps to decorrelate the spacecraft's motion from the central moon's, but does not necessarily reduce the correlations between the moon's own state components. This still results in lower state uncertainties for the normal points. However, it does not automatically reduce the correlations level in the second step of the decoupled estimation, the spacecraft's states being excluded anyway (see Section 3.3, Table 3). In the decoupled case, the reduction in the moons' state uncertainties thus remains limited by the normal points' high correlations, especially in the radial and tangential directions. On the contrary, the correlations between the spacecraft's and moon's states are directly included in the coupled model, since their dynamics are concurrently estimated. The coupled solution therefore fully benefits from longer tracking arcs, explaining why a more significant improvement is achieved for the radial and tangential positions than what the decoupled model allows.

Discussion: main strengths and challenges of both methods
This paper assesses the relative performance of two state estimation strategies, applied to the JUICE test case. To put the obtained formal errors and correlations into context, we discuss practical considerations related to data analysis in the following, and how they affect the coupled and decoupled models in different ways. Data processing challenges are detailed in Section 5.1, while model-related issues will be addressed in Section 5.2. As a possible mitigation strategy, we finally suggest a possible alternative state estimation approach in Section 5.3.  Figure 7. The gravity coefficients are plotted along the x-axis: they are grouped by degree and plotted by increasing order m, with the cosine coefficients first, followed by the sine terms. The label im indicates where the cosine coefficients of degree i start, and the order m of these coefficients progressively increases. They are followed by the sine coefficients of degree i (ordered similarly), until the start of the next group of coefficients (degree i + 1).  Figure 10: Ratios of the formal errors in position for different arc durations over JUICE's orbital phase (in the RTN frame: radial, tangential, normal). The results obtained with both the coupled and decoupled methods are displayed for the four Galilean moons (with no noticeable differences for Callisto). For each method and each moon, we provide the ratio of the formal uncertainties obtained with 1 day-arcs over 2-days arcs, and 1 day-arcs over 1 week-arcs.

Data processing considerations
One of the major differences between the coupled and decoupled techniques lies in the way the available data are processed. The decoupled method can theoretically treat each mission, and each mission phase, independently and generate as many normal points as required by the mission design. These normal points can later be combined with those determined from other missions and/or with different observation sets (e.g. optical, astrometric data). For the coupled solution, however, all data need to be processed at once.
For the Galilean moons test case in particular, reaching a global solution would ideally require to include spacecraft data from various missions, as well as Earth-based optical astrometric observations. More precisely, JUICE's imbalanced data set with its strong focus on Ganymede would be efficiently complemented by the Europa Clipper and Juno missions, with over 40 flybys at Europa planned for the former (e.g. Verma and Margot, 2018;Tarzi et al., 2019), and an orbital phase at Jupiter lasting since 2016 for the latter, which might be combined with a few crucial flybys at Io during the extended mission phase. In addition to radiometric data, spacecraft-based optical observations captured with navigation cameras can also be useful to help constraining the ephemeris solution. These can be direct imaging of other moons (e.g. JANUS data for JUICE, Dirkx et al. 2017), as well as eclipses and Sun or stellar occultations observed from the spacecraft (Andreoli and Zannoni, 2018).
Additionally, Earth-based photo-and astrometric observations of the Galilean moons have been collected over centuries. They include absolute and differential astrom-etry, already performed since the 17th century (e.g. data starting from late 19th century used in Lainey et al., 2009), as well as eclipses and occultations (e.g. Arlot and Emelyanov, 2019). More recently, mutual approximations were identified as interesting observables (Morgado et al., 2019b;Fayolle et al., 2021) and the first stellar occultation by the moon Europa was observed in 2017, with a remarkable accuracy of 0.80 mas (i.e. 2.55 km at Jupiter's distance) (Morgado et al., 2019a). Interestingly, the GAIA catalogue will facilitate the observations of such occultations in the future.
Merging all above-mentioned data sets and processing them in a single step does not only make the coupled estimation process slower, but also substantially increases its complexity. It first implies to carefully weigh all different observations to obtain a statistically balanced and realistic solution. The JUICE radiometric data, which led to formal state uncertainties below the meter level for the Galilean moons (see Figure 8), would indeed need to be properly combined with Earth-based optical astrometric data whose current accuracy remains larger than 1km (even stellar occultations). Furthermore, optical astrometric and radiometric tracking data are typically processed by different estimation tools, while the coupled estimation model would require a single software able to handle and process all observation types concurrently, which imposes significant practical constraints (many different observation types required, with, for each of them, suitable error and dynamical models, etc.).
Moreover, Earth-based observations are collected over longer periods of time compared to the spacecraft data, which are condensed over the planetary missions' timeline. This requires dynamical models to be consistent over both short and long timescales, and represents a major challenge for the reconstruction of a coupled solution, which will be further discussed below.

Modelling-related considerations
For the JUICE mission specifically, our covariance analysis indicates that Ganymede's formal state uncertainties get lower than the meter level during the orbital phase (see Figure 8). For such position errors to be meaningful, major modelling efforts would however be essential. Model-related issues are therefore expected to occur when real JUICE data become available, but their influence on the estimation can unfortunately not be easily quantified in a simulation analysis. Dirkx et al. (2016) analysed which effects would likely be detectable from a long-term (several years) signature in the dynamics. However, they did not address the observability and relevance of short-term periodic variations. Such potentially mismodelled dynamical effects (both long-and short-term) are however crucial to discuss here, as they would have different impacts on the coupled and decoupled solutions.
Among possible sources of inaccuracies, the models currently used for the dissipation occurring inside the moons rely on simplified analytical formulations (Lainey et al., 2009;Lari, 2018). These approximations circumvent the need for perfectly consistent rotational and translational models, still unavailable for natural satellites, but are based on time-averaging assumptions. The dissipation inside the moons however directly influences the orbital evolution of the moons themselves and plays a key role in the long-term dynamical history of the Jovian system, such that its accurate modelling is crucial.
Additionally, due to the coupling between the moons' translational and rotational motions, issues can also originate from mismodelled librations. Ganymede's libration will likely be observable in the dynamics of JUICE itself during the orbital phase (Cappuccio et al., 2020). The coupled model thus presents an interesting oppurtunity, since the libration's signature on the spacecraft and moon dynamics is different but would need to be fitted concurrently, possibly yielding a better constrained solution. On the other hand, modelling inconsistencies in Ganymede's librational motion would more critically degrade the residuals of the coupled solution.
For JUICE specifically, properly modelling the moons' rotations would furthermore require to reconcile what the different instruments are sensitive to. JANUS and GALA (navigation camera and altimeter, respectively) will observe the rotational motion of the moon's surface shell, which might be decoupled from the full body inertial rotation sensed by the radiometric data. Additionally, temporal variations in the central planet's gravity field are also not yet perfectly understood, and are suspected to be responsible for small, unmodelled time-dependent accelerations detected at periapsis in Juno  and Cassini data . It should nonetheless be stressed that an imperfect dynamical model can still achieve the required accuracy if properly parametrised (e.g. if time-dependent librational and gravitational variations are adequately adjusted), provided the relevant free parameters are incorporated in the estimation, and sufficiently decorrelate from the other parameters.
The fundamental differences between the coupled and decoupled approaches (see Figure 2) significantly influence how the above-mentioned modelling issues might affect the estimation solutions. Compared to the coupled case, the decoupled approach indeed estimates more state parameters, determined more locally (Table 3), which directly increases the ability of this model to absorb dynamical modelling inaccuracies.
Additionally, the long-term moons' dynamics are only reconstructed in the second step of the decoupled estimation, and they are adjusted to the normal points (i.e. arcwise covariances of the moons' states). For the JUICE flyby phase for example, the decoupled model is thus fitting formal state uncertainties generally ranging from tens up to hundreds of meters (the last flybys at Callisto getting closer to the meter level for the radial and tangential positions). The coupled approach, on the other hand, is directly adjusting the parameters to the radiometric data, with expected accuracies of 20 cm and 15 µm/s for the range and range-rate of JUICE with respect to the Earth (Section 3.3). It is therefore significantly easier to obtain flat residuals with the decoupled estimation strategy (i.e. zero-mean residuals, noise within expected observation errors), as observables with accuracies up to 10 1 − 10 2 m bring sufficient flexibility to (at least partially) absorb dynamical modelling inaccuracies.
The fact that the coupled method estimates the moons' dynamics in a single arc also drastically reduces its ability to absorb such model errors. It indeed implies that the estimation model cannot compensate for modelling inaccuracies by a local moon's state variation without affecting other arcs, and possibly conflicting with their own observational constraints. Furthermore, the decoupled estimation leaves out the arc-wise spacecraft states when determining the moons' long-term dynamics (see Section 2.3.1). The coupled model, on the contrary, imposes perfect consistency between the spacecraft's and moons' state estimation solutions, which further reduces its degrees of freedom.
In practice, the above entails that the concurrent estimation of the spacecraft's and moons' states would need accurate and consistent dynamical models over both short and long timescales. The current modelling fidelity has for example not yet allowed a coupled solution to be achieved from Cassini and Juno data (e.g. Durante et al., 2019). Furthermore, even if a global solution can still be reconstructed for the natural bodies' dynamics, modelling errors are expected to manifest themselves in high, incompressible residuals, due to the coupled strategy's lack of flexibility. On the contrary, in the decoupled case, modelling issues are more likely to be absorbed in the final estimated states and to thus remain unnoticed. Reflecting back on the normal points obtained for JUICE, with errors largely below one meter, more realistic uncertainties in agreement with the available dynamical models would significantly raise this error level. It would therefore degrade the quality of the decoupled estimation, but not prevent the obtention of a viable solution. On the other hand, the (large) post-fit residuals obtained with the coupled estimation can be indicative of the magnitude of the dynamical modelling inaccuracies, and help to interpret the decoupled solution's true uncertainty.

Possible alternative strategy
Despite the promising results obtained with the coupled model, Sections 5.1 and 5.2 highlighted crucial challenges that would need to be addressed before a reliable and statistically consistent coupled solution can be reconstructed for the Galilean moons' dynamics. Hybrid approaches, halfway between fully coupled and decoupled strategies, could therefore prove promising. In such hybrid scenarios, the coupled model can be applied more locally, rather than on the entire time period of interest. For the JUICE test case, a coupled solution might typically be attainable over the GEO/GCO5000 and/or GCO500 phases. Some flybys could also be processed concurrently (e.g. the two flybys at Europa). It would also be possible to reconstruct a coupled solution over the entire JUICE mission, but to combine it with other mission/data sets in a separate step, to mitigate the data merging issues highlighted in Section 5.1.
This would efficiently mitigate the effects of long-term modelling inconsistencies and thus make a coupled solution achievable locally. Such local coupled solutions can then be treated as normal points or a priori information for other analyses, and combined with those generated for different missions or with optical astrometric data, as discussed in Section 5.1. This would also eliminate the need to process astrometric and radiometric data in a single estimation tool as in the fully coupled strategy (see Section 5.1), by splitting the data analysis steps. This hybrid approach could represent an interesting alternative: it would potentially allow the estimation solution to benefit from a higher sensitivity to the system's full dynamical coupling, while still guaranteeing that a viable solution can be achieved.

Conclusions
We provided the complete formulation for a coupled, concurrent state estimation of the spacecraft and natural bodies from planetary missions' data. Such a coupled model has already been used in past studies (i.e. Dirkx et al., 2018;Lari and Milani, 2019), but, to the best of our knowledge, was not explicitly described in the literature. We then performed a detailed covariance analysis comparing the decoupled and coupled estimation strategies for the upcoming JUICE mission. The realism of the formal errors given by the coupled model was verified by running a deterministic simulation and comparing the least-squares estimation errors with the formal uncertainties.
The JUICE mission will make us face both unique opportunities and challenges due to the unusual mission profile and unprecedented accuracy of the radiometric data, used to reconstruct the strongly coupled dynamics of the Galilean moons. Our study primarily assessed how a coupled solution, if attainable, would affect the accuracy of the Galilean moons' ephemerides. The results of the decoupled estimation approach, on the other hand, indicate the uncertainty level that would be achievable in case the coupled model failed to reconstruct a viable solution for the moons' dynamics (see discussion in Section 5). It must be stressed that we do not anticipate the conclusions of our study to depend on the CReMA (i.e. trajectory) version used for the JUICE spacecraft. While the absolute uncertainty levels for the estimated parameters might vary, the comparison between the two state estimation methods is expected to yield similar results.
We first showed that selecting appropriate a priori values for the moons' states is critical for the decoupled solution. The state solution improvement must be accounted for by updating the a priori covariance from one normal point to the next. Discarding the information gained in previous arcs indeed leads to poorly constrained normal points (limited to kilometer level accuracy). This eventually drives the decoupled solution to be one or two orders of magnitude less accurate than the coupled one, for all moons and in all directions (e.g. tens of meters against ∼ 500 m for Ganymede's position uncertainties with and without a priori update, respectively, see Table 4).
Furthermore, the flybys-based results already highlighted notable differences between the coupled and decoupled state solutions. The coupled method achieves lower radial position uncertainties (∼ 10 2 m for Io and Europa, ∼ 10 1 m for Ganymede, ∼ 10 0 m for Callisto), which are one or two orders of magnitude smaller than in the decoupled case. In the tangential direction, the formal errors given by the decoupled approach are however slightly lower: ∼ 10 m against ∼ 30 m, keeping Ganymede's initial position uncertainties as an example. As discussed in Section 2.4, the coupled model is nonetheless considered to provide a statistically accurate representation of the estimation problem, which might imply that the tangential position uncertainties obtained with the decoupled model are slightly optimistic.
Our results also proved that JUICE orbital phase is crucial for the quality of the ephemerides of the four Galilean moons. Assuming a perfect dynamical model, Ganymede's position uncertainties would even reach submetric levels during the orbit. Interestingly, it seems that Io's and Europa's solutions benefit more from the orbital phase when adopting a coupled approach (Europa's radial position also determined with an accuracy better than one meter, see Figure 4.2.1), while the improvement appears stronger for Ganymede in the decoupled case. Including the orbital phase in our simulations actually enhances the differences between the two estimation strategies, and indicates that the dynamical coupling between the Galilean moons is better captured by the coupled model.
Concerning dynamical parameters, the adopted state estimation method is not so influential for the gravity field coefficients' estimates. This indicates that the normal points strategy is perfectly adapted for gravity field determination from JUICE radiometric data (e.g. Magnanini, 2021), and generally for most local parameters studies as well. Rotational and tidal dissipation parameters, on the other hand, were omitted from our simulated estimations. Future studies should include these parameters and try determining the dissipation in Jupiter at the moons' frequency (Lainey et al., 2009(Lainey et al., , 2020. It would then be crucial to investigate whether the decoupled model can achieve realistic uncertainties despite a slightly optimistic estimation of the moons' along-track motion.
As a preliminary step towards an improved solution for the Galilean moons' ephemerides, we limited ourselves to an analysis based on simulated data only. As extensively discussed in Section 5, many issues are however expected to arise when real JUICE radiometric data become available and are fed to the estimation process, as happened with Cassini and Juno data .
For the JUICE test case, the coupled estimation model nonetheless appears to yield promising improvements with respect to the decoupled strategy: it overall reconstructs a more balanced solution for the three Galilean moons in resonance, with significantly reduced uncertainties in the moons' radial positions. Despite the challenging issues that may emerge, this motivates future efforts to achieve a coupled state estimation for the Galilean moons. The coupled model being an extension of the classical formulations used by the decoupled estimation, the added value of performing a separate deterministic verification for the decoupled model appears limited. Additionally, implementing and running the entire least-squares estimation with the decoupled model would be more demanding, since it relies on two consecutive estimation steps ( Figure  2b) and in practice also requires multiple sequential estimations to determine all normal points, due to the chosen a priori update strategy (see Section 2.3.3).

Appendix B.2. Approach and settings
When performing an iterative weighted least-squares inversion, the variation in parameters values ∆q i , at each iteration i, is given by (e.g. Montenbruck and Gill, 2000) ∆q i = (B.1) where ∆z i is the vector containing the observations residuals at iteration i, and ∆q 0,i is the difference between the current parameters estimates and their a priori values. The rest of the notation follows that of Equation 9. In our simulations, true errors are directly defined as the difference between the true parameters values (i.e. values assumed in our simulated environment) and their estimated values. To keep our simulation as realistic as possible, the observations are modelled with the bias and noise levels given in Section 3.3: noise levels of 20 cm, 15 µm/s and 0.5 nrad for range, Doppler and VLBI observations, respectively, with biases of 1 m for range measurements and 0.5 nrad for VLBI data. We also applied small initial perturbations to the a priori values of the estimated parameters, to mimic imperfect knowledge of the true parameters values. The way the initial parameters' offsets were defined is detailed below: • moons' initial states: initial perturbation ∆q 0 calculated such that, once propagated over the duration of the JUICE mission, the difference between the perturbed and unperturbed moons' trajectories would be of the order of magnitude of the current ephemerides' accuracy (see Section 3.3). This led to initial offsets of 100 m in position and 5 mm/s in velocity.
• JUICE's arc-wise states: offset of 100 m in position and 0.01 m/s in velocity (in all inertial directions); • other parameters (gravitational parameters, gravity coefficients, accelerometer calibration biases, range and VLBI observation biases): perturbation set to 5% of their true values.
Some differences should be noted with respect to the nominal estimation settings used for our covariance analyses (Section 3.3). Io's initial state was removed from the list of estimated parameters, due to the absence of Io's flyby, resulting in a lack of observations at this moon. The least-squares inversion could indeed not converge when trying to estimate its state alongside with those of the other Galilean moons. This already provides valuable insights about the challenges one will face when conducting real data analysis (as further discussed in Section 5). It should be noted that this convergence issue is only caused by JUICE's imbalanced data set, and that many more difficulties are to be expected when dynamical and observational models will be confronted with real data (see dis-cussion in Section 5). It is also important to highlight that this was only observed when perturbing the a priori values of the estimated parameters. When no initial perturbation is applied, the least-squares estimation could converge while solving for Io's state along with the other moons'. All true errors were then found to be smaller than 3σ, σ designating the corresponding formal errors.
Furthermore, the elevation and occultation checks normally performed to verify the viability of an observation were switched off, thus assuming constant link during the tracking arcs. One specific tracking arc did not have enough viable observations for the estimation problem to converge otherwise. While this issue could have also been fixed with tighter a priori constraints, we adopted a simpler approach. We must stress that removing the visibility requirements, thus adding a few observations, does not lessen the relevance of our verification analysis, since it focuses first and foremost on the consistency between true and formal uncertainties provided by the coupled model, and not on absolute error values.
We conducted this iterative least-squares inversion for the flyby phase case only. This was mostly motivated by the high computational load required by such deterministic runs. We do expect similar results and conclusions when including the orbital phase at Ganymede.

Appendix B.3. Results
After 5 iterations, the least-squares inversion reached convergence, and the final observations residuals are provided in Figure B.1. As expected, the residuals follow a Gaussian distribution with almost zero mean, and a standard deviation close to the observations noise level (Section 3.3). This can most clearly be observed for the Doppler residuals, due to the larger number of observations available compared to range and VLBI data.
The histogram of the ratios between true and formal errors is displayed in Figure B.2, and all ratio values are between 0 and 9. Precisely quantifying which true to formal errors ratio should be expected is far from straightforward, as it is both data-and parameter-dependent. While true errors 2 to 3 times larger than formal errors can be typically expected for planetary ephemerides (Jones et al., 2015(Jones et al., , 2020, for Saturn's ephemeris derived from VLBA observations), larger true to formal errors ratios have been found for dynamical parameters estimated from spacecraft tracking data (gravity fields, rotational parameters, etc.) (e.g. Konopliv et al., 2011;Mazarico et al., 2015). These results were however all based on real observations, such that inaccurate dynamical or observations error models significantly contributed to the true/formal errors discrepancy.
In our simulations, however, our models are still deemed perfect, and true errors should be comparable to formal uncertainties if the estimation is able to converge towards the correct parameters values. In other JUICE simulation studies, true errors were indeed found to be similar or even slightly lower than formal ones for Callisto's gravity field in Di , although it is unclear if the parameters values were initially perturbed, and by how much. In Lari and Milani (2018), similar results were obtained, but true errors were averaged over 10 experiments with different initial parameters perturbations (smaller than those we applied).
In our simulation, for most parameters (97% of the 432 estimated parameters), the true error is smaller than 3σ (i.e. true to formal errors ratio smaller than 3). A significant fraction (62%) of the true errors are actually smaller than 1σ. These results confirm that, for the vast majority of estimated parameters, the coupled estimation model is able to properly converge towards the true parameter value, and that errors can be expected to be in agreement with 3σ uncertainties. All parameters exhibiting a true to formal errors ratio larger than 3 were actually found to be Europa-related properties (gravity field coefficients, initial position along the x-axis, etc.). Due to the limited amount of observations collected at Europa (only two flybys, see Figure 1), these parameters are either highly correlated (second degree and order gravity field coefficients with Europa's state), or cannot be estimated beyond their a priori constraints (c q close to 1 for Europa's other gravity coefficients, implying a so-called biased estimation). This might at least partially explain why unexpectedly large true to formal errors ratios were obtained for these parameters.
Overall, this deterministic simulation succeeded in verifying that the 3σ uncertainties provided by the coupled estimation are good indicators of the expectable estimation errors. This increased confidence in our covariance analyses. The effect of inaccuracies in dynamical or observational modelling was however not investigated (see discussion in Section 5). Interestingly, and despite still relying on ideal models, this simulation nonetheless highlighted issues arising when trying to complete the leastsquares estimation, due to the high condition number resulting from the lack of observations at Io and Europa.  Figure B.2: Distribution of the ratio between true and formal errors for all estimated parameters. A value smaller than 1 indicates that the true estimation error is lower than the associated formal uncertainty. 432 parameters were estimated in total.