Theory of Figures to the Seventh Order and the Interiors of Jupiter and Saturn

In the Mogrop code (Nettelmann 2017), the constant K is adjusted to fit the mass M. The mean radius R_m is adjusted to fit R_eq. The radial grid, for this application, is split into N grid points, out of which N/2 are equally distributed over 0–0.95 R_m , and the other half equally over 0.95–1 R_m . Such a choice was found to give a better match to the analytic solution than a split at 0.9 R_m or deeper. Indeed, with MOGROP, we find that the accuracy increases the farther out the separation is made, with a difference up to an order of magnitude compared to a flat distribution. The integrals in Equations (A7, A8) are converted into integrals over density by partial integration and solved by the simple trapezoidal rule. The integration of the equations of mass conservation, dm/dl = 4π l² ρ(l), and hydrostatic equilibrium, (1/ρ(l))dP/dl = dU/dl, is performed by the Runge–Kutta fourth-order method. The J_2n are computed using Equation (5) and denoted by "4,7/Ne" and shown as green curves in Figure 1.

The second code we use in our n = 1 polytrope comparison test case is an independent implementation of the ToF algorithm using the same coefficients but otherwise unrelated to the Mogrop code. The two codes are therefore expected to reproduce very similar solutions if given the same conditions. TOF-planet has previously been applied in a Bayesian study of Saturn's possible interior (Movshovitz et al. 2020). Since for that purpose it was necessary to run tens of millions of density models to draw representative statistical samples, the code had to be optimized for speed and memory usage. An optional feature allows the shape functions to be explicitly calculated on a subset of level surfaces, while the shape of the rest can be spline interpolated in the radial direction. This "skip-n-spline" trick can provide a significant speed advantage when high-resolution density profiles are needed. We find that, even when a very high resolution of the density profile is required to accurately calculate integrals over density, there is no advantage in calculating the shape functions for more than a few hundred level surfaces. The speed advantage of this optimization applies mainly to high-resolution ToF7 calculations. For lower resolution, and for most ToF4 runs, the interpolation overhead ruins the effort. (ToF7 is much slower than ToF4 for a given N owing to the many more terms appearing in each of the shape function equations.)

To validate TOF-planet with both ToF4 and ToF7 coefficients, we use it to reproduce the n = 1 polytrope test of Wisdom & Hubbard (2016). To make a direct comparison in a consistent way, some care is needed. The mass and equatorial radius are taken as in Wisdom & Hubbard (2016) and remain fixed for the duration of the calculation. (The mass is taken from the reported GM and with G = 6.6738480 × 10⁻¹¹ m³ kg⁻¹ s⁻².) However, in Wisdom & Hubbard (2016) the rotation state is given by the parameter q_rot, whereas the ToF algorithm needs the related parameter m_rot. The conversion needs the ratio R_eq/R_m, which is only available after the equilibrium shape is solved. To obtain a self-consistent solution, we fix the planet's rotation frequency ω using the value q_rot = 0.089195487. We then proceed with a guess for R_eq/R_m and therefore m_rot, solve for the shape function and gravity field, integrate the hydrostatic equilibrium equation to solve for pressure everywhere, update the density everywhere to match the polytropic relation, renormalize the level radii grid to match the reference equatorial radius, renormalize the density to match the reference mass, recalculate m_rot for the updated R_eq/R_mean ratio, and rerun all the steps until a self-consistent solution is found.

In this test, for both ToF7 and ToF4, we compute all integrals with the trapezoidal rule and constant grid spacing. With TOF-planet, we experimented and found that different integration schemes and grid spacing schemes did not reduce substantially the number of grid points required for a given precision. This should not discourage, however, future users of our ToF7 tables from optimizing their grids and integration schemes for their particular cases.

The resultant J_2n values appear in blue and are denoted by "4,7/Mo" in Figure 1.

As in previous work (Ni 2020), we apply the CEPAM code (Guillot & Morel 1995) to calculate the gravity field and shape using ToF5 (Zharkov & Trubitsyn 1975). Here we have expanded the code to address the case of the rotating n = 1 polytrope.

For the n = 1 polytropic EOS P = K ρ², the constant K is determined in terms of mass conservation, and the mean radius R_m is adjusted to reproduce the equatorial radius R_eq = 71,492 km. The initial density distribution is first given by that of a nonrotating n = 1 polytrope $\rho (z)={\rho }_{c}\sin \pi z/\pi z$. The figure function s_2k (z) and total potential U(z) are computed using the ToF5 as described in Zharkov & Trubitsyn (1975) and Ni (2019). In its original version, CEPAM uses an automatic grid refinement method that distributes the grid points in a way that a distribution function of the variables pressure, temperature, mass, radius, and luminosity changes by a constant amount at the grid points. This method requires smooth behavior of the variables and their derivatives. B-splines are used as the interpolating polynomials, which exhibit the desired properties. However, for number of grid points larger than 10³, we did not obtain stable solutions with CEPAM. Therefore, for a higher number of grid points we switched to our own solution of the pressure profile using the trapozoidal rule

$$\begin{eqnarray}&&P({z}_{j})=P({z}_{j-1})+0.5[\rho ({z}_{j})+\rho ({z}_{j-1})][U({z}_{j})-U({z}_{j-1})]\end{eqnarray} \tag{ 6 }$$

with the outer boundary condition P(R_m ) = 0 Mbar. In this case, and in view of the fact that gravitational harmonics show greater sensitivities to the external levels of a planet, more radial grid points are taken for the outer part: N/2 equally distributed over 0.85–1 R_m and the other half equally over 0–0.85 R_m . With CEPAM, we find that this choice yields a modest optimum in accuracy.

Finally, the new density distribution is obtained from the n = 1 polytropic EOS $\rho ({z}_{j})=\sqrt{P({z}_{j})/K}$. This procedure is iteratively performed until all changes in the density distribution are reduced to within a specified tolerance.

In the work of Ni (2020), the gravitational zonal harmonics are calculated as weighted integrals over the internal density distribution ρ(z) using the resulting figure functions s_2k (z),

$$\begin{eqnarray}&&{J}_{2i}=-\displaystyle \frac{2\pi }{{{MR}}_{{eq}}^{2i}}{\int }_{-1}^{+1}d\cos \theta {P}_{2i}(\cos \theta )T({R}_{m},\theta ),\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&T({R}_{m},\theta )={\int }_{0}^{{R}_{m}}{dl}\ \rho (l,\theta ){r}^{2i+2}\left(\displaystyle \frac{{dr}}{{dl}}\right).\end{eqnarray} \tag{ 8 }$$

Using the scaled mean radius z = l/R_m and abbreviating the level surface Equation (4) as r = zR_m [1 + Σ(z, θ)], one can express the function T(R_m , θ) as

$$\begin{eqnarray}\begin{array}{rcl}T({R}_{m},\theta ) & = & \displaystyle \frac{3M\,{R}_{m}^{2i}}{4\pi }\left\{\displaystyle \frac{{[1+{\rm{\Sigma }}(1,\theta )]}^{2i+3}}{2i+3}\ \displaystyle \frac{\rho (1,\theta )}{\bar{\rho }}\right.\\ & & -\left.\,{\displaystyle \int }_{0}^{1}{dz}\ \displaystyle \frac{{[1+{\rm{\Sigma }}(z,\theta )]}^{2i+3}}{2i+3}\,{z}^{2i+3}\ \displaystyle \frac{d\rho (z,\theta )/\bar{\rho }}{{dz}}\right\}.\end{array}\end{eqnarray} \tag{ 9 }$$

Table 1 shows a comparison of the even harmonics obtained from Equation (7) and from Equation (5) for a typical number of grid points N = 2000. The numerical accuracy in J₂ is almost the same for both of them. But the results from Equation (7) are in better agreement with the analytic Bessel-function-based solution for J₄–J₁₀, where the numerical accuracy for Equation (7) is about 1–2 orders of magnitude better than that for Equation (5).

Table 1. Comparison of the J_2n Obtained from Equation (7) and from Equation (5)

