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

Interior modeling of Jupiter and Saturn has advanced to a state where thousands of models are generated that cover the uncertainty space of many parameters. This approach demands a fast method of computing their gravity field and shape. Moreover, the Cassini mission at Saturn and the ongoing Juno mission delivered gravitational harmonics up to J 12. Here we report the expansion of the theory of figures, which is a fast method for gravity field and shape computation, to the seventh order (ToF7), which allows for computation of up to J 14. We apply three different codes to compare the accuracy using polytropic models. We apply ToF7 to Jupiter and Saturn interior models in conjunction with CMS-19 H/He equation of state. For Jupiter, we find that J 6 is best matched by a transition from an He-depleted to He-enriched envelope at 2–2.5 Mbar. However, the atmospheric metallicity reaches 1 × solar only if the adiabat is perturbed toward lower densities, or if the surface temperature is enhanced by ∼14 K from the Galileo value. Our Saturn models imply a largely homogeneous-in-Z envelope at 1.5–4 × solar atop a small core. Perturbing the adiabat yields metallicity profiles with extended, heavy-element-enriched deep interior (diffuse core) out to 0.4 R Sat, as for Jupiter. Classical models with compact, dilute, or no core are possible as long as the deep interior is enriched in heavy elements. Including a thermal wind fitted to the observed wind speeds, representative Jupiter and Saturn models are consistent with all observed J n values.


