Measuring Ganymede’s Librations with Laser Altimetry

Abstract

Jupiter’s moon Ganymede might be in possession of a subsurface ocean located between two ice layers. However, from Galileo data it is not possible to unambiguously infer the thickness and densities of the individual layers. The upcoming icy satellite mission JUICE (JUpiter ICy moons Explorer) will have the possibility to perform more detailed investigations of Ganymede’s interior structure with the radio science experiment 3GM and the GAnymede Laser Altimeter (GALA). Here we investigate the possibility to derive the rotational state of the outer ice shell by using topography measured by laser altimetry. We discuss two different methods to invert synthetic laser altimetry data. Method 1 is based on a spherical harmonics expansion and Method 2 solves for B-splines on a rectangular grid. While Method 1 has significant limitations due to the omission of high degrees of the global expansion, Method 2 leads to stable results allowing for an estimate of the in-orbit measurement accuracy. We estimate that GALA can measure the amplitude of Ganymede’s librations with an accuracy of 2.5–6.6 $\mathsf{\mu }$rad (6.6–17.4 m at the equator). This allows for determining the thickness of an elastic ice shell, if decoupled from the deeper interior by a subsurface ocean, to about an accuracy of 24–65 km.

1. Introduction

Three of the Galilean satellites, Io, Europa, and Ganymede, are locked in a 1:2:4 orbital resonance. This Laplace resonance forces the satellites into eccentric orbits making them subject to strong, time-variable, gravitational forces exerted by Jupiter. These forces trigger a response of the satellites dependent on their interior structure. One consequence is tidal heating, which allows the sustaining of subsurface oceans making Europa and Ganymede compelling targets for future exploration. Measuring the tidal deformation and resulting changes in the gravity field provides direct evidence of the subsurface oceans and can reveal structural and rheological characteristics of the outer ice shell [1,2]. In addition, the potential decoupling of the ice shell from the solid interior by an internal ocean would trigger a response on the rotational dynamics of the ice shell [3,4]. Since the long axis is pointed towards the empty focus of the moons’ orbit around Jupiter, librational torques arise, which might force the ice shell to librate about their mean rotation axis [5]. The amplitude of this physical libration depends on the interior structure of the respective satellite [4].

For Ganymede, our current understanding originates mainly from Galileo spacecraft data which determined the gravitational constant $GM$ and the second-degree gravity field [6]. However, since only near-equatorial flybys have been used, only the value for $C_{22}$ could be determined by radio science tracking. Assuming that Ganymede is in hydrostatic equilibrium, the Darwin–Radau relation [7] allows determination of the value for $J_2=10C_{22}/3$. These two gravity field parameters give the normalized polar moment of inertia $C/M{R}^2=0.3105$ [6] which indicates that Ganymede is a highly differentiated object. Together with the mean density of $\rho =1936$ kg/m${}^3$, the resulting internal structure as we understand it today is as follows: Ganymede is likely composed of an inner iron core responsible for the generation of the intrinsic magnetic field, surrounded by a silicate mantle [8]. A putative liquid ocean revealed by induced magnetic field measurements [9] would then be sandwiched between a high-pressure ice layer and the outer ice I layer [10]. Since the high-pressure ice layer might consist of ice III, IV, or VI, it is also conceivable that multiple oceans are found between the individual phases of high-pressure ices [11]. However, the measurement of the induced magnetic field does not allow constraining the thickness of the ocean or the outer ice shell. Also, the small density contrast between water and ice does not allow for a determination of the ice/water ratio.

Therefore, additional measurements of yet unknown geophysical parameters have been suggested. In particular, the combination of tidal Love numbers and the measurement of the rotation state can unambiguously confirm or disprove the existence of Ganymede’s subsurface ocean and constrain the ice shell thickness. Among the geodetic measurements, the most promising investigation is the measurement of the radial and gravitational Love numbers, $h_2$ and $k_2$, which is well suited to detect an ocean and may constrain the thickness and rheology of the overlaying ice cover [12,13]. The peak-to-peak amplitude of the radial tidal deformation is around 6–7 m at the equator if Ganymede possesses an interior ocean and only a few centimeters for an entirely solid satellite [12]. Consequently, the detection of the tidal response is one of the major scientific objectives targeted by the JUpiter Icy moons Explorer (JUICE) at Ganymede [14] as it is for the Europa Clipper mission at Europa [15]. The potential to measure tides has been investigated in previous studies for gravitational and radial tidal deformations in the case of Europa [16,17] and for Ganymede [13,18]. However, these studies also show that the tidal Love numbers $h_2$ and $k_2$ further depend on the ice shell rheology and to some extent on the deeper interior. Also, a linear combination of both Love numbers [2,13,19] still shows a wide range of ambiguities and is highly model dependent. Hence, the analysis of the libration amplitudes may impose additional independent constraints, as demonstrated recently by Thomas et al. [20] for Saturn’s moon Enceladus. If the shell is decoupled by a global ocean and if it is assumed to be rigid, then the amplitude on the surface can further reveal the ice shell thickness and density [21]. This comes with the caveat that models by Van Hoolst et al. [4], which include elastic effects and additional pressure and gravitational coupling of the surface and interior, show that the libration amplitude can be significantly decreased in the presence of tides. Also, in the case of Europa, a thin ice shell can lead to libration amplitudes heavily dependent on elasticity [22].