Method	J2 × 106	J4 × 106	J6 × 106	J8 × 106	J10 × 106
Bessel	13 988.51	−531.8281	30.118 32	−2.13212	0.174 07
Equation (7)	13 988.54	−531.8292	30.119 89	−2.13048	0.17446
$| {\rm{\Delta }}{J}_{2i}/{J}_{2i}^{\mathrm{Bessel}}|$	2.42 × 10−6	1.99 × 10−6	5.22 × 10−5	7.67 × 10−4	2.23 × 10−3
Equation (5)	13 988.55	−531.8207	30.135 06	−2.10486	0.19555
$| {\rm{\Delta }}{J}_{2i}/{J}_{2i}^{\mathrm{Bessel}}|$	2.52 × 10−6	1.40 × 10−5	5.56 × 10−4	1.28 × 10−2	1.23 × 10−1

Note. The numerical values in this table are for a typical number of grid points N = 2000 and using the ToF5 coefficients of Zharkov & Trubitsyn (1975). The Bessel solution is taken from Wisdom & Hubbard (2016)

Download table as:  ASCIITypeset image

In Figure 1 we show the relative deviations of the calculated even harmonics from the analytic solutions of Wisdom & Hubbard (2016) as a function of the number N of radial grid points.

With the CEPAM code and ToF5, denoted by 5/Ni in the figure, the numerical accuracy of all the calculated J_2i shows good convergence with an increased number of grid points. When the number of grid points is increased beyond ∼ 10³, the numerical accuracy in J₄–J₁₀ falls below the current observational uncertainty (Juno D20) reported in Durante et al. (2020). Moreover, the accuracy in all the harmonics J₂–J₁₀ is better than the CEPAM-WH16 results from Guillot et al. (2018), who reportedly applied ToF4, by a factor of roughly 5–100.

Using ToF4 and ToF7 in conjunction with the Mogrop code, denoted by 4/Ne and 7/Ne in the figure, the accuracy significantly improves with denser grid points. Apparently, this code requires a factor of 100 more radial grid points than CEPAM to obtain the same accuracy in J₂ and J₄. For these low-order harmonics, ToF7 versus ToF5 provides a negligible improvement in accuracy. The situation changes with J₆. Here, CEPAM with ToF5 levels off at a relative uncertainty of ∼ 5 × 10⁻⁵, while the higher accuracy of ToF7 versus ToF5 becomes evident as N increases beyond 20,000. However, the typical number of grid points used for planet interior models ranges between 2000 and 4000. For such N values, the numerical accuracy in J₆ with 7/Ne is about the same as the current observational uncertainty reported in Durante et al. (2020). ToF7 begins to pay off with J₈ and higher, even with Mogrop, where J₈ becomes an order of magnitude better than the observational one, 2.5 orders of mag in J₁₀, and 3 orders of mag in J₁₂. We conclude that ToF7 is sufficiently accurate to address the influence of the winds on J₆ and higher, given current observational uncertainties.

Using the independent TOF-planet code of Movshovitz et al. (2020), we obtained similar J_2n values to those with the Mogrop code; compare the blue and green lines in Figure 1. In particular, the results for J₄–J₈ with ToF4 after convergence with grid point number N agree, indicating that the remaining errors ΔJ_n /J_n of 2.5 × 10⁻⁴ in J₄, 10⁻² in J₆, and 10⁻¹ in J₈ are due to the truncation of the expansion of ToF4. The ToF7 values also agree when convergence is reached, though this applies only to J₁₀ and J₁₂ for high values of N > 10,000. Before convergence with N is reached, the ToF7 errors deviate by a factor of a few, suggesting that the radial grid discretization error matters, which can differ between different implementations even if they use the same trapezoidal rule.

The similar accuracy of Mogrop and TOF-planet is due to similar methods for the numerical integration over density (trapezoidal rule) and of the differential equations dm/dr and dP/dr (Runge–Kutta). There is room for improvement. As an example, we extrapolate the J_2n values using linear regression on the three solutions for N = 1000, 2000, and 4000 for each J_2n . In Figure 1, the resulting accuracy is conservatively compared to the result for N = 9000, corresponding to the computational cost that scales linearly with N and a small offset for each run. Apparently, the gain in accuracy amounts to two orders of magnitude in J₂ and J₄, and it is still better than compared to using 10⁵ grid points. This suggests that methods other than simply skyrocketing the number of grid points may help to improve the accuracy of numerical J_2n computations.

We note that at present ToF is entirely outperformed by the CMS method with regard to accuracy. Militzer et al. (2019) reported relative inaccuracies of only 7.3 × 10⁻⁹ in J₂, 2.1 ×10⁻¹⁰ in J₄, 3.6 × 10⁻⁸ in J₆, 4.2 × 10⁻⁸ in J₈, 1.1 × 10⁻⁷ in J₁₀, and 6.7 × 10⁻⁹ in J₁₂ for N = 2¹⁷ = 131,072 CMS layers, out of which only 512 are treated explicitly, while the shape of intermediate ones is obtained by interpolation.

The models of this work assume a four-layer structure. By Y_i we denote the helium mass fraction in layer number i with respect to the H/He system. Layer 4 is a compact rocky core. Layer 1 is an atmosphere with a helium mass fraction of Y₁ = 0.238 as measured by the Galileo entry probe. Layers 2 and 3 have the same helium abundance (Y₂ = Y₃), which is adjusted to yield an overall helium mass fraction Y = 0.2700(4). A possible He-rain region in Jupiter is represented by a jump in helium abundance between layers 1 and 2. Transition pressures of 1–4 Mbar for P₁₂ are considered, which are typical pressures where the Jupiter adiabat reaches the closest point to the H/He demixing boundary of H/He phase diagrams for protosolar H/He ratios as predicted by first-principles simulations (Lorenzen et al. 2011; Morales et al.2013; Hubbard & Militzer 2016; Schoettler & Redmer 2018). Adjusting the local He abundance to the local P–T conditions along the phase boundary yields an approximately linear increase in Y; however, the gradient and width of the He-rain region depend on the temperature profile assumed therein (Nettelmann et al. 2015), which may range in Jupiter from adiabatic to modest superadiabaticity (Mankovich & Fortney 2021). A recent analysis of reflectivity data obtained for H/He samples that were pre-compressed to 2–4 GPa in diamond anvil cells and further shock-compressed to 60–180 GPa using the OMEGA layer indicates that an even larger portion of the Jupiter adiabat may intersect with the H/He phase boundary, as at the highest pressure where evidence of demixing is seen, 150 GPa, the measured temperatures were 10,000 K (Brygoo et al. 2021). Assuming a flat T(P) phase curve at Mbar pressures, such a temperature corresponds to ∼8 Mbar along the Jupiter adiabat (Hubbard & Militzer 2016).

Although He droplets may carry specific elements such as Ne with them downward (Wilson & Militzer 2010) and affect the metallicity between the He-depleted outer and He-enriched inner regions, we assume constant heavy-element mass fractions across that boundary (Z₁ = Z₂). Finally, between layers 2 and 3, the heavy-element mass fraction is allowed to change. We use either a constant Z₃ value, implying a jump in Z at a transition pressure P₂₃, or a Gaussian-Z₃ profile that starts with Z(P₂₃) = Z₂ and smoothly increases toward a maximum ${Z}_{3,\max }$ at P = 38 Mbar near the core. The choice of 38 Mbar is arbitrary and was taken to be just above the usual core–mantle boundary pressures, which are found to be around 40 Mbar in Jupiter. The two free parameters in that setup to adjust J₂ and J₄ are Z₁ and Z₃ or ${Z}_{3,\max }$.

We employ the CMS-2019 EOSs for H and He (Chabrier et al. 2019) and mix them with the water EOS H2O-REOS with respect to density only. The T–P profile is that of the H/He adiabat, which begins at T = 166.1 K at 1 bar. We construct curves of constant entropy (adiabats) by using the specific entropy values s_H(P, T) and s_He(P, T) provided in the tables for hydrogen and helium (Chabrier et al. 2019) after adding a composition-dependent entropy of mixing term ${s}_{\mathrm{mix}}({X}_{{{\rm{H}}}_{2}},{X}_{\mathrm{He}})$. For the concentrations, we assume that helium is nonionized and that hydrogen is either molecular or ionized, taking the degree of ionization as in Nettelmann et al. (2008). Since we found these H/He adiabats to be too dense to yield Jupiter models with nonnegative atmospheric metallicity, we also perturb that adiabat toward lower densities as described in Section 4.3.

In Figure 2 we show the even J₂–J₁₀ values from rigidly rotating Jupiter models. Models adjusted to the Juno observations of J₂ and J₄ (Durante et al. 2020) are shown in bluish color, while models adjusted to the wind-corrected J₂ and J₄ values by the corrections of Kaspi et al. (2018) applied to the J₂, J₄ values of Durante et al. (2020) are shown in reddish color.