Introduction
Since the era of the Voyager 1 and 2 gravity field determinations and shape measurements of the outer planets, only two methods have extensively been employed to calculate the shape and gravity field from interior models to compare with the data. These methods are the theory of figures (ToF; Zharkov & Trubitsyn 1978) and the Concentric Maclaurin Spheroid (CMS) method (Hubbard 2012(Hubbard , 2013. ToF has served that purpose before the advent of accurate gravity data from Juno at Jupiter and from the Cassini Grand Finale Tour at Saturn. Beforehand, only the gravitational moments J 2 , J 4 , and J 6 were measured, and the smallest given uncertainty in Jupiter's J 6 of 10% was still rather large (Jacobson 2003).
The low-order gravitational harmonics are important observables, as they constrain the density profile about midway into the planetary interior. They are expansion coefficients of the external planetary gravity field evaluated at a reference radius in the equatorial plane, R eq , which encompasses the planet's total mass. They are defined as integrals over density ρ(r) in the planet's interior, J M R d r r r P 1 cos , 1 n n n n eq 3 ò r where P n are the Legendre polynomials and ϑ is colatitude. Thanks to the Juno and Cassini missions, the observational accuracy in the even harmonics J 2n has seen significant improvement (Iess et al. 2018(Iess et al. , 2019Durante et al. 2020). For both Jupiter and Saturn, the uncertainties in the low-order harmonics reduced to a level that can be considered exact from the perspective of adjusting internal density distributions to reproduce the data. However, significant spread in the deep interior density distributions is still possible (see Movshovitz et al. 2020 for Saturn), as the sensitivity of the J 2n toward the center fades with r R n eq 2 ( ) (see Equation (1)). This spread is a residual uncertainty related to other causes such as the positioning of internal helium gradients due to uncertainty in the H-He phase diagram, the temperature profile in stably stratified regions, the H-He equations of state (EOSs), or the positioning of heavy-element gradients due to uncertainties in planet formation and evolution.
ToF to the fourth order (ToF4) has been deemed sufficiently accurate for computation of J 2 and J 4 usually used to constrain the density distribution (Nettelmann 2017). Higher-order moments beyond J 4 have been provided to a precision in J 6 to J 10 of better than (0.01-0.1) × 10 −6 for Jupiter (Durante et al. 2020) and (0.1-1) × 10 −6 for Saturn (Iess et al. 2019), but as the order increases, so does the influence by the zonal flows on the harmonics. At present, this is where the limitations of the ToF method become evident. ToF is an expansion method. An nth-order expansion (ToFn) allows us to compute up to J 2n and to an error of the order of q n rot 1 + , where q R GM rot 2 eq 3 w = is the ratio of centrifugal to gravitational force at the equatorial radius. The highest presented ToF order so far is 5 (Zharkov & Trubitsyn 1975). It has recently been applied to compute Saturn's J 2 -J 10 values (Ni 2020); however, its accuracy has not been validated yet.
Although the CMS method is also an expansion method, it can conveniently be carried out to the order of 15-20 or higher (Hubbard 2013). Therefore, it enables high-accuracy computation of the high-order J 2n up to the order of the measurements (Wahl et al. 2017b;Militzer et al. 2019). The CMS method provides further advantages, such as its expansion to 3D to account for tidal shape and gravity field perturbations (Wahl et al. 2017a) and brevity in its formulation (Hubbard 2013). Its only drawback is that the CMS method goes along with high computational cost even in its accelerated version . This is because the CMS method explicitly solves for the 2D planetary shape not only by taking the sum over radial spheroids but also by integrating over latitude. Even if making use of Gaussian quadrature, typically several tens of angular grid points are required. If the integrand were a polynomial, only N lat = n + 1 grid points would be required to evaluate the integrals over P 2n , which for J 12 (n = 6) amounts to only N lat = 7, or even 4 points when accounting for hemispheric symmetry. However, the integrands are functions of the nonpolynomial shape itself. In practice, 48 grid points (Wahl et al. 2017a) are used. Obtaining the shape to sufficient accuracy at these grid points is the most timeconsuming part of the CMS method. In contrast, ToF solves for the shape explicitly only in the equatorial plane, while the shape at higher latitudes is obtained by spherical harmonics expansion, and the required precision of the shape is only the one wanted for the J 2n . One may thus see a benefit in using ToF for computation of the high-order moments.
Here we introduce ToF7 tables, which allow one to calculate up to J 14 . In Section 2, we give an overview of the ToF method, while for further details we refer to the Appendix. In Section 3, we assess the accuracy of the ToF method by comparing to the analytic n = 1 polytrope solution. In Section 4, we apply the new tables to Jupiter models, and in Section 5 we apply them to Saturn. In Section 6, we connect representative interior models to thermal wind models to predict the wind decay depth profiles. Observing notoriously low atmospheric metallicities of our Jupiter models, we discuss further influences in Section 7. Section 8 concludes the main body of the paper. In Appendix A.3 we introduce the ToF7 tables for public usage.

Theory of Figures
The ToF is described in Zharkov & Trubitsyn (1978), and the coefficients up to the third order are presented therein. Nettelmann (2017) followed their notation and calculated the fourth-order coefficients. We note that fifth-order coefficients were presented in Zharkov & Trubitsyn (1975) and adopted by Ni (2020) for application to Saturn. Building on the work of Nettelmann (2017), we here conduct the expansion of ToF to the seventh order, meaning that the even harmonics up to J 14 can be calculated.
Both the ToF and CMS methods assume that surfaces of equal potential U exist on which density and pressure are constant. One can show that this assumption holds for planets in hydrostatic equilibrium that rotate along cylinders, e.g., if their rotation rate can be expressed as s e w w = w ( )   with axis distance s. Rigid rotation and no rotation meet this condition. For rotation along cylinders, the odd harmonics J 2n+1 disappear. However, the Juno measurements at Jupiter (Iess et al. 2018;Durante et al. 2020) and Cassini Grand Finale at Saturn (Iess et al. 2019) revealed that these planets' odd harmonics are nonzero. Instead, they are of the order of 0.1 × 10 6 , comparable to the values of J 10 , J 12 , and Saturn's uncertainty in J 6 . Based on the commonly used approach of using the thermal wind equation (TWE) to infer the density anomalies (Kaspi et al. 2010;Kaspi 2013), the depth of the windinduced deviation from cylinder rotation has been inferred to be about 3000 km ( ∼ 0.035 R Jup ) in Jupiter (Kaspi et al. 2018) and 9000 km ( ∼ 0.14R Sat ) in Saturn (Galanti et al. 2019). These results are consistent with the tangent-cylinder model of Dietrich et al. (2021), which goes beyond the TWE simplification by including not only the wind-induced density perturbation but also the associated gravitational perturbation, which has long been argued to be significant (Zhang et al. 2015). Taking into account also constraints from the observed secular variation of the magnetic field on deep flows (Moore et al. 2019) suggests a somewhat steeper decay function for the winds, with the zonal flow extending inward on cylinders almost barotropicaly to a depth of about 2000 km on Jupiter and 8000 km on Saturn, and then the winds decay abruptly within the next 1000 km . While the nonasymmetric gravity field is important for these width depth issues, here we have to neglect the asymmetries, as otherwise neither the ToF nor the CMS method could be applied.
In the absence of tides, the problem at hand is axisymmetric and thus 2D: r, ϑ. In ToF, the description is further reduced to 1D by introducing the mean radius coordinate l. Spheres of radius l are defined by the condition that they enclose the same volume as the equipotential surface r l (ϑ), On the surface of the planet, with R(π/2) = R eq . In ToF, the potential is thus constant on spheres. Both the total potential U(l) and the axisymmetric 2D shape r l (ϑ) are expanded into Legendre polynomials. One can write while the shape of an equipotential surface, also called level surface, is given by where the s 2k (l) are the figure functions. The condition that U(l) is constant on spheres of radius l implies that A 2k = 0 for k > 0. These expansions are carried out up to an order O. In the absence of tides, U is a superposition of only the gravitational potential V and the centrifugal potential Q so that By definition, the gravitational harmonics J 2n can be obtained in the ToF as where the integrals S 2i , not to be confused with the figure functions s 2i , are given by Equation (A7) in the Appendix. However, using Equation (A7) for the S n and the ToF expansion coefficients to calculate the functions f n on which the S n depend (see Equations (A7) and (A9) in the Appendix) implies that only information on the equatorial radius r l (π/2) enters the computation of the J 2n , while information on the full shape r(l, θ) is reduced to the order of the expansion, that is, up to P 14 . An alternative method is to calculate the integrals over latitude explicitly. In Section 3.3 we compare both methods. For details on how the ToF coefficients are computed and for an example of the machine-readable ASCII tables that contain their values for public usage, see Appendix A.3. Moreover, to facilitate the application of our ToF7 tables, we share computer routines for read-in of the tables and documentation at https://doi.org/10.6084/m9. figshare.16822252.

Validation against the n = 1 Polytrope
The n = 1 polytropic planet is specified by a number of conditions. First, the polytropic EOS P = K ρ 2 . Furthermore, the gravity field of the rotating polytrope depends on the values of q rot , equatorial radius R eq , and planet mass M. The density profile ρ(r) is not known in advance but is obtained from solution of the equation of hydrostatic equilibrium, ∇P/ρ = ∇U. In ToF, the radial coordinate r is transformed to the level surface l, and the equation of hydrostatic equilibrium reduces to dP/dl = ρdU(l). The internal m-l relation is obtained by integrating the equation of mass conservation, m l dl l l 4 l 0 2 ò p r = ( ) ( ) . The latter is a source of numerical inaccuracy. We employ three different codes to compute the solution to the rotating polytrope. All polytropic models use q rot = 0.089195487 and GM = 126686536.1 × 10 9 km 3 s −1 as in Wisdom & Hubbard (2016).
Before we compare the results of our application of three different codes and different orders of expansion of ToF to the analytic, Bessel-function-based method of Wisdom & Hubbard (2016) in Figure 1, we describe each of the three employed methods in Sections 3.1-3.3.

Polytrope with MOGROP
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 2 ρ(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.

Polytrope with TOF-PLANET
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 −11 m 3 kg −1 s −2 .) 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.

Polytrope with CEPAM
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ρ 2 , 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 z z z sin c r r p p = ( ) . 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 3 , 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 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 z P z K j j r = ( ) ( ) . 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),    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 2 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 4 -J 10 , where the numerical accuracy for Equation (7) is about 1-2 orders of magnitude better than that for Equation (5).

Comparison of the J 2n of the Uniformly Rotating Polytrope
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 3 , the numerical accuracy in J 4 -J 10 falls below the current observational uncertainty (Juno D20) reported in Durante et al. (2020). Moreover, the accuracy in all the harmonics J 2 -J 10 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 2 and J 4 . For these low-order harmonics, ToF7 versus ToF5 provides a negligible improvement in accuracy. The situation changes with J 6 . Here, CEPAM with ToF5 levels off at a relative uncertainty of ∼ 5 × 10 −5 , 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 6 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 8 and higher, even with MOGROP, where J 8 becomes an order of magnitude better than the observational one, 2.5 orders of mag in J 10 , and 3 orders of mag in J 12 . We conclude that ToF7 is sufficiently accurate to address the influence of the winds on J 6 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 4 -J 8 with ToF4 after convergence with grid point number N agree, indicating that the remaining errors ΔJ n /J n of 2.5 × 10 −4 in J 4 , 10 −2 in J 6 , and 10 −1 in J 8 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 10 and J 12 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 2 and J 4 , and it is still better than compared to using 10 5 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.

Jupiter Models
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 1 = 0.238 as measured by the Galileo entry probe. Layers 2 and 3 have the same helium abundance (Y 2 = Y 3 ), 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 12 are considered, which are typical pressures where the Jupiter Table 1 Comparison of the J 2n Obtained from Equation (7) and from Equation (5) Method J 2 × 10 6 J 4 × 10 6 J 6 × 10 6 J 8 × 10 6 J 10 × 10 6 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) 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 precompressed to 2-4 GPa in diamond anvil cells and further shockcompressed 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 1 = Z 2 ). Finally, between layers 2 and 3, the heavy-element mass fraction is allowed to change. We use either a constant Z 3 value, implying a jump in Z at a transition pressure P 23 , or a Gaussian-Z 3 profile that starts with Z(P 23 ) = Z 2 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 2 and J 4 are Z 1 and Z 3 or Z 3,max .
We employ the CMS-2019 EOSs for H and He ) 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 ) after adding a compositiondependent entropy of mixing term s X X , mix H He 2 ( ) . 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/  Guillot et al. (2018;yellow), models with CMS-19 H/He EOS and P 12 = 1, 3, 4 Mbar (triangles), P 12 = 2 Mbar and P 23 varied from 5 to 20 Mbar, and constant Z 3 (diamonds) or Gaussian Z 3 (circles). All J 2n values are multiplied by 10 6 .
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.