Further on, the obliquity can provide additional information about the ice thickness and coupling processes in the interior [3]. The detection of a difference in the libration and/or obliquity, when measured at the surface and within the gravity field, could suggest a differentially orientated solid interior and hence provide additional evidence for a global ocean.

The measurement of tides and rotational state is envisioned for the upcoming missions to the Jovian system, for Ganymede primarily by the JUICE spacecraft [14] which is expected to enter the moon’s orbit in 2033. The spacecraft carries ten instruments, one being the GAnymede Laser Altimeter (GALA) [23,24]. The most important mission phase for GALA will begin after the transition from an initial high elliptic into the circular GCO-500 orbit at 500 km altitude, which will be maintained for 132 days. During this phase GALA aims at completing its two main scientific goals: (1) obtaining a global topographic coverage of Ganymede and (2) detect its radial tidal deformations [25]. Further important objectives include the determination of the rotational state of Ganymede including rotation rate, librations, and obliquity. The aim of this work is to investigate to which accuracy the measurement of the rotational state and in particular Ganymede’s libration amplitudes as well as the tidal Love number $h_2$ can be used to constrain the interior structure of the satellite.

For this purpose, we investigate two independent non-linear inversion routines. One is applied to a spherical harmonic expansion of the global topography and the second one focuses on the inversion of cubic B-splines on a rectangular grid [26,27]. In the following we will call these two inversion strategies Method 1 and Method 2, respectively. We assessed to what accuracy the libration amplitudes and $h_2$ can be inferred by both methods, considering the current mission scenario and instrument as well as spacecraft performance. Section 2.1 and Section 2.2 give an overview of the state of knowledge of Ganymede’s geodetic and geophysical parameters derived from former missions and how they are implemented in this work. Whereas in Section 2.3 we explain how the synthetic measurements are created. Section 2.4 then explains in detail the two different inversion strategies investigated here. In Section 3 the results and possible errors are discussed. Finally, in Section 4 we discuss the implications of our findings for revealing the interior structure of Ganymede.

2. Methods

2.1. Model of Ganymede’s Rotation State

The rotation state of a planetary body is commonly described by three time-dependent angles, which are defining its orientation with respect to the International Celestial Reference Frame (ICRF). The angles are the right ascension $\alpha \left(t\right)$, describing the orientation of the equatorial plane with respect to the x-axis of the ICRF, the declination $\delta \left(t\right)$, which is the inclination of Ganymede’s equatorial plane towards the celestial equatorial plane, and the orientation of the prime meridian $\omega \left(t\right)$ (Figure 1) for a given time t. The time series of these three Euler angles are given by

$$\begin{array}{ccc}\hfill \alpha & =& {\alpha }_0+{\alpha }_{1}t+{\alpha }_2{t}^2\hfill \end{array}$$

(1)

$$\begin{array}{ccc}\hfill \delta & =& {\delta }_0+{\delta }_{1}t+{\delta }_2{t}^2\hfill \end{array}$$

(2)

$$\begin{array}{ccc}\hfill \omega & =& {\omega }_0+{\omega }_{1}t+{\omega }_{\mathrm{lib}}\left(t\right),\hfill \end{array}$$

(3)

where ${\alpha }_0$ and ${\delta }_0$ are the orientation of the rotation pole in right ascension and declination at $t=0$, while ${\alpha }_{1,2}$ and ${\delta }_{1,2}$ are their respective secular variation. ${\omega }_0$ denotes the orientation of the prime meridian at $t=0$ (prime meridian constant) and ${\omega }_1$ is the rotation rate of the body. The last term ${\omega }_{lib}$ includes the longitudinal physical libration terms under investigation

$${\omega }_{\mathrm{lib}}\left(t\right)=\sum _{j=1}^J{\gamma }_{j}sin({\omega }_{j}t+{\varphi }_j),$$

(4)

with the libration amplitudes ${\gamma }_j$, the libration frequencies ${\omega }_j$ and the phases ${\varphi }_j$.

The time evolution of Ganymede’s rotational state is described by the International Astronomical Union (IAU) report [28]. Additional periodic librations due to external torques by Jupiter change the rotation rate periodically. Besides the main libration forcing period of 7.16 days, resulting from Ganymede’s orbit about Jupiter, diverse additional libration periods emerging from perturbations of Ganymede’s orbit by other Galilean moons have been calculated by Rambaux et al. [21]. In addition to the 7.16 days period, we considered the 50.14 days period as well as the 482.06 period. Long periodic libration terms with periods of more than 5 years cannot be determined due to the comparatively short GCO-500 phase of 132 days and have been neglected.