Zoom In Zoom Out Reset image size

Due to imperfect fit to the J₂, J₄ values, the scatter in model J₂ and J₄ values is larger than the observational uncertainty. For J₄, the scatter ΔJ₄/J₄ is about ±2 × 10⁻⁴ and of the same size as the relative uncertainty due to using ToF7 in the Mogrop code, while for J₂ the latter relative uncertainty is with 7 × 10⁻⁵ overwhelming. Still, these relative deviations are too small to matter for the inferred metallicities. Guillot et al. (2018) allowed for a similarly wide scatter in J₂ model values of ±3.4 × 10⁻⁵ and a wider scatter in J₄ of ±10⁻³ relative deviations. They found that, nevertheless, the high-order moments J₈ versus J₆ and J₁₀ versus J₈ were strictly confined to a straight line. We confirm that behavior.

Notably, the model ∣J₈∣ values are higher than the observed value, and a trend in that direction is also seen for J₁₀, although the model J₁₀ values are still within the 3σ observational uncertainty. Wind models assuming rotation along cylinders indeed predict a slight decrease in ∣J₈∣ and J₁₀ if the observed wind profile of the southern hemisphere is applied to the entire surface, while they predict an enhancement if the wind profile of the northern hemisphere is used (Hubbard 1999). The wind model by Kaspi et al. (2018) that was adjusted to explain the odd moments of Jupiter observed by Juno yields a correction qualitatively in the direction as predicted for the southern winds and seen in the model values for rigid rotation, albeit quantitatively stronger by a factor of two. This deviation may have many reasons; clearly, further exploration of the connection between interior and wind models is desirable.

For J₆, the uncertainties from the application of ToF7 in the Mogrop code, from observations, and from the wind contribution are all small and of the same size. In contrast, model assumptions such as the location of layer boundaries not only have a larger influence on J₆ but also yield a scatter around the observed value. Hence, we conclude that Jupiter's J₆ is unique in that it neither is adjusted nor seems to be significantly influenced by the winds, and therefore it offers an additional parameter to further constrain interior models. We find that models with a Gaussian Z₃ and an abrupt He-poor/He-rich transition at P₁₂ = 1–2 Mbar yield J₆ = ( 34.11–34.18) x10⁻⁶, slightly below the observed value, while models with that transition deeper inside at P₁₂ = 2.5–3 Mbar yield J₆ = ( 34.23–34.28) x 10⁻⁶, slightly above the observed value. Constant-Z₃ profiles and P₁₂ = 2 Mbar stretch from J₆ = ( 34.18 to 34.28) times 10⁻⁶ around the observed value (34.200 7 ± 0.0067) x 10⁻⁶ upon shifting P₂₃ from 5 Mbar deeper down to 18 Mbar. Taken at face value, Jupiter's observed J₆ value indicates that the He-depleted/He-enriched transition occurs at around 2–2.5 Mbar.

In Figure 3 we show the radial heavy-element distribution of some of the models with Gaussian Z in layer No. 3. Models with unmodified H/He adiabat appear in light blue and are described in Section 4.3.1, while models with modified H/He adiabat appear in red and are described in Section 4.3.2.

Zoom In Zoom Out Reset image size

4.3.1. Unmodified H/He Adiabat

4.3.2. H/He Adiabat Modification

Models with negative Z values are, needless to say, not considered a viable solution. There are two obvious ways how negative Z values in the atmosphere and outer envelope can be circumvented. One possibility is to assume a superadiabatic region above the region where J₄ is most sensitive, which—for a polytropic Jupiter model—is in the molecular envelope at around 50 GPa (0.9 R_J ). A superadiabatic temperature profile may result from stable stratification. Christensen et al. (2020) showed that meridional flows in a stably stratified, slightly conducting region slow down the strong zonal flows and suggested the existence of such a region in Jupiter as an explanation for the truncation of the zonal flows, which, according to recent combined analysis of magnetic field and gravity field data, occurs rather sharply at 0.97 R_Jup (Galanti & Kaspi 2021).

However, stable stratification does not necessarily result in a superadiabatic temperature profile. Depending on its origin, stable stratification can also be accompanied by a subadiabatic temperature profile. For instance, in the absence of alkali metals the opacity of the H/He fluid becomes sufficiently low for the intrinsic heat to be transported by radiation along a subadiabatic radiative gradient (Guillot et al. 1994, 2004), leading to subadiabatic stable stratification according to the Schwarzschild criterion. Superadiabatic gradients are predicted in a Ledoux-stable, inhomogeneous medium of upward-decreasing mean molecular weight. Clouds formed by condensibles of higher molecular weight than the background composition has can induce Ledoux stability if their abundance is high enough. Such a scenario has been proposed for the presumably water-rich atmospheres of the ice giants (Leconte et al. 2017). Water clouds may occur in Jupiter at 100 bars, and silicate clouds at 1000 bars, both below the level that so far could be probed by the Galileo entry probe (22 bars) and Juno MWR remote sensing (Li et al. 2020). Therefore, clouds are candidate causes for superadiabatic stable stratification.

Another possibility to avoid negative metallicities is to perturb the H/He adiabat toward lower densities. The CMS-19 hydrogen EOS shows excellent agreement with a variety of experimental data ranging from shock compression experiments for H and D at various initial conditions to isentropic compression (Chabrier et al. 2019). At 50 GPa, the H EOS is even slightly stiffer than the experimental data. Only the helium EOS shows significantly higher densities in the 20–150 GPa area than inferred from the shock compression experiments (Chabrier et al. 2019). Although the good agreement between the theoretical P–ρ relations and the experiments, as well as between CMS-19 H/He adiabats and MH13 EOS adiabats, is far from suggesting that the CMS-19 H/He EOS would significantly overestimate the density along the Jupiter adiabat, we here perturb it toward lower densities. We conduct a three-parameter study where we vary the maximum deviation $\delta {\rho }_{\max }$, the pressure entry point P_start, and the pressure exit point P_end. Between P_start and P_end, δ ρ adopts its maximum at the logarithmic mean pressure value and is otherwise linearly interpolated as $\delta \rho (\mathrm{log}P)$. We explore P_start values between 1 and 50 GPa and P_end values between 50 and 150 GPa. The smaller P_start and P_end, the lower the ∣J₄∣ and ∣J₄∣/J₂ ratio. We find P_start ≤ 30 GPa necessary in order to have a noticeable influence on J₄. Conversely, higher P_end values lead to a stronger reduction on J₂. The question we ask is, for what values of $\delta \rho$, P_start, and P_end can a model be found with a 1 × solar homogeneous Z?

For a homogeneous, unperturbed adiabat at Y₁ = 0.238 and Y₂ adjusted to meet Y = 0.27, both ∣J₄∣ and J₂ turn out to be significantly too large. This is in contrast to the result by Debras & Chabrier (2019), who could match J₂ at 1 × Z_Gal. Thus, we need P_start to be sufficiently low for J₄ and P_end sufficiently high for J₂. We find such an optimized solution for P_start = 26 GPa, P_end = 150 GPa, and δ ρ = −0.1257, i.e., a maximum reduction of the H/He adiabat by 12.57%. The resulting P–ρ profile is shown in Figure 4(c), while the ensemble of models in the J₂–J₄ and J₆–J₄ space is shown in Figures 4(a) and (b). Notably, among the wide spread of intermediate models in the J₂–J₄–J₆ space, the one model (red cross) that meets the Juno J₂, J₄ values yields J₆ = 34.20 x 10⁻⁶, in excellent agreement with the Juno observation. For this model, the total Z amounts to 15.6 M_E, in good agreement with the formation model of Lozovsky et al. (2017). This exploration suggests that winds on Jupiter have a negligible influence on J₆.

Zoom In Zoom Out Reset image size

4.3.3. Models for Enhanced 1 bar Temperature

In this section, we present Jupiter models for T_1bar = 175 K and T_1bar = 180 K. Here, we do not modify the adiabat or EOS, and we adjust the J₂, J₄ model values to the wind-corrected observed values using the corrections of Kaspi et al. (2018).