Results for Jupiter's Even Harmonics
In Figure 2 we show the even J 2 -J 10 values from rigidly rotating Jupiter models. Models adjusted to the Juno observations of J 2 and J 4 (Durante et al. 2020) are shown in bluish color, while models adjusted to the wind-corrected J 2 and J 4 values by the corrections of Kaspi et al. (2018) applied to the J 2 , J 4 values of Durante et al. (2020) are shown in reddish color.
Due to imperfect fit to the J 2 , J 4 values, the scatter in model J 2 and J 4 values is larger than the observational uncertainty. For J 4 , the scatter ΔJ 4 /J 4 is about ±2 × 10 −4 and of the same size as the relative uncertainty due to using ToF7 in the MOGROP code, while for J 2 the latter relative uncertainty is with 7 × 10 −5 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 2 model values of ±3.4 × 10 −5 and a wider scatter in J 4 of ±10 −3 relative deviations. They found that, nevertheless, the high-order moments J 8 versus J 6 and J 10 versus J 8 were strictly confined to a straight line. We confirm that behavior. Notably, the model |J 8 | values are higher than the observed value, and a trend in that direction is also seen for J 10 , although the model J 10 values are still within the 3σ observational uncertainty. Wind models assuming rotation along cylinders indeed predict a slight decrease in |J 8 | and J 10 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 6 , 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 6 but also yield a scatter around the observed value. Hence, we conclude that Jupiter's J 6 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 3 and an abrupt He-poor/He-rich transition at P 12 = 1-2 Mbar yield J 6 = ( 34.11-34.18) x10 −6 , slightly below the observed value, while models with that transition deeper inside at P 12 = 2.5-3 Mbar yield J 6 = ( 34.23-34.28) x 10 −6 , slightly above the observed value. Constant-Z 3 profiles and P 12 = 2 Mbar stretch from J 6 = ( 34.18 to 34.28) times 10 −6 around the observed value (34.200 7 ± 0.0067) x 10 −6 upon shifting P 23 from 5 Mbar deeper down to 18 Mbar. Taken at face value, Jupiter's observed J 6 value indicates that the He-depleted/He-enriched transition occurs at around 2-2.5 Mbar.