The rotational elements currently adopted by the IAU Working Group for Cartographic Coordinates and Rotational Elements (WGCCRE) date back to the 1970s [29] and do not reflect the improved knowledge in the orbits of the satellites obtained in recent decades [30,31]. Nonetheless, to model Ganymede’s rotational state, a series expansion can be derived from ephemerides of the satellites under the assumption of a perfect spin-orbit resonance. To obtain the resonant rotation model we followed the procedure of Stark et al. [32]. We first obtained the osculating orbital elements of the satellite in the ICRF, applying a sampling time step of 12.5 min, corresponding to about 0.1% of Ganymede’s orbital period and using the full time interval available in the satellite ephemerides (200 years). Then we decomposed each orbital element in a time series composed of a secular trend and a sum of periodic terms by using a Fourier transformation. For the resonant rotation state we assumed that the satellite occupies a Cassini State 1 with zero obliquity. Furthermore, we assumed that the free precession period is much smaller than any periods in the orbit orientation variations. With these assumptions the satellite’s rotation axis is precisely following the instantaneous orientation of the orbital plane. This assumption is supported by the fact that amplitudes of very short variations of the orbital orientation are typically very small and can be neglected in practice. Finally, we linearized the precession and nutation terms for the time of GALA observations. As a consequence of this linearization we obtain large precession terms, which reflect the long-period precession of Ganymede’s orbital plane. We want to stress that the derived rotation parameters listed below are only valid in the timeframe from January 2033 until June 2033. Following these assumptions, the resonant rotation parameters for Ganymede are

$$\begin{array}{cc}\hfill \alpha \left(t\right)& ={268.682303}^{\circ }-{2.360616}^{\circ }t/\mathrm{cy},\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill \delta \left(t\right)& ={72.017968}^{\circ }-{46.494653}^{\circ }t/\mathrm{cy}+{70.706538}^{\circ }{(t/\mathrm{cy})}^2,\hfill \end{array}$$

(6)

where cy is a Julian century (36,525 days) and t is the time with respect to the J2000.0 epoch measured in days at the TDB time scale (ephemeris time). The prime meridian angle $\omega \left(t\right)$ is given by

$$\omega \left(t\right)={42.7011}^{\circ }+{50.31754547}^{\circ }t+{\omega }_{\mathrm{lib}}\left(t\right),$$

(7)

where the rotation rate parameter reflects the mean motion and the pericenter precession averaged over the time frame of GALA observations. The forcing terms of the longitudinal libration ${\omega }_{\mathrm{lib}}\left(t\right)$ we again follow the assumption of the perfect spin-orbit resonance and obtain these terms by measuring the angular separation between the line joining the centers of mass of Jupiter and Ganymede and the direction of the x-axis for a uniformly rotating Ganymede (with the rotation rate and rotation axis orientation as noted above). A frequency analysis of the angular separation gives

$$\begin{array}{cc}\hfill {\omega }_{\mathrm{lib}}\left(t\right)=& {0.237779}^{\circ }sin\left({360}^{\circ }t/7.156596\mathrm{days}+{42.083401}^{\circ }\right)\hfill \\ & +{0.070105}^{\circ }sin\left({360}^{\circ }t/7.051491\mathrm{days}+{225.924113}^{\circ }\right)\hfill \\ & +{0.018412}^{\circ }sin\left({360}^{\circ }t/6.262927\mathrm{days}+{66.710268}^{\circ }\right)\hfill \\ & +{0.005389}^{\circ }sin\left({360}^{\circ }t/12.476995\mathrm{days}+{107.957696}^{\circ }\right)\hfill \\ & +{0.000833}^{\circ }sin\left({360}^{\circ }t/50.141824\mathrm{days}+{75.978031}^{\circ }\right)\hfill \\ & +{0.010478}^{\circ }sin\left({360}^{\circ }t/482.056267\mathrm{days}+{14.639341}^{\circ }\right).\hfill \end{array}$$

(8)

In addition, we have included the forcing terms with periods of 50 and 483 days based on the work by Rambaux et al. [21]. Ganymede’s actual libration amplitudes, i.e., the response to the forcing, are however dependent from unknown interior properties. The response of Ganymede would correspond to the here given forcing amplitudes only in the unrealistic case of the free libration period being much smaller than the forcing periods, i.e., Ganymede exactly follows the evolution of the forcing angle. The actual libration amplitudes are certainly much smaller as Ganymede’s free libration period is likely to be between 20 to 140 days (see Section 4).

2.2. Synthetic Topography

Due to the lack of detailed topographic information of Ganymede we used a synthetic topography model for simulation purposes. The model has been parameterized in spherical harmonics and can therefore be expressed for a given latitude and longitude ($\mathsf{\Theta },\lambda$) as

$$r_{\mathrm{topo}}=\sum _{n=0}^{N_{\max }}\sum _{m=0}^n{P}_{nm}(sin\mathsf{\Theta })\left[{\overline{C}}_{nm}cos\left(m\lambda \right)+{\overline{S}}_{nm}sin\left(m\lambda \right)\right],$$