Such warmer models are not preferred, first, because these 1 bar temperatures significantly exceed the Galileo entry probe measurement of 166.1 K. This would not have posed a problem if a mechanism had been studied that predicted a superadiabatic region underneath the 22 bar region, wherein temperatures and thus entropy would rise to the level corresponding to these or even higher 1 bar temperatures. Clouds may have a warming effect if they stabilize the region of condensation (Leconte et al.2017); however, latent heat release from condensation opposes this effect and leads to a cooler interior underneath the cloud region, as has been discussed for ice giant atmospheres (Kurosaki & Ikoma 2017). Second, Jovian adiabats for different H/He EOSs tend to intersect with H/He demixing curves at best in a small region at 1–3 Mbar. At present, only the rather cool MH13 EOS-based Jupiter adiabat for 166.1 K shows a clear intersection by about 450 K (Hubbard & Militzer 2016) with the state-of-the-art first-principles-based H/He demixing curve of Morales et al. (2013), while the intersection with the lower demixing curve T_dmx(P) of Schoettler & Redmer (2018) is only marginal (Mankovich & Fortney 2021). Since enhancing T_{1 bar} from the Galileo value of 166.1 K by only 14 K leads to an enhancement by ∼350 K at 1 Mbar and even by 460 K at 2 Mbar according to our CMS-19 EOS-based Jupiter adiabats, higher surface temperatures might let the demixing region in Jupiter entirely disappear. We stress that although our CMS-19 Jupiter adiabat for T_{1 bar} = 166.1 is rather dense, it is also rather warm and, with 5700 K at 1 Mbar and 6840 K at 2 Mbar, outside of the first-principles-based demixing regions (Morales et al. 2013; Schoettler & Redmer 2018).

On the other hand, the recent experimentally predicted phase boundary inferred from an observed upward jump in reflectivity at 0.93 Mbar and downward jump at 1.5 Mbar, which are interpreted as the entry and exit of the compressed H/He sample in and out from the demixing region (Brygoo et al. 2021), suggests high demixing temperatures of 10,000 K. Primarily it is this finding that motivates us to allow for higher surface temperatures and for higher transition pressures, which we allow to reach the maximum where the core mass disappears, or for practical reasons drops below 1 M_E.

In Figure 5 we show the resulting outer envelope metallicity Z_atm and J₆ values. Obtaining 1× solar metallicity requires T_1bar of 180 K (orange curve) or higher. While not negligible, additional uncertainties in the atmospheric helium abundance and J₄ are small and not considered here. Notably, for T_{1 bar} = 180 K and at the transition pressure between the He-poor and He-rich regions at P₁₂ = 6 Mbar, we obtain 1× solar metallicity throughout the interior down to ∼ 0.4R_J, thus a largely solar-metallicity envelope. At 6 Mbar, the temperature amounts to 10,400 K and is thus at the upper limit of the experimentally inferred demixing temperature (Brygoo et al. 2021). For that model, the static J₆ value is consistent with the observed value and its small wind correction according to Kaspi et al. (2018). As is well known (Nettelmann et al. 2012), Z₁ rises with P₁₂.

Zoom In Zoom Out Reset image size

If there were no uncertainties in the H/He EOS, these models would suggest that the internal Jupiter adiabat lies at higher entropy than the observed adiabat down to 22 bars, and that Jupiter's envelope metallicity is not much higher than 1 × solar.

The Saturn models of this work are built in the same manner as the Jupiter models described in Section 4.1, although in the real planets, helium rain may induce a dichotomy (Mankovich & Fortney 2021). We fit the Saturn models to the observed J₂, J₄ values without accounting for the wind corrections. We assume that a rotation rate of 10:32:45 hr, as suggested by Helled et al. (2015), would yield a best match of interior models to the observed pre-Cassini Grand Finale gravity and Pioneer and Voyager shape data. Within the given uncertainty of 46 s, this value is consistent with the more recently suggested rotation rates of 10:33:34 hr (Militzer et al. 2019), using the Cassini Grand Finale gravity and the same shape data, and with the rotation rate of 10:33:38 hr${}_{-1m19s}^{+1m52s}$ inferred from the comparison of Saturn ring wave frequencies observed by Cassini with theoretical predictions for f-mode frequencies as a function of the planet's rigid-body rotation rate (Mankovich et al. 2019). We set the 1 bar surface temperature to 135 K in accordance with the Voyager measurement of 135 ± 5 K (Lindal 1992). The outer boundary is placed at a reference radius for the J_n of 60,330 km, which corresponds to the 0.1 bar level.

Saturn's atmospheric He abundance can be considered poorly known, as different estimates only agree in finding depletion compared to the protosolar value Y_proto ∼ 0.27 but disagree about the level of depletion. The lowest estimate Y₁ = 0.06 ± 0.05 stems from a combined Voyager radio occultation and infrared spectra analysis (Conrath et al. 1984), while the highest estimate Y₁ = 0.18–0.25 stems from a reanalysis of only the Voyager infrared data (Conrath & Gautier 2000). More recent Cassini-data-based estimates fall in between, ranging from Y₁ = 0.075–0.13 from Cassini infrared remote sensing (Achterberg & Flasar 2020) to Y₁ = 0.158–0.217 from Cassini stellar occultations and infrared spectra (Koskinen & Guerlet 2018). A low value of 0.07 ± 0.01 is independently inferred from shifting the most recent H/He phase diagram of Schoettler & Redmer (2018) to reproduce the He abundance measurement by the Galileo entry probe on Jupiter, in conjunction with the MH13 H/He EOS and adiabat (Mankovich & Fortney 2021). When the same procedure is applied to the H/He phase diagram of Lorenzen et al. (2011) in conjunction with the SCvH-H/He EOS, which both are now outdated, the yield is Y₁ = 0.13–0.16 (Nettelmann et al. 2015), in between the most recent observational estimates (Achterberg & Flasar 2020; Koskinen & Guerlet 2018). Here we construct models for the two moderate depletion values Y₁ = 0.14 and Y₁ = 0.18 using the CMS-19 EOS and allow for a wider spread of 0.06–0.18 when using the modified CMS-19 adiabat. For comparison, Galanti et al. (2019) used Y₁ = 0.18 ± 0.07, Militzer et al. (2019) used Y₁ = 0.18–0.26, and Ni (2020) used Y₁ = 0.12–0.23. Lower Y₁ values yield higher atmospheric and higher maximum deep interior metallicities (Militzer et al. 2019; Ni 2020).

Saturn thermal evolution models with H/He phase separation and helium rain predict that in Saturn helium rains down to the core, forming an He-rich shell of 0.90–0.95 mass percent helium atop the core (Püstow et al. 2016; Mankovich & Fortney 2021). If the helium abundance between the onset pressure of H/He phase separation and the He shell follows an H/He phase diagram, it increases gradually with depth. However, these thermal evolution models assume for simplicity a constant low envelope metallicity. How deep He droplets can sink in a high-Z and thus higher-density deep interior, as predicted by Saturn models constrained by gravity and ring seismology data (Mankovich & Fuller 2021), remains to be investigated. For simplicity, we here represent the He gradient by a jump in the He abundance at a pressure of P₁₂ = 2–4 Mbar deep within the He-rain region, for which predicted onset pressures in present Saturn range from 0.8 Mbar (Militzer et al. 2019) to 2 Mbar (Mankovich & Fortney 2021). Below P₁₂, we keep the He/H ratio constant with depth.

Typical core-envelope pressures are around 15 Mbar. We vary P₂₃ between 4 and 7 Mbar and let the Gaussian-Z₃ profiles adopt its maximum at 12 Mbar if the core is compact (Z = 1) or farther out at 6–12 Mbar if the core is dilute and thus more extended. Dilute cores are created by setting Z = 0.6–0.7 in the core, the remaining constituent being the H/He/Z mix from envelope layer No. 3 above.

Figure 6 shows the even J_n values from our Saturn models, from the two Saturn models of uniform rotation (UR) in Iess et al. (2019), and from the Cassini Grand Finale observed values (Iess et al. 2019). Unlike the case of Jupiter, Saturn's observed even ∣J_n ∣ values are clearly enhanced over the model predictions for n ≥ 6. The enhancement can be explained by rotation along cylinders that rotate approximately but not exactly with the observed speeds of clouds in the equatorial region and midlatitude region up to ±40° (Militzer et al. 2019). The observed J_n can also be explained by a thermal wind if a little deviation of the deeper wind speeds from the cloud speeds is allowed (Galanti et al. 2019). The enhancement of the even J_n by the zonal winds is quite substantial. Already for J₆, we find a 4.2%–5.3% influence, although it diminishes to 2.4%–3.7% for the modified adiabat. Iess et al. (2019) obtain slightly lower UR model J₆ values and a 5.5% effect, while Galanti et al. (2019), who allow for a much wider scatter in model J₄ values of ±40 × 10⁻⁶, obtain UR model J₆ values up to 87 × 10⁻⁶, which encompasses the observed value. However, the mean of their distribution for fast, uniform rotation lies at 82 × 10⁻⁶, implying a 5.5% influence of the winds on J₆, consistent with this work. The strong influence of the winds seen in J₆ suggests that also J₄ and J₂ are affected by the winds. In Section 6, we investigate whether the observed winds are consistent with a smaller influence on J₆ than found in previous work (Iess et al. 2019).