Z Profiles for Jupiter
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.  (2016) (orange), the Z(m) profiles of the formation models A and C (dark blue) of Lozovsky et al. (2017) at the final stage of mass accretion, after settling but before possible homogenization by mixing, and Z(m) profiles of Müller et al. (2020) for their envelope accretion models assuming a hot-compact or a cold-extended state at the onset of gas accretion (green).

Unmodified H/He Adiabat
All our models with CMS-19 EOSs that fit J 2 and J 4 have negative Z 1 values between −0.005 and −0.020 (−0.33 to −1.33 × solar). This is consistent with Hubbard & Militzer (2016), who obtained −0.6 M E of heavy elements in the molecular region for their MH13 EOS-based model DFT-MD7.15, which, with J 4 = 587 x 10 −6 , is the one that comes closest to the Juno value of -586.61 x10 −6 . For a conservative estimate of their Z 1 value we take 1 Mbar, the entry of their Jovian adiabat into the H/He demixing boundary of Morales et al. (2013), or 2 Mbar, the pressure medium in their H/Herain region. With a corresponding molecular envelope mass of ∼30 M E and ∼53 M E , respectively, we obtain Z 1 between −0.011 and −0.02 for model . In contrast, Debras & Chabrier (2019) found a variety of models for nonnegative atmospheric Z values of 1 × Z Gal = 0.0167 using CMS-19 H/He EOSs. We cannot reproduce the results of Debras & Chabrier (2019) quantitatively.
The deeper the layer boundary for heavy elements is placed, the higher will the deep interior heavy-element abundance become, and the smaller the core mass (Nettelmann et al. 2012). With CMS-19 EOSs, the response of M core to P 23 is comparably weak, so that P 23 can be placed as deep as 20 Mbar before the compact core disappears. The thick blue model in Figure 3 is for P 23 = 21 Mbar and has a core mass of only 0.25 M E .
The nonexclusive compact core mass values of our models range from 0.2 to 4.8 M E for P 12 = 2 Mbar and Gaussian-Z envelopes, from 0.6 to 6.0 M E for constant-Z envelopes, and from 1.2 to 3.8 M E for wind-corrected models.
For Gaussian-Z deep envelopes, Z 3,max can become quite large toward the center. We obtained Z 3,max values up to 0.5, although larger values may be possible if the maximum of the Gaussian curve is placed at the center, while we placed it slightly off at 38 Mbar.
The mass of heavy elements in the deep interior below the negative-Z envelope amounts to 7.5-10.1 M E . Assuming a 1 × solar instead of negative-Z envelope would add another 3.8M E of heavy elements. A total of 11.3-13.9 M E of heavy elements is consistent with the Jupiter core accretion formation models A and C of Lozovsky et al. (2017). These models assume solid surface densities of 6 and 10 g cm −2 and planetesimal sizes of 100 and 1 km, respectively. Lozovsky et al.  Figure 3, models A and C are shown after settling of heavy elements but before possible convective mixing. Settling takes place if the partial pressure of ablated incoming material exceeds its vapor pressure. The resulting Z(m) profiles after formation resemble our interior models with Gaussian Z 3 , although the Z gradient in the post-formation models begins farther out at ∼0.5M J than at ∼0.2-0.3M J as in our models. On the other hand, a shallow, primordial compositional gradient has been found to erode and to be erased in present Jupiter if vigorous convection takes place in the envelope (Müller et al. 2020), while a steep compositional gradient may still persist within 0.2M J . The present-state models for Jupiter of Müller et al. (2020) are similar to our models with either Gaussian Z 3 or constant Z 3 when the heavy-element-enriched deep interior (or dilute core) is assumed to begin deep inside at > 15 Mbar, except that our models underestimate the outer envelope metallicity while the evolution models of Müller et al. (2020) overestimate it, as they yield too small a present-day radius.

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 4 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 .
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(Guillot et al. , 2004, leading to subadiabatic stable stratification according to the Schwarzschild criterion. Superadiabatic gradients are predicted in a Ledoux-stable, inhomogeneous medium of upwarddecreasing 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 ). 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 . 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 . 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 threeparameter study where we vary the maximum deviation max dr , 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 P log dr ( ). 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 4 | and |J 4 |/J 2 ratio. We find P start 30 GPa necessary in order to have a noticeable influence on J 4 . Conversely, higher P end values lead to a stronger reduction on J 2 . The question we ask is, for what values of dr, P start , and P end can a model be found with a 1 × solar homogeneous Z?
For a homogeneous, unperturbed adiabat at Y 1 = 0.238 and Y 2 adjusted to meet Y = 0.27, both |J 4 | and J 2 turn out to be significantly too large. This is in contrast to the result by Debras & Chabrier (2019), who could match J 2 at 1 × Z Gal . Thus, we need P start to be sufficiently low for J 4 and P end sufficiently high for J 2 . 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 2 -J 4 and J 6 -J 4 space is shown in Figures 4(a) and (b). Notably, among the wide spread of intermediate models in the J 2 -J 4 -J 6 space, the one model (red cross) that meets the Juno J 2 , J 4 values yields J 6 = 34.20 x 10 −6 , 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 6 .

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 2 , J 4 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 6 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 4 are small and not considered here. Notably, for T 1 bar = 180 K and at the transition pressure between the Hepoor and He-rich regions at P 12 = 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 6 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 1 rises with P 12 .
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.