(9)

with $P_{nm}$ describing the associated Legendre Polynomials of degree n and order m. To adopt a realistic long wavelength topography we used the ellipsoid 1 of Zubarev et al. [33] for the degree two terms, with ${\overline{C}}_{00}=\mathrm{2,632,973}$ m, ${\overline{C}}_{20}=-713$ m, and ${\overline{C}}_{22}=292$ m. Higher order topography coefficients up to degree and order L are randomly generated (Figure 2) under the assumption that the topography follows a power law [34]

$${\sigma }_n^2=\frac{P^2}{n^b}=\sum _{m=-n}^n{\overline{C}}_{nm}^2.$$

(10)

Above, ${\sigma }_n^2$ is the variance of the topography at degree n and the parameters P and b describe the roughness and decay of signal power. The variance of a single topography spherical harmonic coefficient ${\overline{C}}_{nm}$ is given by

$${\sigma }_{nm}^2=\frac{{\sigma }_n^2}{2n+1}.$$

(11)

When using the spherical harmonics expansion, Method 1, it is sufficient to generate topography up to degree $L=200$. For Method 2, when inverting B-splines on a rectangular grid, we generated the topography up to $L=7999$. Overall we tested three different power laws with $P=800$ m and $b=1$, $P=1549.19$ m and $b=2$, as well as $P=2144.76$ m and $b=3$, respectively. The resulting topography for the version with $b=3$ is depicted in Figure 2. Power law models with $b=1$ give a rather rough topography, which means by consequence that we are left with more topography that is not modeled at the end of the analysis, increasing the uncertainty of the overall solution. For comparison, for Earth and Venus $b\approx 2$ [35]. Only limited topographic data on Ganymede is available to confidently state which power law applies. However, the power law with $b=1$ yields rather rough topography with amplitudes exceeding 20 km. Since such topographic variations are not observed on Ganymede, we considered $b=1$ only for Method 1 and discarded this option for Method 2.

As a non-rigid body, Ganymede is affected by tidal deformations. The radial deformation $u_r$ is usually parameterized by the dimensionless Love number $h_2$ and the tidal potential V

$$u_r=\frac{h_{2}V(\lambda ,\theta ,t)}{g_0},$$

(12)

where $g_0$ is the gravitational acceleration on the surface. The topography vector at a given position on the surface is then $r_s=r_{\mathrm{topo}}+u_r$.

2.3. Generation of Synthetic Laser Observations

To test the performance of the inversion for librational amplitudes, we generated synthetic laser altimetry data. We considered groundtracks from whenever JUICE is closer than 1300 km to Ganymede’s surface including dropouts due to downlink times (approximately 8 h/day considering ground station visibility). The considered data includes the elliptical orbit and the polar, circular orbit with 500 km altitude. This final phase, called GCO-500 will be the primary science phase for GALA and will be maintained for 132 days. For the type of measurement investigated here, a good temporal and spatial coverage is needed to sample multiple libration periods on the one hand and to obtain a dense net of groundtracks for a stable geodetic frame on the other hand. However, since for the latter the main limiting factor is the distance between individual tracks, a high along-track resolution is less important for this specific measurement. GALA’s nominal shot frequency is 30 Hz and leads to a ground separation of consecutive footprints of about 50 m in the GCO-500. However, the cross-track distance between neighboring profiles is much higher (in the order of a few km). To reduce the computational load in this simulation we used only every 10th GALA shot. We expect these results to improve when using the nominal shot frequency.

For the given trajectory, the spacecraft’s position within the ICRF with respect to Ganymede’s center of mass is the initial input for the forward model. Each laser observation has a time stamp which is referenced to the spacecraft clock (Figure 3). To simulate range measurement errors, we used the numerical performance model presented in [36]. GALA is operating with 17 mJ laser and a 25 cm (diameter) telescope. The return pulses are sampled on board with a frequency of 200 MHz and fitted to assess their centroid. This technique allows for a subsample range resolution of 10 cm. The ranging accuracy is in the order of 1 m, slightly dependent on the slope and the local albedo of the surface which influence the signal-to-noise ratio and hence the accuracy of the fit. For the simulation of spacecraft and pointing errors we also employ a similar strategy to the one used in previous works [13,36]. The orbit error is modeled assuming deviations in along-, cross-track, and radial direction on the orbital frequency due to residuals in the orbit determination. Pointing errors reflecting the uncertainties in the transmitter alignment and in the orientation of the spacecraft with respect to nadir are randomized assuming a Gaussian distribution with random azimuth. The assumed amplitudes are 5 m along-track, 2 m cross-track, and 0.5 m radial [37], and the pointing knowledge error is assumed to be 20 arcseconds (all values 1-sigma). The synthetic measurements are in the ICRF with respect to Ganymede’s center of mass and are therefore independent of the assumed rotational state or global topography. While the radial position error and the instrument measurement error directly distort the range measurement, the errors along and cross-track lead to observations being badly positioned and hence affecting the ability to separate between the topography and the rotation state.