Zoom In Zoom Out Reset image size

Saturn interior models allow for higher atmospheric metallicities than Jupiter interior models when using the same H/He EOS. For instance, Nettelmann et al. (2013) obtains 1.5–6 × solar for Saturn and fast rotation of 10:32:00 hr, but only 0–2.5 × solar for Jupiter when applying H-REOS.2 EOS. Wahl et al. (2017b) obtained only 0–0.7 × solar for Jupiter using MH13 EOS, while Militzer et al. (2019) obtained 1–4 × solar for Saturn, consistent with Ni (2020), who obtained 0–6× for Saturn by considering a wide range of atmospheric He abundances, rotation rates, and wind corrections.

Here we obtain Z₁ = 0.02–0.06 (1.5–4 × solar) for nominal He abundances Y₁ of 0.14–0.18, and we obtain compact nonexclusive core masses of 5–8.6 M_E, meaning that solutions with lower core masses are expected to be possible for deeper transition pressures than considered here. Representative Z and ρ profiles are shown in Figure 7. Compact rocky cores yield higher central densities than suggested by the 16th–84th percentile probability range of models of Movshovitz et al. (2020), which are constrained by the gravity data. The latter models agree well with density distributions with inhomogeneous Z profiles constrained by Cassini ring seismology data (Mankovich & Fuller 2021). However, we were not able to find a model with a low-density core and the original H/He EOS, as such cores extend far out and yield J₂ values that are too large. Applying the same modification to the Saturn adiabat as for the Jupiter-optimized adiabat (see Section 4.3.2), we are able to obtain Saturn models with extended, low-density cores, in agreement with the likelihood distributions. As we mix H/He, with little addition of ice, into the rocky core region, we need rather high H/He amounts of 30%–40% by mass (Figure 7, right panel) to reduce the core densities to 6–7 cm⁻³ (left panel). This is consistent with Mankovich & Fuller (2021), who need 30%–40% of H/He for rocky cores but only 0%–10% for icy cores. When also allowing the atmospheric He abundance to decrease down to 0.06, we obtain up to 7 × solar atmospheric Z.

Zoom In Zoom Out Reset image size

Even in our dilute core models, the high-Z part is concentrated in the innermost 0.4 R_Sat at pressures above 4 Mbar. For comparison, seismic constraints are required (Mankovich & Fuller 2021) to extend the Z-gradient zone, depending on the assumed functional form of Z(r), out to 0.6–0.7 R_Sat, where the pressure is around 1 Mbar. While we did not explore such extended Z gradients in order to keep them separated from the He-gradient zone, which we placed at at 2–4 Mbar, this comparison suggests that in the real Saturn heavy-element and helium gradients overlap. Together, the Z and He abundance profiles of our models suggest that the compositional gradient required to explain the ring seismology data could be due to both a diffuse core (inner region) and helium rain (outer region).

In comparison to Jupiter, the total amount of heavy elements is clearly higher in Saturn. Models with the unperturbed H/He adiabat yield M_Z ∼ 12.6–13.6M_E for Saturn and 7.5–10.1M_E for Jupiter (Section 4.3.1). Lower densities along the H/He adiabat and warmer interior temperatures would increase these values.

We did not include Z in the computation of the adiabatic P−T profile. This leads to a slight overestimation of the temperatures along the adiabat. Mixing first the EOSs H/He-REOS and H2O-REOS linearly and then computing the adiabats as a function of Z using thermodynamic integration described in Nettelmann et al. (2012) shows that 10× solar water would lower the temperatures by only −100 K in the 10–100 GPa region relevant for J₂ and J₄. Conversely, an adiabat more rich in atomic helium would be warmer. Considering molecular volatiles in the entropy calculation would tend to make the adiabat slightly cooler and denser and therefore lead to even lower envelope metallicities, but our estimate shows that this effect should be small.

Inspired by the recent experimentally derived H/He demixing boundary that extends over a large region from ∼0.9 Mbar to ∼10,000 K (8–10 Mbar) in Jupiter (Brygoo et al.2021), one could consider helium abundances that increase over a wide region, allowing for more heavy elements to replace helium. However, our variation of the H/He adiabat showed that reduced densities are needed near the top and beyond (∼20 GPa) the demixing region. Exploration of the helium abundance profile deep inside may thus have too little influence to solve the low atmospheric metallicity problem in Jupiter.

Guillot et al. (2018) constrained the maximum amplitude of a deep wind that would extend along cylinders all the way to the center and be consistent with the even J_n to < 10 m s^–1. Kong et al. (2018) found that a flow with 1 m s^–1 down to 0.8 R_J can explain the odd J_n , but including the influence of the induced magnetic field through ohmic heat dissipation bounded by the total convective power (Wicht et al. 2019), Li et al. (2020) are able to limit this depth to only 0.96R_Jup results. Moore et al. (2019) even constrain the flow velocity to a few millimeters per second at depths of 0.93–0.95 R_J by explaining the observed secular variation of the magnetic field with advection by the flow.

In contrast, in order to lift Jupiter's atmospheric metallicities substantially, a much stronger and retrograde deep wind in the interior where J₂ and J₄ are sensitive would be needed. This deep wind must not be seen in the high-order gravity data, in the secular variation of the magnetic field, or in the System III rotation period derived from magnetic field observations. It would be seen in the moment of inertia, the static Love numbers, and the shape. At present, there is no indication of a strong (>10 m s⁻¹) wind in Jupiter's deep interior.

With both the CMS and ToF methods the interior models are derived from a self-consistent, static solution between the gravity field and the shape; however, the shape and the gravity field of the planet can be influenced by various dynamic effects.

For instance, Kong & Zhang (2020) propose that the winds are shallow while convective motions could induce a zonal flow disjunct from the surface winds. They find a dynamic influence of 1 × 10⁶ in J₂ and 0.2 × 10⁶ for J₄. While small, this effect on J₄ could be noticeable in the interior models. However, this estimate of the dynamic contribution due to convective motions is based on an Ekman number of 5 × 10⁻⁵, about 10 orders of magnitude larger than in the real Jupiter and Saturn. It is therefore possible that the dynamic contribution from convective motions on the low-order J_2n is smaller in the real planets.

For Saturn, the uncertainty in its deep rotation rate maps on an uncertainty in equipotential shape of about 120 km (Helled & Guillot 2013), far outweighting other influences like from the winds, which lift the dynamical height above a reference isobar to no more than ∼20 km (Buccino et al. 2020). Moreover, the zonal flows on Saturn are symmetric enough to be described by rotation along cylinders up to midlatitudes (Militzer et al. 2019). In that case, equipotential theory still applies. We do not suggest that dynamic effects play a major role for the uncertainty in Saturn's shape and gravity field.

For Jupiter, the uncertainty in rotation rate is tiny, so that it is the influence of the winds of 2–4 km (Buccino et al. 2020) against which further effects must be compared. Such are the tidal bulges from the Galilean satellites. Nettelmann (2019) estimated a maximum elongation of 28 km in the direction of Io from static tidal response. The tidal flows around Jupiter are a dynamic perturbation and subject to Coriolis force (Idini & Stevenson 2021; Lai 2021). The flow and the Coriolis force acting upon it lead to dynamic contributions to the Love numbers k_nm. Juno measurements revealed a deviation by 1%–7% from the static k₂ value (Idini & Stevenson 2021). Approximating the corresponding shape deformation h₂ by h₂ = 1 + k₂ yields a tentative estimate of a (1%–7%) × 28 km ∼ (3–21) km additional shape deformation due to dynamic tidal response, which exceeds the wind effect. The possible importance of (periodic) perturbations on Jupiter's interior structure inference remains to be investigated.