Saturn Models
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 2 , J 4 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 ), using the Cassini Grand Finale gravity and the same shape data, and with the rotation rate of 10:33:38 hr m s m s 1 19 1 52 -+ 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 1 = 0.06 ± 0.05 stems from a combined Voyager radio occultation and infrared spectra analysis (Conrath et al. 1984), while the highest estimate Y 1 = 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 1 = 0.075-0.13 from Cassini infrared remote sensing (Achterberg & Flasar 2020) to Y 1 = 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 1 = 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 1 = 0.14 and Y 1 = 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 1 = 0.18 ± 0.07, Militzer et al. (2019) used Y 1 = 0.18-0.26, and Ni (2020) used Y 1 = 0.12-0.23. Lower Y 1 values yield higher atmospheric and higher maximum deep interior metallicities 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 12 = 2-4 Mbar deep within the He-rain region, for which predicted onset pressures in present Saturn range from 0.8 Mbar  to 2 Mbar (Mankovich & Fortney 2021). Below P 12 , we keep the He/H ratio constant with depth.
Typical core-envelope pressures are around 15 Mbar. We vary P 23 between 4 and 7 Mbar and let the Gaussian-Z 3 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.  (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° . 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 6 , 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 6 values and a 5.5% effect, while Galanti et al. (2019), who allow for a much wider scatter in model J 4 values of ±40 × 10 −6 , obtain UR model J 6 values up to 87 × 10 −6 , which encompasses the observed value. However, the mean of their distribution for fast, uniform rotation lies at 82 × 10 −6 , implying a 5.5% influence of the winds on J 6 , consistent with this work. The strong influence of the winds seen in J 6 suggests that also J 4 and J 2 are affected by the winds. In Section 6, we investigate whether the observed winds are consistent with a smaller influence on J 6 than found in previous work (Iess et al. 2019).

Z and ρ Profiles for Saturn
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 1 = 0.02-0.06 (1.5-4 × solar) for nominal He abundances Y 1 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),  (Iess et al. 2019;golden). This work's models: blue upward-pointing triangles: adiabatic, Y 1 = 0.14, Gaussian Z 3 ; blue downward-pointing triangles: adiabatic, Y 1 = 0.14, constant Z 3 ; blue circles: Y 1 = 0.18, superadiabatic; red circles: modified H/He adiabat, dilute core, Y 1 = 0.18; red diamonds: modified H/He EOS and Y 1 between 0.16 and 0.06, dilute core, Gaussian Z 3 . 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 2 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 −3 (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.
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 heavyelement 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.