2.4. Inversion Strategy

The implemented routines aim at estimating the mean radius and global topography of Ganymede, tidal Love number $h_2$, rotation rate, pole orientation, and the libration amplitudes for the frequencies derived in Section 2.1. The inversion process of these parameters must account for the fact that the initial locations of the laser spots are unknown due to the initially limited knowledge of Ganymede’s rotational state. This leads to a highly non-linear inversion problem which needs to be solved iteratively. Within this problem, the presence of topography is essential since on a perfect sphere no rotational state could be determined. However, for this study, the solution for the topography itself is considered to be a byproduct.

Following the synthetic laser altimetry data, we use two different approaches for the inversion. Method 1 is a least-squares inversion based on a global spherical harmonic expansion of the topography while Method 2 is an adaption of the method developed by Koch et al. [27], parameterizing the surface topography into cubic B-splines on top of a rectangular grid.

Both inversions are implemented to fit the unknown parameterized topography and rotation model $\mathbf{F}\left(\mathbf{p}\right)$ to the synthetic observations, which are the ranges ($r=r_{sc}-r_s$) as measured by the laser. Here, $r_{sc}$ is the distance of the spacecraft from Ganymede’s center of mass. The best approximation of the inversion problem and at the same time most likely solution vector $\mathbf{p}=(p_1,p_2,\cdots ,p_k)$ is determined as the minimum of the weighted misfit-function. For a linear model, this result is given by the normal equation

$$\mathbf{p}={\left({\mathbf{F}}^T{{\mathbf{C}}_{\mathbf{d}}}^{-1}\mathbf{F}\right)}^{-1}{\mathbf{F}}^T{{\mathbf{C}}_{\mathbf{d}}}^{-1}\mathbf{r},$$

(13)

with ${\mathbf{C}}_{\mathbf{d}}$ being the covariance matrix of the observation vector. However, the given problem is highly non-linear and therefore the problem is first linearized for a given a-priori parameter vector ${\mathbf{p}}^0$. Using a Taylor series expansion, the problem can be formulated as

$$\mathbf{F}({\mathbf{p}}^0+\mathsf{\Delta }\mathbf{p})=\mathbf{F}\left({\mathbf{p}}^{\mathbf{0}}\right)+\mathsf{\nabla }\mathbf{F}\left({\mathbf{p}}^{\mathbf{0}}\right)\left[\mathbf{r}-\mathbf{F}\left({\mathbf{p}}^{\mathbf{0}}\right)\right].$$

(14)

$\mathsf{\nabla }\mathbf{F}\left({\mathbf{p}}^0\right)$ is the Jacobian matrix, or design matrix $\mathbf{A}$, which is used to calculate a parameter update

$$\mathsf{\Delta }\mathbf{p}={\left({\mathbf{A}}^T{{\mathbf{C}}_{\mathbf{d}}}^{-1}\mathbf{A}\right)}^{-1}{\mathbf{A}}^T{{\mathbf{C}}_{\mathbf{d}}}^{-1}\mathsf{\Delta }\mathbf{r}$$

(15)

at each iteration. As a result, of gaps between laser profiles due to the given JUICE ground-track spacing, the resolution of the topography is limited. Moreover, the polar orbit of JUICE leads to a non-uniform distribution of laser spots on Ganymede’s surface, with a high point density at the poles and a lower point density in the equatorial region. The effect of the resulting irregular sampling is mitigated by weighting the observations with respect to the data density. This strategy has been proven useful in previous studies [26].

Due to the gaps between the GALA tracks the series of spherical harmonics in Method 1 must be cut off at an adequate degree $N_{\mathrm{est}}$ to avoid an ill-posed problem. However, since the surface topography is consequently only expressed by a finite expansion of spherical harmonics, errors will arise in the solution even in absence of measurement noise. The result is an omission error which depends on the discrepancy between the degree of the spherical harmonics expansion used for solving and the one for simulating the topography.

In non-linear least-squares inversions, regularization can improve the convergence performance of the inversion. While well-constrained a-priori information about Ganymede’s rotational state exists for small signal parameters such as the Love number $h_2$ and the libration amplitudes, no such information is available. Therefore, the strategy is as follows. Within the first iteration only the global topography (C and S coefficients) is solved, the rotational parameters are all kept constant. In the second iteration the rotation rate and spin axis orientation are estimated in addition to the refinement of the topography coefficients. For the subsequent iterations, updates are then computed for the full parameter vector $\mathbf{p}$.

However, as already noted by Koch et al. [26,27], solving for spherical harmonics comes with a significant computational load. In practice, this means that in our simulation it is not practical to solve the topography above degree and order 60 using Method 1. Method 2 has been implemented to circumvent this problem and we will investigate the impact on the results in Section 3.