Juno MWR data revealed that the ammonia abundance below the cloud level shows strong vertical and latitudinal variation (Guillot & Fletcher 2020). It is therefore possible that Jupiter's atmosphere is not everywhere well mixed where observations were taken. Consequently, the abundances and temperatures measured by the Galileo entry probe in a hot spot may not be representative of Jupiter's global atmosphere. On the other hand, analysis of Voyager 1 and 2 radio occultation data spanning a broad range of latitudes between 70° south and the equator yielded a 1 bar temperature of 165 K +/−5 K (Lindal et al. 1981), consistent with the Galileo measurement of 166.1 K in the hot spot. Present data therefore do not indicate that hot spots, in which deeper layers are exposed that appear brighter than surrounding regions at higher altitudes, were particularly cool regions, allowing us to suppose warmer global average temperatures. Rather, it is possible that the hot-spot temperature gradient is steeper than the global one since it is close to a dry adiabat (Seiff et al. 1998), whereas moist regions above the water cloud level may follow a less steep P–T profile (Kurosaki & Ikoma 2017), implying an even cooler interior below the cloud base. A colder and thus denser interior would strain the low-metallicity models even more.

The first step of the Legendre polynomial expansion of the centrifugal potential $Q=1/2\,{\omega }^{2}\,{r}^{2}{\sin }^{2}\vartheta$ is to write Q = − 1/3 ω² r²(P₀ − P₂(μ)). Its full expansion is obtained by replacing r with r_l (ϑ). Q(l, ϑ) can then be written in the form

$$\begin{eqnarray}\begin{array}{rcl}Q(l,\vartheta ) & = & -\,\displaystyle \frac{\mathrm{GM}}{{R}_{m}}\ {\left(\displaystyle \frac{l}{{R}_{m}}\right)}^{2}\displaystyle \sum _{k=0}^{O}{A}_{2k}^{(Q)}{P}_{2k}(\mu )\quad {\rm{with}}\\ {A}_{k}^{(Q)} & = & \displaystyle \frac{{m}_{\mathrm{rot}}}{3}\ \displaystyle \sum _{i}{c}_{i\,0\,k}.\end{array}\end{eqnarray} \tag{ A1 }$$

The gravitational potential $V({\bf{r}})=-G\int {d}^{3}r^{\prime} \ \rho ({\bf{r}}^{\prime} )/| {\bf{r}}^{\prime} -{\bf{r}}|$ can be separated into an external contribution D from the mass density interior to a sphere of radius r, i.e., to which r is exterior ($r\gt r^{\prime}$), and an internal contribution $D^{\prime}$ from the mass density exterior to a sphere of radius r, i.e., to which r is interior ($r\lt r^{\prime}$). V then reads

$$\begin{eqnarray}\begin{array}{rcl}V(r,\vartheta ) & = & -G\displaystyle \sum _{n=0}^{\infty }({r}^{-(n+1)}{D}_{2n}(r)\\ & & +\,{r}^{n}{D}_{2n}^{\prime} (r))\ {P}_{2n}(\mu )\end{array}\end{eqnarray} \tag{ A2 }$$

with

$$\begin{eqnarray}\begin{array}{rcl}{D}_{n} & = & {\displaystyle \int }_{r^{\prime} \lt r}{d}^{3}r^{\prime} \rho \ \ r{{\prime} }^{\,n}\ {P}_{n}(\mu ^{\prime} ),\\ {D}_{n}^{{\prime} } & = & {\displaystyle \int }_{r^{\prime} \gt r}{d}^{3}r^{\prime} \rho \ \ r{{\prime} }^{\,-(n+1)}\ {P}_{n}(\mu ^{\prime} ).\end{array}\end{eqnarray} \tag{ A3 }$$

Although this multipole expansion is valid only for spheres, in ToF, the radial coordinate r in Equations (A2) and (A3) is simply replaced by the nonspherical equipotential surface r_l (ϑ) and the interior/exterior criterion is transferred to l. In the CMS method, this expansion is also used, but the expression for the external potential is only applied to spheroids of level surface r_i (μ) at or interior (i ≥ j) to a point B of radial distance r_B that resides on a level surface r_j (μ). Since all spheroids share the same center but extend outward to different level surfaces r_i (μ), where i = 0 denotes the surface of the planet and i = N the center, a point B on r_j (μ) is also located exterior to the mass of the spheroids of index i < j but only as far as the radius r_B . This is taken care of in the CMS method by adding the external gravitational potential of the spheres of densities δ_i,i<j interior to B from the spheroids i < j. This improvement of the CMS method over the ToF method is still limited by the deformation and spacing of the spheroids. Rapid rotation, or dense spacing, could lead to an overlap of the sphere of radius r_j with the spheroid r_j−1(μ). Kong et al. (2013) developed the full solution to the Poisson equation and demonstrated that the CMS method converges as long as the flattening (R_eq/R_pol − 1) remains sufficiently small. Similarly, Hubbard et al. (2014) showed that ToF converges for sufficiently small flattening and toward the correct solution.

By replacing all powers of r by r_l (ϑ) in Equation (A2), one obtains

$$\begin{eqnarray}\begin{array}{rcl}V(l,\vartheta ) & = & -\displaystyle \frac{\mathrm{GM}}{{R}_{m}}{\left(\displaystyle \frac{l}{{R}_{m}}\right)}^{2}\\ & & \times \,\displaystyle \sum _{n=0}^{O}({\left(1+{\rm{\Sigma }}\right)}^{-(n+1)}{S}_{2n}(z)\\ & & +{\left(1+{\rm{\Sigma }}\right)}^{n}{S}_{2n}^{\prime} (z))\ {P}_{2n}(\mu )\end{array}\end{eqnarray} \tag{ A4 }$$

with z = l/R_m , z [0, 1], and, with the help of the transformation,

$$\begin{eqnarray*}&&{r}^{n}\ {dr}={dl}\ {r}^{n}\displaystyle \frac{{dr}}{{dl}}=\displaystyle \frac{1}{n+1}\ {dl}\ \displaystyle \frac{d}{{dl}}\ {r}^{n+1},\end{eqnarray*}$$

$$\begin{eqnarray}\begin{array}{rcl}{S}_{n}(z) & = & \displaystyle \frac{3}{2(n+3)}\ \displaystyle \frac{1}{{z}^{n+3}}{\displaystyle \int }_{0}^{z}{dz}^{\prime} \ \displaystyle \frac{\rho (z^{\prime} )}{\bar{\rho }}\ \displaystyle \frac{d}{{dz}^{\prime} }\\ & & \times \left[z{{\prime} }^{\,n+3}{\displaystyle \int }_{-1}^{1}d\mu ^{\prime} \ {\left(1+{\rm{\Sigma }}\right)}^{n+3}\ {P}_{n}(\mu ^{\prime} )\right],\end{array}\end{eqnarray} \tag{ A5 }$$

$$\begin{eqnarray}\begin{array}{rcl}{S}_{n}^{{\prime} }(z) & = & \displaystyle \frac{3}{2(2-n)}\ \displaystyle \frac{1}{{z}^{2-n}}{\displaystyle \int }_{z}^{1}{dz}^{\prime} \ \displaystyle \frac{\rho (z^{\prime} )}{\bar{\rho }}\ \displaystyle \frac{d}{{dz}^{\prime} }\\ & & \times \left[z{{\prime} }^{\,2-n}\ {\displaystyle \int }_{-1}^{1}d\mu ^{\prime} \ {\left(1+{\rm{\Sigma }}\right)}^{2-n}\ {P}_{n}(\mu ^{\prime} )\right].\end{array}\end{eqnarray} \tag{ A6 }$$

This can further be written as

$$\begin{eqnarray}&&{S}_{n}(z)=\displaystyle \frac{1}{{z}^{n+3}}{\int }_{0}^{z}{dz}^{\prime} \ \displaystyle \frac{\rho (z^{\prime} )}{\bar{\rho }}\ \displaystyle \frac{d}{{dz}^{\prime} }\left[z{{\prime} }^{\,n+3}\ {f}_{n}(z)\right],\end{eqnarray} \tag{ A7 }$$

$$\begin{eqnarray}&&{S}_{n,}^{\prime} (z)=\displaystyle \frac{1}{{z}^{2-n}}{\int }_{z}^{1}{dz}^{\prime} \ \displaystyle \frac{\rho (z^{\prime} )}{\bar{\rho }}\ \displaystyle \frac{d}{{dz}^{\prime} }\left[z{{\prime} }^{\,2-n}\ {f}_{n}^{{\prime} }(z)\right],\end{eqnarray} \tag{ A8 }$$