Including the Zonal Wind Profiles
Our Jupiter models with a modified or unmodified H/He adiabat allow for a largely homogeneous interior down to 0.4 R J with a small rock core ( Figure 3) and for J 6 unaffected by the winds (Figure 2). Our Saturn models allow for larger static J 6 values in the range (82.5-84.0) × 10 −6 than previous work (81.8 × 10 −6 ); see Figure 6.
Here we investigate whether such models for Jupiter and Saturn are consistent with the observed wind profiles and odd and even J n values. We pick one representative Jupiter model (unmodified adiabat, P 12 = 2, P 23 = 20 Mbar, compact core of mass 1.26 M E ) and the Saturn model highlighted in red in Figure 7 (modified H/He adiabat, dilute core, Y = 0.14).
The wind modeling approach is the same as described in , except that the Saturn rotation period to which the wind speeds refer is adjusted to the value of the interior model, 10h32m45s. We use either only the gravity data to constrain the wind decay depth profiles, as in Kaspi et al. (2018) and Galanti et al. (2019), or both the gravity and magnetic field data combined (grav+MHD), as in . As an uncertainty in the static J 2n values we take the Tof7 errors produced by MOGROP for N = 4000.
With this approach, we are able to find fits that reproduce all the observed J n within the observational uncertainties. The wind profiles and decay depths are shown in Figure 8. This implies that extending the wind profiles, roughly as they appear at the cloud level, gives a good match to the difference between the Juno and Cassini measurements and our preferred models. Nonetheless, there is enough freedom in these solutions that other wind profiles with small shifts to the wind profiles can give fitting solutions as well . For these wind profiles, grav and grav+MHD yield similar solutions. Jupiter's wind profile is slightly less well matched (red and blue lines more strongly deviate from the observed profile (gray) than does the black dashed line), while for Saturn, the shoulders at 20°-40°latitude are somewhat better matched than in previous work (Galanti et al. 2019;Iess et al. 2019). The wind decay depths for grav only are slightly steeper than for previous interior models and thus closer to the grav+MHD solutions.
For each of Jupiter and Saturn, we picked only one specific interior model to calculate the wind contribution and optimize for the agreement with the observed wind velocities and gravitational harmonics. The fact that these two interior models allowed for solutions within the observed values and the ToF7/ MOGROP uncertainty strongly suggests that there are further interior models for which such a fit can be obtained. This means that the joint interior and wind solutions are not unique, given the uncertainties we allowed for. In addition, alternative interior models which fit all the J n when combined with a wind model may be possible for different EOSs and wind models, such as the MH13 H/He EOS and wind models that account for the oblate shape (Cao & Stevenson 2017) or solve for the gravo-thermal wind equation that accounts for the dynamic self-gravity of the flow (TGWE; Kong et al. 2018;Wicht et al. 2020). We note that Galanti et al. (2017) find that these modified wind models introduce corrections that are an order of magnitude smaller for most J n , while Dietrich et al. (2021) obtain corrections of ∼60% and ∼20% for J 3 and J 5 , respectively, when including the dynamic self-gravity for polytropic models and an additional correction of ∼40% and ∼10% when accounting for Jupiter-model specific background density and gravity profiles.
Internal flow structures, which are decoupled from the observed cloud-level winds, can also be found to fit the J n (Kaspi et al. 2018;Kong et al. 2018) and thus lead to nonuniqueness of the solutions (Kong et al. 2018). Here we conclude for nonuniqueness because of uncertainties in the interior models, the high-order J n to be fitted, and the wind profile.