The second method used in this paper is described in detail in Koch et al. (2010) [27] and has been refined by Thor et al. (2017) [38]. Instead of modeling the surface by spherical harmonics, it is parameterized on a rectangular grid. The topography in each grid cell is then described by a set of local 2D cubic B-spline basis functions. Each spline function is centered at one grid point on the map but is also non-zero in the two adjacent grid intervals in longitude and latitude in both directions. The significant advantage of the parameterization using 2D cubic B-splines is that the parameters can be kept local (instead of the global spherical harmonics expansion). The resulting normal equation matrix ${\mathbf{A}}^T{{\mathbf{C}}_{\mathbf{d}}}^{-1}\mathbf{A}$ is highly sparse, allowing for an inversion at significantly lower computational cost. Apart from the topography, Method 2 solves for $h_2$, ${\omega }_1$, and focuses on the short term libration amplitudes ${\gamma }_{7.16}$, ${\gamma }_{7.05}$, ${\gamma }_{6.26}$, ${\gamma }_{12.48}$. Long term librations have been omitted in this method, since they only contribute marginally to the understanding of Ganymede’s interior. Again, a regularization aims at keeping the topography smooth in the first iteration, while not yet solving for rotational parameters. Then, in the second iteration, the regularization is relaxed to allow for a quick and unbiased convergence.

3. Results

For the inversion based on a spherical harmonic expansion we solved for the 7-day, 50-day, and 483-day period since similar frequencies (e.g., 7.05 d and 7.16 d) cannot be distinguished by the inversion routine. We find that the differences of the true and recovered libration amplitudes are highly dependent on the degree of the expansion. The root-mean-squared errors of the libration amplitudes show a strong hyperbolic decay with increasing inversion degree. The retrieved uncertainties when simulating the topography up to degree and order 200 and solving for degree 60 (case 1) are given in

Table 1. Derived measurement uncertainties for Ganymede’s libration amplitudes by using a spherical harmonics approach.

Parameter		Error Case 1	Error Case 2
Rotation Rate	${\omega }_1$ [$\mathsf{\mu }$rad/day]	$\pm 30.5$	$\pm 0.4$
Right Ascension	${\alpha }_0$ [$\mathsf{\mu }$rad]	$\pm 547.0$	$\pm 8.1$
Right Ascension rate	${\alpha }_1$ [$\mathsf{\mu }$rad/day]	$\pm 4.2$	$\pm 0.1$
Declination	${\delta }_0$ [$\mathsf{\mu }$rad]	$\pm 168.0$	$\pm 1.7$
Declination rate	${\delta }_1$ [$\mathsf{\mu }$rad/day]	$\pm 1.4$	$\pm 0.1$
Libration 7.16 d	${\gamma }_7$ [$\mathsf{\mu }$rad]	$\pm 83.2$	$\pm 1.3$
Libration 50 d	${\gamma }_{50}$ [$\mathsf{\mu }$rad]	$\pm 114.0$	$\pm 1.2$
Libration 483 d	${\gamma }_{483}$ [$\mathsf{\mu }$rad]	$\pm 2740$	$\pm 29.5$

. The result for the rotational state is shown in comparison to the case where the topography has been simulated and solved up to degree 60 (case 2). As can be seen from the results, the artificial elimination of the omission error for the case 60-60 increases the precision of the retrieved libration amplitudes by two orders of magnitude.

Furthermore, the inversion method does not yield a statistically significant solution for the tidal Love number $h_2$. The determined uncertainties are in the same order as the Love number itself or at least much higher than estimates from cross-over methods [13]. In summary we also find that the uncertainties resulting from the here modeled distorted spacecraft positioning and misalignment are of minor importance in comparison to the omission error which arises from the finite spherical harmonics expansion of the topography.

Also, for Method 2, using cubic B-splines on a rectangular grid, we find that the result has a dependence on the resolution of the grid (Figure 4). With increasing resolution, the topography is parameterized using more coefficients and the grid spacing becomes denser than the cross-track distance. At this point, which is reached at a resolution of about 4 grid cells per degree (equivalent to spherical harmonic degree 360), unconstrained topography coefficients will soak up rotational signals and thereby bias the libration amplitudes and rotation rate. A certain range measurement can be explained by both topography (radial effect) and rotational state (lateral effect), leading to convergence problems if not treated properly. On the other hand, a fine resolution yields a lower standard deviation. Therefore, the inversion favors a relatively low resolution of about two grid cells per degree, equivalent to spherical harmonic degree 180, or a resolution of 23 km at the equator.

Table 2. Derived measurement uncertainties for Ganymede’s rotation parameters and $h_2$ by using an equi-rectangular grid approach (Method 2).

Parameter		Error Case 1	Error Case 2	Error Case 3
		($2.4\times {10}^6{\mathbf{m}}^2·{\mathit{n}}^{-2}$)	($6.4\times {10}^5{\mathbf{m}}^2·{\mathit{n}}^{-2}$)	($4.6\times {10}^6{\mathbf{m}}^2·{\mathit{n}}^{-3}$)