with

$$\begin{eqnarray}\begin{array}{rcl}{f}_{n}(z) & = & \displaystyle \frac{3}{2(n+3)}{\displaystyle \int }_{-1}^{1}d\mu \ {P}_{n}(\mu )\ {\left(1+{\rm{\Sigma }}\right)}^{n+3},\\ {f}_{n}^{{\prime} }(z) & = & \displaystyle \frac{3}{2(2-n)}{\displaystyle \int }_{-1}^{1}d\mu \ {P}_{n}(\mu )\ {\left(1+{\rm{\Sigma }}\right)}^{2-n}\quad (n\ne 2),\\ {f}_{2}^{{\prime} }(z) & = & \displaystyle \frac{3}{2}{\displaystyle \int }_{-1}^{1}d\mu \ {P}_{n}(\mu )\ \mathrm{ln}(1+{\rm{\Sigma }}).\end{array}\end{eqnarray} \tag{ A9 }$$

The ToF coefficients c_ink are of the form

$$\begin{eqnarray*}&&{c}_{\mathrm{ink}}={q}_{\mathrm{ink}}\ \displaystyle \prod _{j=1}^{O}{s}_{2j}^{{p}_{2j,\mathrm{ink}}},\end{eqnarray*}$$

where the q_ink are rational numbers and the exponents p_j are small natural numbers including 0.

Since the number of coefficients rises with the order of expansion faster than quadratically, it becomes impractical to write down all the coefficients. We present them in the form of five online tables. Tables tab_Sn and tab_Snp contain the coefficients c_ink and ${c}_{\mathrm{ink}}^{\prime}$ in front of the S_n and ${S}_{n}^{{\prime} }$ in the ${A}_{k}^{(V)}$, respectively, so that

$$\begin{eqnarray*}\begin{array}{rcl}{A}_{k}^{(V)} & = & \displaystyle \sum _{i=1}^{{N}_{0\,k}}{c}_{i\,0\,k}{S}_{0}+\displaystyle \sum _{i=1}^{{N}_{2\,k}}{c}_{i\,2\,k}{S}_{2}+\ldots \\ & & +\,\displaystyle \sum _{i=1}^{{N}_{14\,k}}{c}_{i\,14\,k}{S}_{14}+\ \displaystyle \sum _{i=1}^{{N^{\prime} }_{0\,k}}{c^{\prime} }_{i\,0\,k}{S}_{0}^{{\prime} }\\ & & +\,\displaystyle \sum _{i=1}^{{N^{\prime} }_{2\,k}}{c^{\prime} }_{i\,2\,k}{S}_{2}^{{\prime} }+\ldots +\displaystyle \sum _{i=1}^{{N^{\prime} }_{14\,k}}{c^{\prime} }_{i\,14\,k,}{S}_{14}^{{\prime} }.\end{array}\end{eqnarray*}$$

Table tab_m contains the summands in the ${A}_{k}^{(Q)}$ so that ${A}_{k}^{(Q)}={m}_{\mathrm{rot}}/3\ {\sum }_{i=1}{c}_{i\,0\,k}$ with ${m}_{\mathrm{rot}}={\omega }^{2}\,{R}_{m}^{3}/\mathrm{GM}$. Tables tab_fn and tab_fnp contain the summands of the functions f_n and ${f}_{n}^{{\prime} }$ so that f_n = ∑_i c_in , ${f}_{n}^{{\prime} }={\sum }_{i}{c^{\prime} }_{i\,n\,}$, respectively. In Table 2, we give an example of the read-in ASCII table tab_Sn. All five tables have the same format.

Table 2. Example of One of the Five Online-only ASCII Tables

n	k	Nnk
Order	p2	p4	p6	p8	p10	p12	p14	qink, i = 1—Nnk	Comment
0	0	24							n = 0, k = 0, next 24 rows
0	0	0	0	0	0	0	0	1.000000000000000e+00	Order = 0, c1 0 0 = 1
2	2	0	0	0	0	0	0	4.000000000000000e-01	Order= 2, ${c}_{2\,0\,0}=0.4\,{s}_{2}^{2}$
⋮
7	2	1	1	0	0	0	0	2.157842157842158e-01	Order = 7, ${c}_{24\,0\,0}={q}_{24\,0\,0}\ {s}_{2}^{2}{s}_{4}^{1}{s}_{6}^{1}$
...
4	8	17							n = 4, k = 8, next 17 rows
4	2	0	0	0	0	0	0	2.937062937062937e+00	Order = 4, ${c}_{1\,4\,8}={q}_{1\,4\,8}\ {s}_{2}^{2}$
⋮

Note. This example is for Table tab_Sn, which contains the coefficients in front of the S_n in the A_k . Since the ${A}_{2k}^{(Q)}$ have no dependence on index n, index n is set to 0 in Table tab_m. Since the f_n , ${f}_{n}^{{\prime} }$ have no dependence on index k, k is set to 0 in Tables tab_fn and tab_fn'.

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

The A_2k in Equation (3) are of the form 0 = A_2k (l) = −s_2k (l)S₀(l) + B if k ≠ 0. We rewrite this as an expression for direct, iterative computation of the figure functions in the form s_2k = B/S₀, where the functions B depend on the {s_2k } from the previous iteration step. We omit the summand − s_2k S₀ from Table tab_Sn.

To facilitate the application of our ToF7 tables by external users for their own planetary models, we share routines for read-in of the tables in Matlab, Python, and C++. For the latter variant, we also provide functions that can be used to easily access the coefficient values. The archived service routines and descriptions can be found at https://doi.org/10.6084/m9.figshare.16822252.

In the binomial expansion of (1 + x)^−m for m > 0,

$$\begin{eqnarray}&&{\left(1+x\right)}^{-m}=\sum _{i=0}^{\infty }\left(\begin{array}{c}-m\\ i\end{array}\right){x}^{i},\end{eqnarray} \tag{ A10 }$$

it is sufficient to expand to i = 7 because x = Σ and the minimum order of Σⁱ is i. The binomial expansion of ${\left(1+x\right)}^{m}$ for m > 0,

$$\begin{eqnarray}&&{\left(1+x\right)}^{m}=\sum _{i=0}^{m}\left(\begin{array}{c}m\\ i\end{array}\right){x}^{i},\end{eqnarray} \tag{ A11 }$$

is carried out to m ≤ 7. Products P_n P_m occurring in Σⁱ and in (1 + Σ)ⁱ P_j are expanded as ${\sum }_{k=0}^{n+m}{b}_{k}{P}_{k}$. It becomes evident that all terms can linearly be expanded in Legendre polynomials and that numbers in the expansion coefficients are rational numbers q = n _e /n _d . The natural numbers n _e , n _d can be represented exactly on a computer, although size limitations may apply. For ToF7, the vast majority of numbers could be decomposed into prime numbers that individually do not exceed 3 million. However, in rare cases this was not possible and larger prime numbers would have been required, perhaps indicating an error in the code used. In such a case, the given number is not decomposed into prime numbers. In any case, the enumerators and denominators are computed as exact numbers in all coefficients. They are cast to real numbers of 15 digits only for the purpose of printing the tables.

We calculate s₀ to the seventh order with the help of the defining integral Equation (2), which, with z = l/l₁, can be written as

$$\begin{eqnarray}&&\displaystyle \frac{4\pi }{3}{1}_{1}^{3}=2\pi \,{l}_{1}^{3}{\int }_{-1}^{1}d\mu \ {\int }_{0}^{1}{dz}\ {z}^{2}{\left(1+{\rm{\Sigma }}(z,\mu )\right)}^{3}.\end{eqnarray} \tag{ A12 }$$

Because of hemispheric symmetry and $1/3={\int }_{0}^{1}{dz}\ {z}^{2}$, the comparison of integrands yields

$$\begin{eqnarray}&&1={\int }_{0}^{1}d\mu \ {\left(1+{\rm{\Sigma }}\right)}^{3},\end{eqnarray} \tag{ A13 }$$