Discussion
Our Jupiter and Saturn models exhibit a strong trend toward low envelope metallicities that extend deep into Saturn's interior to ∼ 0.4 R S (Figure 7) or are negative in Jupiter (Figure 3). We have attributed these model properties to possible uncertainties in the H/He EOS; however, one may think of further processes.

Z and the Adiabatic P-T Profile
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 Figure 8. Wind profiles (top) for Jupiter (left) and Saturn (right) using the constraints from observed gravity data only (blue) or also from MHD (red). The black dashed lines for Jupiter are taken from Kaspi et al. (2018;grav only), while for Saturn they are adjusted from grav only). Solutions for grav +MHD from that previous work are not shown since the solutions are close to the ones from this work. The gray lines in the top panels show the observed profiles (Tollefson et al. 2017). 10-100 GPa region relevant for J 2 and J 4 . 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.

H/He Demixing?
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 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.

Deep Internal Flows in Jupiter?
In contrast, in order to lift Jupiter's atmospheric metallicities substantially, a much stronger and retrograde deep wind in the interior where J 2 and J 4 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 −1 ) wind in Jupiter's deep interior.

Uncertainty in the Shape due to Dynamical Effects?
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,  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 6 in J 2 and 0.2 × 10 6 for J 4 . While small, this effect on J 4 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 −5 , 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 . 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 2 value (Idini & Stevenson 2021). Approximating the corresponding shape deformation h 2 by h 2 = 1 + k 2 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.

A Cold Hot Spot?
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.