Love number	$h_2$	$\pm 0.0488$	$\pm 0.0256$	$\pm 0.0056$
Rotation Rate	${\omega }_1[\mathsf{\mu }\mathrm{rad}/\mathrm{day}]$	$\pm 0.11$	$\pm 0.11$	$\pm 0.04$
7.16-day libration amplitude	${\gamma }_{7.16}\left[\mathsf{\mu }\mathrm{rad}\right]$	$\pm 6.58$	$\pm 6.57$	$\pm 2.45$
7.05-day libration amplitude	${\gamma }_{7.05}\left[\mathsf{\mu }\mathrm{rad}\right]$	$\pm 4.96$	$\pm 4.98$	$\pm 2.17$
6.26-day libration amplitude	${\gamma }_{6.26}\left[\mathsf{\mu }\mathrm{rad}\right]$	$\pm 7.65$	$\pm 7.63$	$\pm 2.71$
12.48-day libration amplitude	${\gamma }_{12.48}\left[\mathsf{\mu }\mathrm{rad}\right]$	$\pm 7.15$	$\pm 7.16$	$\pm 2.76$

presents parameter uncertainties derived from 100 random realizations for this resolution. For Method 2, we investigated all three different power laws for the generation of synthetic topography. Error case 1 uses a power law ${\sigma }_n^2=2.4\times {10}^6{\mathrm{m}}^2·n^{-2}$. Error case 2 uses a power law ${\sigma }_n^2=4.6\times {10}^6{\mathrm{m}}^2·n^{-3}$. Since with smaller exponents of n, more power is contained in the unmodeled topography (above spherical harmonic degree 180), Error Case 1 leads to significantly higher uncertainties.

We consider the results provided by the second method to be more reliable because the synthetic topography could be simulated up to very high degree, which is a significant shortcoming in Method 1. Despite simulating the topography to high degree, we find similar results for the rotation rate as in Error Case 2 of the spherical harmonics solution which is considered to be overly optimistic due to the artificial elimination of the omission error. Other than in Method 1, Method 2 also provides an inversion result for the tidal Love number $h_2$. In Error Case 2 of Method 2, the value for the $h_2$ error is with about $2\%$ comparable to the result that would be obtained from cross-over measurements [13], in Error Case 1 the error would be about twice as high and in Error Case 3 $h_2$ could be determined very precisely with a value below $1\%$.

4. Implications for Ice Shell Thickness

The potentially observable libration amplitudes are in direct relation to the interior structure of Ganymede and depend in particular on the presence of a decoupling liquid layer. To assess the implications of the libration measurement by GALA we investigated three different scenarios. (1) A completely rigid body, implying no subsurface ocean, (2) a Ganymede with a decoupled rigid ice shell and (3) a Ganymede with a decoupled elastic shell.

For a rigid body without a global ocean the entire body librates. Therefore, the free libration period only depends on the $J_2$ and $C_{22}$ coefficients of the gravitational field. For small libration angles ${\gamma }_j$, the amplitudes of the libration can be approximated by the magnitude of driving force $H_j$, and the free libration frequency, ${\omega }_{\mathrm{free}}=n\sqrt{3(B-A)/C}$ (n is the mean motion) as

$${\gamma }_j=H_j\frac{{\omega }_{\mathrm{free}}^2}{{\omega }_{\mathrm{free}}^2-{\omega }_j^2}.$$

(16)

The difference in the equatorial inertia axes $B-A$ is related to the $C_{22}$ coefficient by $B-A=4M{R}^2{C}_{22}$. If ${\omega }_{\mathrm{free}}\gg {\omega }_j$, the system is far from resonance. If the libration response is equal to the forcing amplitude, i.e., ${\gamma }_j=H_j$, then ${\gamma }_j$ does not contain any information about the interior. If ${\omega }_{\mathrm{free}}\ll {\omega }_j$, the libration amplitude diminishes to ${\gamma }_j\approx 0$. Please note that very close to the resonance, when ${\omega }_{\mathrm{free}}\approx {\omega }_j$, Equation (16) is invalid, since the assumption of small angles is violated. With a $GM$ of 9887.83 km${}^3$/s${}^2$ [8], a $C_{22}$ of $38.26\times {10}^{-6}$ and a radius of 2631.2 km we can compute the free libration frequency for a coupled, rigid Ganymede to 0.034 rad/day and obtain a resulting libration amplitude of about 6 $\mathsf{\mu }$rad.

In presence of a global subsurface ocean, and assuming that the solid interior is spherically symmetric, the ice shell is mechanically decoupled from the mantle. Neglecting viscous coupling at the layer boundaries, the shell can be assumed to librate independent of the interior. In that case only the moment of inertia of the rigid ice shell $C_{\mathrm{shell}}$ gives resistance to the libration force and the free libration frequency becomes ${\omega }_{\mathrm{free}}=n\sqrt{3(B-A)/C_{\mathrm{shell}}}$ resulting in about a magnitude higher libration amplitudes.