where ${\rm{\Sigma }}(z,\mu )={\sum }_{0}^{7}{s}_{2i}(z){P}_{2i}(\mu )$ and (1 + Σ)³ = 1 + 3Σ + 3Σ² + Σ³. We calculate Σ³ and safely remove all terms containing P₀ P_n P_m with n ≠ m or ${P}_{0}^{2}{P}_{n}$, since they would contribute nothing to the evaluated integral in Equation (A13),

$$\begin{eqnarray}\begin{array}{rcl}{{\rm{\Sigma }}}^{3} & = & {s}_{0}^{3}{P}_{0}^{3}+{s}_{2}^{3}{P}_{2}^{3}+{s}_{4}^{3}{P}_{4}^{3}\\ & & +\,3{s}_{0}{s}_{2}^{2}{P}_{0}{P}_{2}^{2}+3{s}_{2}^{2}{s}_{4}{P}_{2}^{2}{P}_{4}+3{s}_{2}^{2}{s}_{6}{P}_{2}^{2}{P}_{6}\\ & & +\,3{s}_{2}^{2}{s}_{8}{P}_{2}^{2}{P}_{8}+3{s}_{2}^{2}{s}_{10}{P}_{2}^{2}{P}_{10}\\ & & +\,3{s}_{0}{s}_{4}^{2}{P}_{0}{P}_{4}^{2}+3{s}_{2}{s}_{4}^{2}{P}_{2}{P}_{4}^{2}\\ & & +\,3{s}_{4}^{2}{s}_{6}{P}_{4}^{2}{P}_{6}+3{s}_{2}{s}_{6}^{2}{P}_{2}{P}_{6}^{2}\\ & & +\,6{s}_{2}{s}_{4}{s}_{6}{P}_{2}{P}_{4}{P}_{6}+6{s}_{2}{s}_{4}{s}_{8}{P}_{2}{P}_{4}{P}_{8}.\end{array}\end{eqnarray} \tag{ A14 }$$

Terms containing ${P}_{0}{P}_{n}^{2}$ yield a factor 1/(2n + 1) for the integral in Equation (A13). By ${({P}_{n}{P}_{m})}_{k}$ we denote the summand b_k P_k in the expansion of P_n × P_m . The other terms contribute

$$\begin{eqnarray}\begin{array}{lll}{P}_{0}^{3} & :\, & \to \,1\\ {P}_{2}^{3} & :\,{({P}_{2}{P}_{2})}_{2}{P}_{2}=\displaystyle \frac{2}{7}{P}_{2}^{2} & \to \,\displaystyle \frac{2}{5\cdot 7}\\ {P}_{4}^{3} & :\,{({P}_{4}{P}_{4})}_{4}{P}_{4}=\displaystyle \frac{2\cdot {3}^{4}}{7\cdot 11\cdot 13}{P}_{4}^{2} & \to \,\displaystyle \frac{2\cdot {3}^{2}}{7\cdot 11\cdot 13}\\ {P}_{2}^{2}{P}_{4} & :\,{({P}_{2}{P}_{2})}_{4}\ {P}_{4}=\displaystyle \frac{18}{35}{P}_{4}^{2} & \to \,\displaystyle \frac{2}{5\cdot 7}\\ {P}_{2}^{2}{P}_{6} & :\,{({P}_{2}{P}_{2})}_{6}\ {P}_{6} & \to \,0\\ {P}_{2}^{2}{P}_{8} & :\,{({P}_{2}{P}_{2})}_{8}\ {P}_{8} & \to \,0\\ {P}_{2}^{2}{P}_{10} & :\,{({P}_{2}{P}_{2})}_{10}\ {P}_{10} & \to \,0\\ {P}_{2}{P}_{4}^{2} & :\,{P}_{2}{({P}_{4}{P}_{4})}_{2}=\displaystyle \frac{{2}^{2}\cdot {5}^{2}}{7\cdot 11}{P}_{2}^{2} & \to \,\displaystyle \frac{{2}^{2}\cdot 5}{{3}^{2}\cdot 7\cdot 11}\\ {P}_{4}^{2}{P}_{6} & :\,{({P}_{4}{P}_{4})}_{6}{P}_{6}=\displaystyle \frac{{2}^{2}\cdot 5}{{3}^{2}\cdot 11}{P}_{6}^{2} & \to \,\displaystyle \frac{{2}^{2}\cdot 5}{{3}^{2}\cdot 11\cdot 13}\\ {P}_{2}{P}_{6}^{2} & :\,{P}_{2}{({P}_{6}{P}_{6})}_{2}=\displaystyle \frac{2\cdot 7}{11\cdot 13}{P}_{2}^{2} & \to \,\displaystyle \frac{2\cdot 7}{5\cdot 11\cdot 13}\\ {P}_{2}{P}_{4}{P}_{6} & :\,{({P}_{2}{P}_{4})}_{6}\ {P}_{6}=\displaystyle \frac{5}{11}{P}_{6}^{2} & \to \,\displaystyle \frac{5}{11\cdot 13}\\ {P}_{2}{P}_{4}{P}_{8} & :\,{({P}_{2}{P}_{4})}_{8}\ {P}_{8} & \to \,0.\end{array}\end{eqnarray} \tag{ A15 }$$

Similarly in Σ², terms P_n P_m yield a factor 1/(2n + 1) if n = m and 0 otherwise, so that

$$\begin{eqnarray*}&&{{\rm{\Sigma }}}^{2}={s}_{0}^{2}+\displaystyle \frac{1}{5}\,{s}_{2}^{2}+\displaystyle \frac{1}{{3}^{2}}\,{s}_{4}^{2}+\displaystyle \frac{1}{13}\,{s}_{6}^{2}.\end{eqnarray*}$$

Out of Σ¹ in the integral in Equation (A13), only s₀ P₀ survives. Now, s₀ can be expressed in terms of the s_n,n≥2, so that

$$\begin{eqnarray*}&&{s}_{0}:= {s}_{0}^{(1)}+{s}_{0}^{(2)}+\ldots +{s}_{0}^{(7)},\end{eqnarray*}$$

where O denotes the order of expansion. Sorting terms according to their order, we obtain (note the − sign)

$$\begin{eqnarray}\begin{array}{rcl}-{s}_{0}^{(1)} & = & 0\\ -{s}_{0}^{(2)} & = & \displaystyle \frac{1}{5}\ {s}_{2}^{2}\\ -{s}_{0}^{(3)} & = & \displaystyle \frac{2}{3\cdot 5\cdot 7}\ {s}_{2}^{3}\\ -{s}_{0}^{(4)} & = & \displaystyle \frac{1}{{3}^{2}}\ {s}_{4}^{2}+\displaystyle \frac{2}{5\cdot 7}\ {s}_{2}^{2}{s}_{4}\\ -{s}_{0}^{(5)} & = & \displaystyle \frac{2}{3\cdot {5}^{2}\cdot 7}\ {s}_{2}^{5}+\displaystyle \frac{{2}^{2}\cdot 5}{{3}^{2}\cdot 7\cdot 11}\ {s}_{2}{s}_{4}^{2}\\ -{s}_{0}^{(6)} & = & \displaystyle \frac{-127}{{3}^{2}\cdot {5}^{3}\cdot {7}^{2}}\ {s}_{2}^{6}+\displaystyle \frac{2}{{5}^{2}\cdot 7}\ {s}_{2}^{4}{s}_{4}+\displaystyle \frac{2\cdot 3}{7\cdot 11\cdot 13}\ {s}_{4}^{3}\\ & & +\,\displaystyle \frac{2\cdot 5}{11\cdot 13}\ {s}_{2}{s}_{4}{s}_{6}+\displaystyle \frac{1}{13}\ {s}_{6}^{2}\\ -{s}_{0}^{(7)} & = & \displaystyle \frac{2}{3\cdot {5}^{3}\cdot 7}\ {s}_{2}^{7}+\displaystyle \frac{2\cdot 41}{{3}^{3}\cdot 5\cdot 7\cdot 11}\ {s}_{2}^{3}{s}_{4}^{2}\\ & & +\,\displaystyle \frac{{2}^{3}}{3\cdot {5}^{2}\cdot 7}\ {s}_{2}^{5}{s}_{4}+\displaystyle \frac{2\cdot 7}{5\cdot 11\cdot 13}\ {s}_{2}{s}_{6}^{2}\\ & & +\,\displaystyle \frac{{2}^{2}\cdot 5}{{3}^{2}\cdot 11\cdot 13}\ {s}_{4}^{2}{s}_{6}.\end{array}\end{eqnarray} \tag{ A16 }$$