Conclusions
We present the expansion of the ToF (Zharkov & Trubitsyn 1978) from formerly fifth order (Zharkov & Trubitsyn 1975) to seventh order. The coefficients are available in the form of five read-in online tables and allow the computation of the even gravitational harmonics J 2 -J 14 and the shape of a rotating fluid body in hydrostatic equilibrium.
We estimate the numerical accuracy of the ToF method carried out to fourth (Nettelmann 2017), fifth (Zharkov & Trubitsyn 1975), and seventh (this work) order by comparing to the analytic Bessel solution of Wisdom & Hubbard (2016) for the rotating n = 1 polytrope and by using three different codes. We find that the CEPAM code (Guillot & Morel 1995) with ToF5 (Ni 2020) has a superior performance in regard to the accuracy in J 2 , J 4 , and J 6 , while for J 8 and J 10 the MOGROP code with ToF7 reaches a similar degree of accuracy for a practical number of radial grid points of a few thousand, although in J 10 the error changes sign between both variants. The accuracy in J 8 , J 10 , and J 12 falls by, respectively, 1, 2, and 3 orders of magnitude below the current 3 × 1σ formal uncertainty of the observational gravity data "halfway through the Juno mission" (Durante et al. 2020). We also apply the CMS-2019 H/He EOS of Chabrier et al. (2019) to interior models of Jupiter and Saturn.
For Jupiter, the high-order J n of the Jupiter models fall along the same line in J n -J n+2 space as in previous work, regardless of detailed model assumptions and the H/He EOS used. We find that J 6 stands out in that it is neither adjusted, as J 2 and J 4 are, nor insensitive to model assumptions, as the J n for n 8 are. We match Jupiter's observed J 6 value by placing the transition pressure between an outer, He-depleted envelope and an inner, He-enriched envelope at P 12 = 2-2.5 Mbar. Transition pressures farther out lead to lower J 6 values, while deeper transitions result in higher J 6 values. The same behavior but with a weaker amplitude is seen for the transition pressure of heavy elements, which we place between ∼5 and 20 Mbar. Gaussian-Z profiles underneath can lead to high metallicities of up to Z = 0.5 at the compact core-mantle boundary. However, the atmospheric heavy-element abundance, represented by an EOS of water, always stays negative (∼−1 × solar) if the adiabat is defined by the 1 bar temperature of 166.1 K as measured by Galileo and extended downward. Alternatively, we set Z 1 to 1 × solar, the 1σ lower limit of the equatorial water abundance measured by Juno , and perturb the adiabat to fit J 2 and J 4 . Such an optimized adiabat was found for a perturbation between 26 and 150 GPa and has a maximum density decrease of 12.6% at a midpoint of 63 GPa. Higher internal temperatures help to decrease the internal density as well. For T 1bar = 180 K and deep transition pressure P trans, He = 6 Mbar we obtain 1× solar metallicity without H/ He EOS modification.
Our Saturn models with CMS-19 H/He EOS are characterized by a few times solar envelope that extends deep down to < 0.4 R Sat and requires a compact core. Its density of ∼20 g cm −3 is higher than the most likely central densities of Saturn that match the gravity field (Movshovitz et al. 2020), which in turn agree with density distributions of a largely stably stratified deep interior with a dilute core (Mankovich & Fuller 2021). By applying the Jupiter-optimized perturbation along the adiabat to Saturn, we are able to obtain density distributions with a dilute core of 30%-40% H/He that reaches out to ∼ 0.4R Sat in the core. This moves the solution in the direction of density distribution constrained by seismic data. Our models suggest that an inhomogeneous central region out to ∼0.6 R Sat (Mankovich & Fuller 2021) is due to both a dilute core and rained-down helium.
Overall, our Jupiter and Saturn models exhibit a strong trend toward low envelope metallicities that extend deep into Saturn's interior to ∼ 0.4 R S (Figure 7) or are negative in Jupiter (Figure 3). We have attributed these model properties to possible uncertainties in the H/He EOS. However, further processes one may think of and that certainly are at play are estimated to be too minor to solve that issue. This work demonstrates that our understanding of the internal heavy-element distribution of Jupiter and Saturn strongly depends on the properties of H and He. We conclude that part of the difficulties of obtaining Jupiter and Saturn models that are consistent with all observational constraints still lies in our imperfect understanding of the material properties. We therefore suggest that further measurements and calculations of the behavior of materials at planetary conditions could improve our understanding of the gas giants. This, in return, will also reflect on the characterization of gaseous planets orbiting other stars.
We thank the IWG members of the Juno Team for discussions. N.N. and J.J.F. acknowledge support through NASA's Juno Participating Scientist Program under grant 80NSSC19K1286. We thank the two anonymous reviewers for the constructive reports and insightful comments.

Appendix A Technical Notes on ToF Coefficient Computations
In the following, we abbreviate the term in brackets in Equation (4) as (1 + Σ) and write cos m J = .
Since any point r = (r, j, ϑ) in and near the planet can be associated with an equipotential surface, we can replace any dependence f (r, ϑ) by f (l, ϑ) using Equation (4).

A.1. Centrifugal Potential Q
The first step of the Legendre polynomial expansion of the centrifugal potential Q r 1 2 sin 2 2 2 w J = is to write Q = − 1/3 ω 2 r 2 (P 0 − P 2 (μ)). Its full expansion is obtained by replacing r with r l (ϑ). Q(l, ϑ) can then be written in the form The gravitational potential V G d r r r r r 3 ò 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 r > ¢), and an internal contribution D¢ from the mass density exterior to a sphere of radius r, i.e., to which r is interior (r r < ¢). V then reads 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 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 ink ¢ in front of the S n and S n ¢ in the A k V ( ) ,  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 k Q 2 ( ) have no dependence on index n, index n is set to 0 in Table tab_m. Since the f n , f n ¢ have no dependence on index k, k is set to 0 in Tables tab_fn and tab_fn'.  . Tables tab_fn and tab_fnp contain the summands of the functions f n and f n ¢ so that f n = ∑ i c in , f c n i 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. The A 2k in Equation (3) are of the form 0 = A 2k (l) = −s 2k (l)S 0 (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 0 , where the functions B depend on the {s 2k } from the previous iteration step. We omit the summand − s 2k S 0 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.

A.4. Powers of the Radius
In the binomial expansion of (1 + x) −m for m > 0, is carried out to m 7. Products P n P m occurring in Σ i and in (1 + Σ) i P j are expanded as b P k n m k k 0 å = + . 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 denote the summand b k P k in the expansion of P n × P m . The other terms contribute P P PP P P P PP P P P P P P P P P P P P P P P P P P P P P P P P P P P P P P P P P P P P P P P P P P P P P P P P P P P P P P · · · · · · · · · · · · · · · · · · · · · ( )