However, in case of an elastic ice shell, Van Hoolst et al. [4] demonstrated that the gravitational torque, and thereby the free libration frequency, can be reduced by a factor $\sqrt{(k_f-k_2)/k_f}$, where $k_f=0.804$ [8] is Ganymede’s fluid Love number. Since for the non-ocean case, the response of a rigid Ganymede without liquid ocean will be similar to the librational response of an elastic Ganymede without liquid ocean. Therefore, the elastic non-ocean case is not further considered here. In case of a subsurface ocean, however, $k_2$ can be in the same order of magnitude as $k_f$, and hence reduce the libration amplitudes. For reference we use a structural model of Ganymede investigated in a previous work [13] and assume a completely elastic outer ice shell using a rigidity of $\mu =3.3$ GPa and a viscosity of $\eta ={10}^{17}$ Pas. The tidal Love number $k_2$ has been computed as a function of ice shell thickness using a matrix propagation method [39]. The libration amplitude as a function of the outer ice shell thickness is shown in Figure 5 for the four most promising forcing frequencies. Please note that gravitational and pressure coupling arising because of the equatorial flattening of the solid interior [4] is not included in our calculations but may affect the presented amplitudes. However, the shape of the inner static bulge is not known and can only be approximated, also inducing uncertainties.

Due to the low gradient at ice shell thicknesses above 50 km, the 6.26 days and 12.48 days forcing periods are not particularly informative about the ice shell thickness. However, librations on the orbital period of 7.16 days have a sufficient dependency on the ice shell thickness to provide a meaningful geodetic constraint. At an ice shell thickness of around 150 km, the gradient of the blue line in Figure 5 is 0.1 $\mathsf{\mu }$rad/km. In the most optimistic Error Case 3 of Method 2, this translates to an error of ±24 km (

Table 3. Constraint on the ice shell thickness from the 7.16-day libration period. Results are presented for the three different error cases as a function of the true ice shell thickness.

Ice Shell Thickness	Error Case 1	Error Case 2	Error Case 3
50 km	$\pm 7.4$ km	$\pm 7.4$ km	$\pm 2.8$ km
100 km	$\pm 29.1$ km	$\pm 29.1$ km	$\pm 10.6$ km
150 km	$\pm 65.3$ km	$\pm 65.2$ km	$\pm 24.3$ km

). However, in the two other error cases of Method 2, the error in ice shell thickness becomes ±65 km and it is questionable if this result would be considered valuable. Generally, for thinner ice shells it becomes easier to constrain their thickness. In case of a 50 km thick shell, the gradient increases to 0.6 $\mathsf{\mu }$rad/km allowing for an error as low as ±3 km assuming Error Case 3.

The results of GALA’s libration measurement will still have ambiguity in terms of different densities and rheologies. Therefore, these measurements will be combined with other approaches in a complete analysis involving several instruments. The result for a 150 km thick ice shell is comparable to the accuracy resulting from the $h_2$ measurement using cross-over points [13] and thereby a valuable independent constraint. In combination with the results from other experiments it would greatly contribute revealing Ganymede interiors. It is expected that the magnetometer J-MAG [40] will provide constraints on the ocean thickness and salinity, and the gravity experiment 3GM [41] provides a highly precise measurement of the gravitational parameter $GM$, the moment of inertia and independent constraints on the rotational state.

5. Conclusions

Based on the conducted analysis, we can state that a solution based on spherical harmonics is rather unpractical due to the necessity to expand the topography up to high degrees. Cutting the expansion will yield a significant omission error dominating all error sources. Since this type of error is generally avoidable it presents a significant flaw. However, it demonstrates that the topography of Ganymede can be reconstructed up to degree and order 60 and a determination of the rotation state should be feasible independently from other instruments. The results for the librations however, depend too strongly on the degree of the topography expansion to state reliable results. Also, a solution for the tidal Love number $h_2$ is not possible within this method and the error is much higher compared to a measurement from cross-over points.

The second method tested, using B-splines on an equi-rectangular grid, allows for stable results for all parameters. However, it also shows a dependency on the resolution of the grid cells. Indeed, if the grid resolution is too low, it yields similar omission errors resulting from the undersampling as in Method 1, while too high resolutions result in the topography absorbing the libration signal. An optimum gridsize has been identified at 2 cells per degree, corresponding to spherical harmonics degree 180. Using Method 2, we estimate that the libration amplitude at the orbital frequency of Ganymede is measurable by laser altimetry to about 2.4–6.6 $\mathsf{\mu }$rad and we infer that this result could constrain the thickness of an elastic ice shell, decoupled by a liquid layer from the deeper interior, with an accuracy of 24 to 63 km. In the more optimistic range of this error interval, the physical librations could be a valuable constraint on the outer ice shell in combination with the measurement of the tidal Love number $h_2$. Moreover, knowledge in the rotational parameters is prerequisite for the realization of Ganymede’s body-fixed frame and consequently for accurate mapping of its surface.

GALA is part of a geophysics package onboard JUICE including the magnetometer J-MAG and the radio science experiment 3GM. Together with the tidal Love number $k_2$ and the ocean salinity, invaluable information can be gained on the H${}_2$O layers and hence on the habitability of the largest ocean world in the Solar System.