{'fraction_non_alphanumeric': 0.03336276129179245, 'fraction_numerical': 0.021487156120066072, 'mean_word_length': 4.262570251030349, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 15, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A complete study of the precision of the concentric MacLaurin spheroid method to calculate Jupiter's gravitational moments

A few years ago, Hubbard (2012, 2013) presented an elegant, non-perturbative method, called concentric MacLaurin spheroid (CMS), to calculate with very high accuracy the gravitational moments of a rotating fluid body following a barotropic pressure-density relationship. Having such an accurate method is of great importance for taking full advantage of the Juno mission, and its extremely precise determination of Jupiter gravitational moments, to better constrain the internal structure of the planet. Recently, several authors have applied this method to the Juno mission with 512 spheroids linearly spaced in altitude. We demonstrate in this paper that such calculations lead to errors larger than Juno's error bars, invalidating the aforederived Jupiter models at the level required by Juno's precision. We show that, in order to fulfill Juno's observational constraints, at least 1500 spheroids must be used with a cubic, square or exponential repartition, the most reliable solutions. When using a realistic equation of state instead of a polytrope, we highlight the necessity to properly describe the outermost layers to derive an accurate boundary condition, excluding in particular a zero pressure outer condition. Providing all these constraints are fulfilled, the CMS method can indeed be used to derive Jupiter models within Juno's present observational constraints. However, we show that the treatment of the outermost layers leads to irreducible errors in the calculation of the gravitational moments and thus on the inferred physical quantities for the planet. We have quantified these errors and evaluated the maximum precision that can be reached with the CMS method in the present and future exploitation of Juno's data.


Introduction
The concentric Maclaurin spheroid (CMS) method has been developed by Hubbard in Hubbard (2012) and Hubbard (2013) (H13). This method consists of a numerical hydrostatic scheme which decomposes a rotating celestial body into N spheroids of constant density (the well-known MacLaurin spheroid would correspond to the case N = 1). It needs two inputs. First, (i) the planet rotational distortion, q = ω 2 a 3 /GM , where ω denotes the angular velocity, a the equatorial radius, M the planet's mass and G is the gravitational constant. The factor q, which is linked to the MacLaurin's parameter, m = 3ω 2 /(4πρG), is also used in the theory of Figures (see for example Zharkov and Trubitsyn (1978) or Chabrier et al. (1992)), and represents the ratio of the rotational over gravitational potentials. Second, (ii) a barotrope (P ,ρ) representing the equation of state (eos) of the planet and its composition.
For a given density profile, we can calculate the gravitational potential of each spheroid (which is constant on the spheroid), and then calculate self-consistently the radius of each spheroid as a function of latitude. The obtained radii give different values for the potential of the spheroids, and successive iterations lead to convergence of both shape and potential of the spheroids (see Hubbard (2012) and H13). With the density and potential of each spheroid, one easily obtains the pressure from the hydrostatic equilibrium condition.
After this first stage of convergence, it is necessary to verify whether the obtained (P ,ρ) profile of the planet is in agreement with the prescribed barotrope. Since generally it is not, an outer loop is necessary to converge the density profile for the given eos.
At last, once the radii, potentials and densities of the multilayer spheroids have been obtained consistently with the required numerical precision, one obtains a discretized profile of these quantities throughout the whole planet. We can then calculate the gravitational moments, to be compared with Juno's observations. The formula is, by additivity: where J ext k is the moment of order k, J i k is the moment of order k of the i th spheroid, λ i is the ratio of the equatorial radius of the i th spheroid over the equatorial radius of the planet and the spheroids have been labeled with index i = 0, N − 1, with i = 0 corresponding to the outermost spheroid and N − 1 to the innermost one. Eq.(1) implies that the uncertainties of each spheroid are added up, so a small error on the radius, potential or density can lead ultimately to a significant error once summed over all the layers.
As explained in H13, for a given piecewise density profile, the evaluation of the gravitational moments is extremely precise, with an error around 10 −13 . The main sources of error, apart from the uncertainties on the barotrope, then arise from the finite number of spheroids and from the approximation of constant density within each spheroid, potentially leading to an incorrect evaluation of the radii and potential of the layers. In this paper, we evaluate the errors due to the discretization of the density distribution, and of the description of the outermost layers, the main aim of the paper is to find the best configuration in the CMS method to safely use it in the context of the Juno mission.
In Section 2, we show analytically that 512 spheroids are not enough to match Juno's error bars, and that the repartition of spheroids is an important parameter in the evaluation of the errors. We confirm these results by numerical calculations in Section 3. We then apply our method with a realistic eos in Section 4. These sections show the need to improve the basic CMS method. In Section 5, we study the impact of the outer boundary condition. Section 6 is devoted to the conclusion.
2 Analytical evaluation of the errors in the CMS method

Evaluation of the errors with 512 spheroids
In order to evaluate the errors analytically, we have considered a polytrope of index n = 1, P ∝ ρ 2 . We call ρ 2 the discretized version of ρ in the CMS model. The various quantities relevant for Jupiter are given in Table 1. If we assume longitudinal and north-south symmetry, only the even values of the gravitational moments are not null. Therefore, we were able to express the difference between the analytical value of the potential and the value obtained with the CMS method on a point exterior to the planet with radial and latitudinal coordinates (r, µ = cos(θ), where θ is the angle from the rotation axis), as: where P 2k is the Legendre polynomial of order 2k and a J is Jupiter's equatorial radius. Our first assumption here was that the CMS method prefectly captures the shapes of the spheroids. In reality, this is the case only for an infinite number of spheroids, but we have neglected this first source of error which is difficult to evaluate analytically. Since the masses of the two models must be the same, the k = 0 term must cancel. In terms of gravitational moments J 2k : By construction we can write ρ 2 (r) = ρ i , r i+1 < r ≤ r i . So the difference on the gravitational moment between the exact and the CMS models from Eq.(3) is: ri+1 (µ) [ρ(r ) − ρ i ] ×r 2k+2 P 2k (µ)dr dµ.
This formulation is exactly the same as in H13 (Eq.(10)), except that we used ρ i instead of δρ i . The analytical expression of the density of a rotating polytrope at first order in m as a function of the mean radius l of the equipotential layers and of the planet mean densityρ is given by (Zharkov and Trubitsyn (1978)): with A = π 3 1 − 2 3 m 1 − 6 π 2 and α = π 1 + 2m where l 0 is the outer mean radius.
At this stage, our method implied an important approximation, namely that ρ 2 followed perfectly the polytropic density profile of Eq.(5) on each layer, if evaluated at the middle of the layer: ρ i = ρ poly (r i + r i+1 ) /2. The errors we calculated are thus smaller than the real ones since we did not take into account those arising from the departure from the exact density profile. In Appendix A, we show how to get rid of the dependence on the angular part, µ, in Eq.(4). As detailed in Appendix B, for a linear spacing of the spheroid with depth, ∆r = a J /N , we obtained: with < P 2k >= 1 0 1 1 + e 2 µ 2 k+1 P 2k (µ)dµ . (8) It is interesting to note that the relative error of each layer is a decreasing function of i: the external layers lead to larger errors on the gravitational moment than the inner ones. Figure 1 displays the contribution to J 2 of the outer layers of the planet in the polytropic case (we have numerically integrated the expression for J 2 from Eq.(3), Eq.(5) and Eq.(33)).
To quantify the error, we focused on J 2 , the case k = 1. As detailed in Appendix C, one gets: With N = 512 (as used in H13), we get |∆J 2 | × 10 6 ∼ 3.6 × 10 −2 . The first data from Juno after three orbits give an uncertainty of 2.72 × 10 −1 on J 2 × 10 6 . We are under these error bars by a factor of eight, but in these analytical calculations, we have made several restrictive approximations. Notably, we have supposed that the numerical method can find perfectly the Table 1: Values of the planetary parameters of Jupiter. R eq and R polar are observed at 1 bar. The value of the pulsation is chosen following Archinal et al. (2011). (a) Cohen and Taylor (1987) Bolton et al. (2017) Parameter Value G a (global parameter) 6.672598 · 10 −11 ± 2 · 10 −17 m 3 kg −1 s −2  Figure 1: Contribution of each part of the planet to the total value of J 2 , normalized to its maximum value, for a polytrope of index n = 1.
shape and the density of each spheroid, and we have neglected the factor α 3 /3 ∼ 2 − 10 in Eq.(57). As mentioned earlier, it is hard to evaluate a priori the errors due to the wrong shapes and densities of the spheroids, but combined with the neglected factor, we expect the real numerical errors to become larger than Juno's error bars. This is confirmed by numerical calculations in the next section. Therefore, even in this ideal case, the CMS method in its basic form is intrinsically not precise enough to safely fulfill Juno's constraints. Calculations with 512 MacLaurin spheroids with a linear spacing cannot safely enter within Juno's error bars.
HM16 proposed a better spacing of the spheroids, with twice more layers above 50% of the radius than underneath. As shown in Appendix C.2, the new error for 512 spheroids becomes: This is about 2.5 times better than with the linear spacing. But once again, when considering the factor of ten neglected in Eq.(57) in nearly the entire planet and the strong approximation of perfect density and perfect shape of the discrete spheroids, it seems very unlikely to reach in reality a precision well within Juno's requirements. Again, this is confirmed numerically in §3.
We note in passing that the difference between the errors obtained with the above two different spacings shows the strong impact of the outermost layers in the method.
In conclusion, these analytical calculations show that, at least with 512 sheroids, the CMS method, first developed in H13 and improved in HM16, leads to a discretization error larger than Juno error bars. Increasing the number of spheroids would give a value of J 2 outside Juno's error bars and would thus require a change in the derived physical quantities (core mass, heavy element mass fraction ...). We investigate in more details the uncertainties on these quantities in §4.3.
It should be mentioned that Wisdom and Hubbard (2016) calculated similar errors by comparing their respective methods, finding (their Eq. (15)): while we get (12) The 0.7 constant term is about 20 times the value obtained in Eq.(10). This confirms that, indeed, we underestimated the errors in the analytical calculations, as mentioned above. These expressions show that one needs to use at least 2000 spheroids to reach Juno's precision, as can already be inferred from Wisdom and Hubbard (2016). These authors, however, suggest that 512 spheroids can still be used to derive interior models, because a slight change in the density of one spheroid could balance the discretization error. A change in one spheroid density, however, will affect the whole planet density structure, because of the hydrostatic condition. Changing the density of each spheroid (within the uncertainties of the barotrope itself) can thus yield quite substantial changes in the gravitational moments.
Even more importantly, if the method used to analyze Juno's data has uncertainties due to the number and repartition of spheroids larger than Juno's error bars, this raises the question of the need for a better precision on the gravitational data, since the derived models depend on the inputs of the method itself. It is thus essential to quantify precisely the uncertainties arising from the discretization errors, as we do in §4 and 5.
Another option would be to use Wisdom's method (Wisdom (1996)), which has virtually an unlimited precision. Table 3 of Wisdom and Hubbard (2016) indeed shows that its accuracy is greater than anything Juno will ever be able to measure. In its present form, however, it is not clear whether this method can easily handle density, composition or entropy discontinuities; it is certainly worth exploring this issue. The CMS method, in contrast, is particularly well adapted to such situations, characteristic of planetary interiors. As for the standard perturbative theory of Figures, it becomes prohibitively cumbersome when deriving high-order (fifth even order for J 10 ) moments (see e.g., Nettelmann (2017)). For these reasons, the CMS method appears to be presently the most attractive one to analyze Juno's data. As mentioned above, it is thus crucial to properly quantify its errors and examine under which conditions it can be safely used to derive reliable enough Jupiter models.

Spacing as a power of k
As shown in Appendix C.3, for a cubic distribution of spheroids as a function of depth, we obtained the error on J 2 : This is the same result as Eq.(9). We checked that we get the same results for a quadratic spacing or with a spacing of power of four (it is probably easy to prove that this holds for any integer power of n). This is interesting for several reasons. First, these types of spacings have more points in the high atmosphere region, where the neglected factor of Eq.(57) is closer to two than ten. The global neglected factor is then smaller than for a linear spacing. Second, we studied where the cubic spacing is equal to the linear one. Using Eq.(72), we got: More generally, with a spacing of power k, ∆r k i is smaller than the linear ∆r for i ≤ N/k 1 k−1 , and larger for larger values of i. That means that above r i = a J 1 − 1/k 1 k−1 the radius, potential and density of the spheroids are better estimated than in the linear case and less well below this radius. And the stronger the exponent in the spheroid distribution, the larger the external gain (internal loss) on the spheroids because the more (less) tight the external (internal) layers.
With a quadratic spacing, the spheroids are closer than in the linear case above 75% of the radius of the planet while inside this limit the spacing is worse than in the linear case. With a cubic spacing the precision in the spheroids repartition is better than in the linear case down to only 80% of the radius, while for a power of 4 the limit is 85%. As shown and explained in §3.1, the cubic spacing ends up being the best compromise.
There is one remaining problem with a cubic spacing: for N > 512, the sizes of the first (outermost) layers are extremely small, well below the kilometer, with densities smaller than 10 −3 kg m −3 . It is thus mandatory to have an eos accurate enough in this regime; if not the errors will be very large. The other solution, that is discussed in §5 is to impose a non zero density as outer boundary condition.

Exponential spacing of the layers
Here, we examine an exponential repartition of spheroids. The details of the calculations are given in Appendix C.4, and the repartition functions of the spheroids are displayed in Figure 2. For the exponential repartition, the value of the error is: where γ ∈ [5, 10] is a parameter which determines how sharp we want the exponential function to be. A large value of γ leads to a steep slope in the innermost layers and a very small ∆r i in the exterior. With N = 512: |∆J exp 2 | × 10 6 ∈ [1.54 × 10 −2 , 1.6 × 10 1 ].
As seen, the smaller γ the smaller the error, but also the larger the range of outer layers which is a problem, as examined in §3. Note also that the neglected factor of Eq.(57) in the analytical calculations is much smaller because of the very large number of spheroids in the outermost part of the envelope (see Figure  2). Compared to Eq.(10), we see that the exponential error is always larger for γ > 3. As discussed in the next section, however, the impact of the first layers becomes very important with this type of spheroid repartition.
In summary, in this section, we have shown analytically that 512 layers spaced linearly is not enough to obtain sufficient accuracy on the gravitational moments of Jupiter to safely exploit Juno's data. Indeed, the theoretical errors are within Juno's error bars but the simplifying assumptions used in the calculations (in particular perfect shape and density of the spheroids) suggest that the real calculations will not match the desired accuracy. As explained in §1, the method, mathematically speaking, is precise up to 10 −13 for a piecewise density profile, but the errors due to the discretization, as derived in this section, are of the order of ∼ 10 −7 . Using a better repartition of spheroids yields, at best, an error of the order of Juno's error bars but, in any case, remains insufficient with 512 spheroids. The impact of such an error on Jupiter's derived physical quantities is discussed in §4.3. According to this analytical study, using a minimum of 1000 spheroids spaced cubically, quadratically or exponentially seems to provide satisfying solutions. That may sound surprising since the spacing used in HM16, with their particular first layer density profile, gives a smaller uncertainty in Eq.(10) than the cubic or square one in Eq.(14). To explain this apparent contradiction, we need to explore more precisely how the error depends on the different spacings by comparing numerical calculations with the analytical predictions. This is done in the next section. In this section, we explore two main issues: how does the error depend (i) on the number of layers and (ii) on the spacing of the spheroids. Another parameter we did not mention in the previous section is the repartition of the first layers. This is the subject of the next section.
To study this numerically, we have developed a code similar to the one developed by Hubbard, as described in H13. The accuracy of our code has been assessed by comparing our results with all the ones published in the above papers. In the case of constant density, we recovered a MacLaurin spheroid up to the numerical precision, while for a polytrope of index n = 1, we recovered the results of Table 5 of H13, as shown in Table 2 (the differences are due to the fact that the published results were not converged to the machine precision. With the very same conditions, our code and Hubbard's agree to 10 −14 ). We have also performed several tests of numerical convergence with different repartitions of spheroids. We have compared our various numerical evaluations of the errors, double precision with 30 orders of gravitational moments and 48 points of Gauss-Legendre quadrature, to a quadruple precision method with 60 orders of moments and 70 quadrature points, from a couple of hundred spheroids to 2000 spheroids. The differences are of the order of 10 −13 , negligible compared with Juno's error bars. Finally, as explained in Hubbard et al. (2014), we have implemented an auto exit of the program if the potential diverges from the audit points.  Figure 3 displays the difference between the expected and numerical values of J 2 for a polytrope of index n = 1, for different spacings with 512 layers. The first thing to note is that, whatever the repartition, the numerical errors are always larger than Juno's error bars: 512 spheroids are definitely not enough to match the measurements of the Juno mission with sufficient accuracy, as anticipated from the analytical calculations of the previous section.
The second and quite important problem, as briefly aluded to at the end of §2, is the difference between calculations using a density ρ = 0, as in HM16, or ρ 0 in the first layer. With 512 spheroids spaced as in HM16, the first layer is 50 km deep. Imposing a zero density in this layer is equivalent to decreasing the radius of the planet by 50 km, which is larger than the uncertainty on Jupiter's radius. As seen in the Figure, if one does so, the error on J 2 is more than 100× Juno's error bar. Specifically, for the HM16 spacing, when the first layer has a non-zero density, its impact on the value of J 2 × 10 6 is ∼ 5 × 10 1 . We discuss that in more details in §3.2. Figures 4 and 5 display the behavior of the error as a function of the number N of spheroids for a cubic, square and a few exponential repartitions, respectively. As mentioned earlier, we see in Figure 4a that only for N > 1000 spheroids do we get an error smaller than Juno's error bar. Importantly enough, we also see that the first spheroid (N = 0) makes no difference on the result when the layers are cubically spaced (Figure 4b). : Absolute value of the difference between the expected and numerical results on J 2 × 10 6 for different spacings of the spheroids with N=512. Lin, squ and cub are repartitions of spheroids following a power law of exponent 1, 2 or 3, respectively; exp γ=8 is the exponential repartition with our favored value γ = 8. HM is the method exposed in HM16, explained in Appendix C.2. The subscript 1 means that the first layer has a non-0 density while the subscript 0 means zero density in the first layer.
This contrasts with the results obtained with the HM16 spheroid repartition, where we see that the first spheroid has a huge impact on the final result ( Figure 3). The cubic spacing with at least 1000 spheroids thus matches the constraints, with both an uncertainty within Juno's error bars and a negligible dependence on the first layer. As mentioned in the previous section (and shown later in Figure 6b), however, with the cubic spacing the densities of the first layers are orders of magnitude smaller than 10 −3 kg m −3 , questioning the validity of the description of the gas properties by an equilibrium equation of state. For this reason, it is preferable to use a square or exponential spacing in the case of a zero pressure outer boundary condition. The effect of changing this boundary condition is discussed in §5.
As seen in Figure 4a, the square spacing requires at least 1500 spheroids to be under the error bars of Juno with reasonable accuracy. For 1500 or 2000 spheroids, the first layers have densities around 10 −2 kg m −3 . In Figure 4b, we confirm that, in that case, having a zero or non-zero density outer layer changes the value under the precision of Juno.
For the exponential spacing, we recall that the higher the γ value the smaller the first layer, which can become a problem (see §3.2). Considering Figures 5a and 5b, the best choices for γ are γ = 8 or γ = 7, because it does not require too many spheroids to converge and at the same time the obtained values for the densities of the first layers are comparable to the ones obtained with the square repartition (slightly larger for γ = 7).
Globally, we conclude in this section that we need at least 1500 and preferentially 2000 spheroids spaced quadratically or exponentially to fulfill our goals of both high enough precision for Juno's data and densities high enough to justify the use of an equilibrium eos. Having two different repartitions that match Juno error bars brings some confidence in results that are obtained with a realistic eos instead of a polytrope. Furthermore, with such choices of repartition, we do not need to impose a first layer with zero density. In §5, however, we show that the cubic repartition can also be used, with a different boundary condition.

Importance of the first layers
As explained in HM16, these authors improved the spacing of their spheroids by adding up a first spheroid at depth a J −∆r/2, instead of ∆r, which has zero density: ρ 0 = 0. Although this could seem quite arbitrary, it has an analytical justification emerging from the linear dependence of the density with radius (Hubbard, priv. com.). The obtained result is indeed closer to the analytical value of Hubbard (1975), but it has two downsides. First, it decreases arbitrarily the radius by 50 km (see discussion above). For example, for a cubic spacing with 512 spheroids, decreasing the radius by 50 km leads to J 2 × 10 6 = 13961.23 with an error of 2.7×10 1 according to the analytical calculation, 100× Juno's uncertainty. Furthermore, there is no guarantee that this repartition will still be valid with a realistic eos, with an exponential dependence of the density upon the radius. In order to assess the robustness of the results and of the CMS method, one needs to confirm that different repartitions yield similar values and that the results are not only correct for one given spheroid repartition.
More generally, a key question is: why are the first layers so important in the CMS method ? Going back to Figure 1, we see that the region outside 90% in radius contributes as much to J 2 as the region within the inner 50%. The problem is the slope of the moment: it is much steeper in the outside region. Consequently, the behavior of a layer in the outermost region yields a much larger error than the contribution from the inner layers. This simply reflects the fact that most of the planet angular momentum is in its outermost part. Another issue is that, in the first layers, the change in density is represents at least 10% of the density itself and can even become comparable to this latter. In contrast, in the internal layers, even a large ∆r yields a relative change of density from one spheroid to another less than 1%. Then, the error due to the inevitable mistake on the true shape of the spheroids is much more consequential in the external layers, where the relative density change is significant. Another source of error stems from the fact that the pressure is calculated iteratively from the hydrostatic equation, starting from the outermost layers. A small error in the first layers will then propagate and get amplified along the density profile.
Equation (10) of HM13 gives the explicit expression of J 2 for each layer. With the notation of the paper: To illustrate our argument, we said that every spheroid has the Difference from analytical value Juno s error bar Juno s error bar (b) Figure 5: Numerical error as a function of the number of spheroids for exponential repartitions. (a): error between the analytical and measured J 2 × 10 6 ; (b): difference on the J 2 × 10 6 when the first layer has 0 or non-0 density. same shape and that δρ i is the same as in the rest case, given by a slight modification of Eq.(63) of HM13: Then, the only different terms for each layer stemmed from the denominator: δρ i × λ 5 i . We plot this value for different spacings of the spheroids in Figure 6.
As seen, the linear spacing is much less precise than the other ones in the external layers because of the very small number of layers; it will thus greatly enhance the numerical errors. The cubic spacing is much better but has so many points in the outermost layers that it is oversampling the extremely high part of the atmosphere. The exponential spacing is the worst one in the lower part of the envelope (see Figure 6), but provides a nice cut-off between linear and cubic spacings in the uppermost layer An appealing solution would be to take a combination of all these spacings to always have the most optimized one in each region of the planet. We have tried this, with a great number of spheroids, and the conclusion remains the same: the limiting effect on the total errors is the precision on the external layers (above ∼ 95% of the radius). Having simply an exponential repartition of spheroids throughout the whole planet or the same exponential in the high atmosphere and then a cubic repartition when this latter becomes more precise, and then square and linear repartitions (such a repartition almost doubles the number of spheroids) yield eventually the same error. To obtain results within Juno's precision (see Figures 4 and 5), we need the external layers to be smaller than 1 km.
To conclude this section, we want to stress that the main limiting factor in the CMS method is the number of spheroids in the outmermost layers (above ∼ 0.95 × a J ) and the proper evaluation of their densities. If we do not want to impose a first layer with 0 density, with consequences not predictable in general cases, we must be able to describe this region of the planet with sufficiently high accuracy. This requires an optimized trade-off between the number of spheroids in the external layers and the value of the smallest densities. As shown in §4.3, the treatment of these layers in the CMS method is the largest source of error in the derived physical quantities of the planet.
One could argue that if the outermost layers can not be adequately described by an eos, a direct measurement of their densities would solve this issue. Unfortunately, beside the fact that, by construction, Mac Laurin spheroids imply constant density layers, there is an additionnal theoretical challenge: in the high atmosphere, the hydrostatic balance barely holds. The isopotential surfaces are not necessarily in hydrostatic equilibrium, as differential flows are dominating. The time variability in these regions is an other problematic issue. At any rate, the treatment of the high atmospheric layers is of crucial importance in the CMS method, because it determines the structure of all the other spheroids. An error in the description of these layers will propagate and accumulate when summing up the spheroids inward, yielding eventually to large errors. In Section 5, we examine the solution used in Wahl et al. (2017) to overcome this issue.

1 bar radius and external radius
So far, Jupiter outer radius in the calculations has been defined as the observed equatorial radius at 1 bar, R 1bar = 71492 km (see Table 1). A question then arises: does the atmosphere above 1 bar have any influence on the value of the gravitational moments ?
When converging numerically the radius at 1 bar to the observed value, we obtained a value for the outer radius R ext = 71505 km. Integrating Eq.(3) from R 1bar to R ext by considering a constant density equal to the 1 bar density, with the simplifications of Appendix A, yields a contribution of the high atmosphere J 2 <1bar × 10 6 1 × 10 −2 .
It thus seems reasonable to use the 1 bar radius as the outer radius with the current Juno's error bars. In §5, however, we examine in more detail this issue and quantify the impact of neglecting the high atmosphere layers in the calculations, as done in Wahl et al. (2017). In the next section, we now examine if the conclusions of §3 still hold when using a realistic eos.

Calculations with a realistic equation of state
The aim of this paper is not to derive the most accurate Jupiter models or to use the most accurate eos, but to verify (i) under which conditions is the CMS method appropriate in Juno's context, (ii) which spheroid repartition, among the ones tested in the previous sections, is reliable when using a realistic eos to define the barotropes, and (iii) whether the high atmosphere is correctly sampled to induce negligible errors on the moment calculations. For that purpose, we have carried out calculations with the Saumon, Chabrier, Van Horn eos (1995, SCvH). This eos perfectly recovers the atomic-molecular perfect gas eos, which is characteristic of the outermost layers of the planet (Saumon et al. (1995)).

Impact of the high atmospheric layers on the CMS method
Concerning the last of the three aforementioned questions, the SCvH eos yields much better results than the polytrope. The polytropic densities in the high atmosphere are much larger than the realistic ones, yielding an overestimation of the contribution of the external layers. Figure 7 illustrates the contribution to J 2 of the outermost spheroids with the SCvH eos. Even though, as mentioned earlier, the very concept of an equilibrium equation of state becomes questionable at such low densities (P < 100 Pa = 10 −3 bar corresponds to ρ < 10 −3 kg m −3 ), this clearly illustrates our purpose. : Contribution of each spheroid to the total value of J 2 × 10 6 and to the truncated J 2 × 10 6 , i.e., the sum of the contributions of the layers from 1 to 10 −5 bar, for a square spacing with the SCvH eos with 1500 spheroids. Figure 7 shows that, when using a realistic perfect gas eos, using only 1 spheroid up to 10 −3 bar seems to be adequate to fulfill Juno's constraints. Unfortunately, the real problem is more complex: if we use only one spheroid of 8 km (a few 10 −3 bar) as a first spheroid, even though its own influence is indeed negligible, as if we divide it into many spheroids, the total value of J 2 ×10 6 is changed by ∼ 1−2, ten times Juno's error ! The reason is that the whole density profile within the planet is modified by the inaccurate evaluation of the first layers. Indeed, remember that the CMS method decomposes the planet into spheroids of density δρ. Therefore, a wrong evaluation of the first layers is propagated everywhere since the density at one level is the sum of the densities of all spheroids above that level. We verified that any change in the size of the first layer, from less than 1 km to 35 km significantly changes the value of J 2 . We have also tried a solution where we converged the 1 bar radius to the observed R 1bar value, in order to avoid this modification in the structure, but the result still holds: it is not possi-ble to reduce the global contribution of the first layers, and the densities of the first layers are too small to be realistically described by an equilibrium equation of state. Clearly, with the CMS method, it is not possible to diminish the size of the first layers without significantly affecting the gravitational moments. This, again, is examined in detail in §5. Table 3 gives the results for an exponential repartition of 1500 spheroids using the SCvH eos, with no core (in that case, the size of the first layer is 0.6 km). Clearly, varying the size of the first layers yields changes of the J 2 value larger than Juno error bars. A quadratic repartition yields similar results with differences on J 2 × 10 6 that are < 3 × 10 −1 . It is not clear, however, whether the induced errors stem from the iterative calculations of the pressure from the hydrostatic equation or from the error on the potential determination. Therefore, we confirm with a realistic eos the results obtained in Section 3.2: a very accurate evaluation of the first layers is mandatory in order to properly converge the CMS method.

Errors arising from a (quasi) linear repartition
In order to assess quantitatively the errors arising from an illadapted repartition of spheroids, we have carried out 3 test calculations. One without a core, where we changed the H/He ratio (i.e., the helium mass fraction Y ) to converge the CMS method (as explained in HM16, that means having a factor β = 1 in H13). A second one where we imposed Y and the mass of the core but allowed the size of the core to vary (we have verified that imposing the mass or the size leads to the same results).
And, finally, one where Y and the mass of the core could change and we converged J 2 to the observed value. These models were not intended to be representative of the real Jupiter but to test the spacings of the CMS explored in the previous sections. All these tests used the observed radius as the outer radius, R eq = 71492 km. The results are given in Table 4. The difference on J 2 × 10 6 between the square and exponential spacings with 1500 spheroids are < 2 × 10 −1 , as expected from Figures 4a and 5a. In contrast, the HM16 method, which is linear in ∆r, though with a different spacing in the high and low atmosphere, respectively, leads to differences 100× larger. These values were obtained with a HM16 spacing with 1500 spheroids; with 512 spheroids the differences from square and exponential are a factor ten larger. This demonstrates that the square and exponential repartitions remain consistent, confirming the results obtained for the polytropic case, whereas the HM16 spacing diverges completely. As expected, the idea of using a first half layer of zero density fails when using a realistic eos.
Another proof is the difference, with the HM16 repartition, between 512 and 1500 spheroids: in the polytropic case the difference on J 2 × 10 6 was found to be around 2 × 10 −1 . Indeed the HM16 repartition with 512 spheroids led to results very close to the theoretical value. In our first test case, the value of J 2 × 10 6 we obtained for the HM16 spacing with 512 spheroids is 15075.52. When we compare to Table 4, the difference is 5 × 10 1 : the discretization error on the gravitational moments with this spacing is way off Juno's precision. Moreover, we see that the order of magnitude of the error is similar to the one in the linear case of Figure 3. Indeed, if one does not use a polytrope, the HM16 spacing behaves as badly as a linear spacing, as expected from the quasi linear shape of this repartition. Figure 8 summarizes graphically the uncertainties arising from using an ill adapted repartition of spheroids. As stated above, the errors on J 2 from the bad handling of the high atmosphere reach almost 100× Juno's error bars, as shown in Figure 8a. On the other hand, when fitting J 2 , Figure 8b shows that the errors on J 4 and J 6 are comparable or inferior to Juno's. Nonetheless, we show in §4.3 that this leads to significant changes in the physical quantities, hampering the determination of interior models in the context of Juno. These changes can be considered as negligible compared to the uncertainties on the various physical quantities in the models but they are significant when one aims at matching Juno's measurements and thus take full advantage of the Juno mission.
What we have shown so far in the present study is that the CMS method, when taking into account the entire atmosphere above 1 bar, relies on uncontrolled assumptions that can change significantly the value of J 2 . The whole purpose of the method, however, is to constrain with high precision Jupiter's internal structure (mass fraction of heavy elements, mass of the core ...). Therefore, we need to assess the impact of these assumptions on such quantities.

Intrinsic uncertainties on the Jupiter models
In this section, we derived interior models that match both J 2 and J 4 from the first Juno data (Bolton et al. (2017)) and determined how they were affected by the number of spheroids and  Figure 8: Comparison of the errors of the CMS method from quadratic and HM16 repartitions with Juno's error bars, shifted to be centered on the observed J n . (a) −J 4 × 10 6 vs J 2 × 10 6 from test 1; (b) J 6 × 10 6 vs −J 4 × 10 6 from test 3, when J 2 is fit.
their repartition. Ideally, one would like the final models not to depend on the method used to derive them. We recall that the purpose of this paper is not to derive the most accurate Jupiter models but to determine the impact of the uncertainties in the CMS method upon these models, in the context of the Juno mission. Therefore, we used the SCvH eos with an effective helium mass fraction (Y ) which includes the metal contribution (see Chabrier et al. (1992)). We supposed that the core was a spheroid of constant density ρ core = 20000 kg m −3 , the typical density of silicates at the pressure of the center of Jupiter. At 50 GPa, we imposed a change in the effective Y at constant temperature and pressure, thus a change in entropy, to mimic either an abrupt metallization of hydrogen or a H/He phase separation. Inward of this pressure, the value of the effective helium mass fraction is denominated Y 2 . All the models were converged to the appropriate J 2 and J 4 values with a precision of 10 −8 , within Juno's error bars. Table 5 gives the values of Y , Y 2 and the core mass for different repartitions and numbers of spheroids. We considered a quadratic repartition with an unchanged first layer or with a first layer of 10 km, and the HM16 repartition. With HM16, the change on the effective helium mass fraction is of the order of 5%, consistent with the expectations from the results of the previous sections. This confirms that the HM16 repartition of spheroids is not suitable to exploit Juno's data. When we used the Y and core mass fraction obtained with 512 spheroids to a case with 1500 spheroids, the value of J 2 × 10 6 became 14647.9, which leads to an error of 5 × 10 1 , a hundred times Juno error bars. As stated in the abstract, this result demonstrates that Jupiter models which have been derived with this spheroid repartition are invalid, in the framework of Juno's constraints. The most recent calculations of Wahl et al. (2017) are examined specifically in §5. Table 5: Effective helium mass fraction above (Y ) and underneath (Y 2 ) the entropy jump, and mass of the core (M c in units of 10 24 kg) for different repartitions and numbers of spheroids. Quadratic is a power law spacing of exponent two, Quadratic 10 is the same with a first layer of 10 km and HM16 corresponds to the repartition used in Hubbard and Militzer (2016)  With the quadratic repartition, using the values obtained with 512 spheroids in the 1500 spheroid calculations yields an error on J 2 × 10 6 of 4 × 10 0 , ten times Juno's error bars (which confirms the necessity to use at least 1500 spheroids and this type of repartition). When correctly fitting J 2 and J 4 with this repartition, the change on Y from 512 to 1500 spheroids is about 0.5 %, about the same as the difference between using a < 1 km or a 10 km first layer. For the core mass, increasing the number of spheroids is useless as the difference is negligible compared to the one due to changing the size of the first layer. Therefore, the uncertainties due to the description of the first layers are the dominant sources of limitation in the determination of Jupiter's physical quantities like Y or the core mass fraction.
In reality, the limitation is even more drastic. First, with a 10 km first layer, the outer density is 3 × 10 −3 kg m −3 . This is too low to use an equilibrium equation of state. If we use instead a 20 km first layer, increasing the density to more acceptable values, the uncertainties on Y are about 1-2% and > 5% on M c . Virtually any model leading to such an uncertainty can fit J 2 and J 4 . Indeed, given the impact of the first layers on the gravitational moments, they could be used, instead of Y and M c , as the variables to fit models with the appropriate gravitational moments. Slightly changing these layers (mean density, size,..) can yield adequate models. Therefore, no model can be derived within better than a few percent uncertainties on Y , thus the heavy element mass fraction, and on the core mass.
The calculations carried out in this section yield three major conclusions: 1. It is not possible to use a linear repartition, especially with a zero density first layer, as the values of the moments would be changed by 100× the value of Juno's error bars when increasing the number of spheroids.
2. The results obtained with a square or an exponential repartition with 1500 spheroids yield consistent errors, no matter the size of the first layer, and can thus be considered as reliable spheroid repartitions. The size of the first layer, however, remains the major source of uncertainty in the CMS method when taking into account the high atmosphere above the 1 bar level.
3. When deriving appropriate internal models of Jupiter, the CMS method such as used in HM16 can not allow the determination of key physical quantities such as the amount of heavy elements or the size of the core to better than a few percents. Changing the size and/or density distribution of the first layers yields a significant change in the values of the gravitational moments. This degeneracy unfortunately hampers the derivation of very precise models.
In the light of these conclusions, the CMS method, such as used so far in all calculations, needs to be improved. To do so, Wahl et al. (2017) propose to neglect the high atmosphere, and define the outer boundary at the 1 bar level. We examine this solution in the next section.
calculations. As shown with a polytropic eos in §3.2, the atmosphere above this level, called here the high atmosphere, has a contribution to J 2 × 10 6 1 × 10 −2 . With a realistic eos, converging the 1 bar radius on the observed value yields a high atmosphere depth of ∼ 70 km but a smaller outer density than the polytrope. The contribution to J 2 is then about the same.
The high atmosphere, however, also alters the evaluation of the potential inside the planet and thus the shape of the spheroids in the CMS method. In Appendix D, we evaluated the first order neglected potential φ neg on each spheroid due to the omission of the high atmosphere. We showed that it is a constant value, which does not depend on depth or angle: where ρ ext is the (constant) density of the high atmosphere,ρ S is the mean density of Jupiter if it was a sphere of radius a 1bar = 71492 km, and ∆a is the depth of the high atmosphere. As mentioned above, ∆a 70 km and we know that ρ S 10 3 kg m −3 . Taking the outer density in a range (0.01 − 1) × ρ 1bar , and the conditions of Jupiter, ρ 1bar 0.17 kg m −3 , yields: This is clearly a first order correction since the potential on each spheroid, given by the CMS method, ranges from GM/a 1bar in the outermost layers to 2GM/a 1bar inside. The neglected potential does not depend on the layer, so the hydrostatic assumption (using the gradient of the potential) remains valid. The main source of change comes from the shapes of the spheroids. We calculated the impact of this change in Appendix D with the use of equation (40) of H13 and obtained the following result : Within the above range of outer densities we get: Since the potential is smallest in the outside layers, we expect the relative change in J 2,i of the i th spheroid to decrease with depth. Nevertheless, we have seen that the CMS method requires a large number of spheroids in the high atmosphere (a point which is validated in the next section). We thus expect the change in the J 2,i in the outer layers, which have the largest impact on J 2 , to be comparable to Eq.(26). If we make the bold assumption that every layer contributes equally and choose for the outer density ρ ext = 0.03 kg m −3 , which is the value that corresponds to mass conservation in the high atmosphere, we get : ∆J ext 2 × 10 6 8 × 10 −6 (J ext 2 × 10 6 ) 1 × 10 −1 . (27) This is approximately 1/2 Juno error bars. Even though this relies on several assumptions, mainly that we can approximate the outer layers by a constant density spheroid and that the error is the same on each spheroid, it gives a reasonable estimate of the impact of neglecting the high atmosphere on J 2 . From the evaluation of J 4 and J 6 in Appendix D, we calculated that: ∆J ext 4 × 10 6 6 × 10 −2 , ∆J ext 6 × 10 6 2 × 10 −2 .
We see that the agreement with Juno's error bars improves with higher order moments, although it stabilizes around 10 −2 from J 6 onward. In order to derive physical information from the high order moments, however, one needs to have confidence in the evaluation of the first moments. Moreover, these estimations rely on many assumptions and give only orders of magnitude. As in Figure 8, we show graphically these uncertainties in Figure 9. We see that we are now able to reach Juno's precision. As mentioned in the abstract, however, the analysis of future, more accurate data will remain limited by these intrinsic uncertainties.
If we were able to perform rigorously the calculations of Appendix D, we would probably find that the errors on the J n 's are correlated, decreasing the aforederived uncertainties. Unfortunately, a rigorous derivation of these errors would require to relax many assumptions and would imply very complicated theoretical and numerical calculations. Until this is not done, the determination of the J n are inevitably blurred by the neglected potential, up to half Juno's error bars, as shown in Figure 9.
In conclusion, we have shown in this section that the neglected potential due to the high (<1 bar) atmosphere region leads to an error of about half Juno error bars. Globally, neglecting the atmosphere above 1 bar in the CMS method leads to an irreducible error on J 2 × 10 6 of the order of 10 −1 . Although such an error level is smaller than Juno's present precision, this will become an issue if the precision improves further.
It must be stressed that this irreducible error is not arising from uncertainties in any physical quantity, such as for example the radius determination or the eos. It stems from neglecting the high atmospheric levels in the CMS method itself, and in that regard cannot be removed; indeed, as we have shown in the previous sections, including the high atmosphere leads to huge uncertainties. Said differently, even with a perfect knowledge of Jupiter radius, internal eos, etc., the error, intrinsically linked to the method, remains. Semi analytical integrations of Eq.(108) shows that, globally, this error is of the order of 10 −2 for higher order moments. Therefore, the J k × 10 6 cannot be determined within better than a few 10 −2 . On J 10 and higher order moments, in particular, this error can be close to the order of magnitude of the moments themselves. One must keep that in mind when calculating values of these moments with the CMS method, as in Wahl et al. (2017).

Error from the finite number of spheroids
Now that we have estimated every possible source of error and showed that the method used by Wahl et al. (2017) is applicable with the current error bars of the Juno mission, we need to examine how the uncertainties depend on the number of spheroids and verify whether the conclusions of §3 are affected or not by the change of outer boundary condition. Figure 10  dependence of the error upon the number of spheroids for different spheroid repartitions in the present case. Interestingly, the exponential γ = 6 repartition happens to provide the best match to the analytical results. More importantly, we show that some of the different repartitions give the same results, a mandatory condition to use the method with confidence with a realistic eos, as mentioned earlier.
We see, again, that 512 spheroids are not enough to satisfy Juno's precision. As mentioned previously, even though the CMS values are precise at 10 −13 , the discretization error dominates. Depending on the spheroid repartition, this error can become quite significant. Unfortunately, Wahl et al. (2017) do not mention which spheroid repartition they use exactly. As shown in the Figure,  lidity. When neglecting the high atmosphere, the calculations are not hampered by any degeneracy due to the first layers since the densities, which are larger than ρ 1bar 0.17 kg m −3 are well described by the eos. Therefore, the only errors on the derived quantities of interest come from the discretization errors. As in §4.3, we give in Table 6 the values of Y , Y 2 and M c obtained when matching J 2 and J 4 in the same conditions. One might notice that these results are different from the ones in Table 5; this is only due to the fact that in Table 5 the outer radius is 71492 km, instead of the radius at 1 bar. We checked that if the external radii are consistent, the derived physical quantities are within the range of uncertainties of §4. The differences from 512 to 1500 spheroids are of the order of 0.3% on the physical quantities. Because of the degeneracy of models leading to the same external gravitational moments (by e.g., changing the pressure of the entropy jump, the amount of metals, etc..) we do not expect to derive a unique model with such a precision on the core mass, so, in this context, 512 spheroids can be considered as sufficient to derive accurate enough models of Jupiter. When using the physical quantities obtained with 512 spheroids in calculations done with 1500 spheroids, we get an error on J 2 × 10 6 of ∼ 1, 10× Juno error bars. Therefore, even though one can derive a precise enough model with 512 spheroids, it is highly recommended, as a sanity check, to verify with 1500 whether Juno precision can be reached with only a slight change of the derived physical quantities.
We still have to evaluate the errors coming from neglecting the high atmosphere. To evaluate the largest possible uncertainty, we added 1 × 10 −1 to J 2 × 10 6 and removed it from J 4 × 10 6 . We stress again that this must be done because the value of both J 2 and J 4 can not be trusted under the aforederived intrinsic irreducible error of the method, due to the omission of the high atmosphere, even when fitting to observed values at much higher precision. The results are given in Table 7. The differences between Table 6 and 7 show that neglecting the high atmosphere yields an uncertainty on the physical quantities of a few parts in a thousand, even up to 2 to 3% on the core mass.  Table 6 except that J 2 ×10 6 and J 4 ×10 6 have been changed by ±1 × 10 −1 . +1-1 means that we add 1 × 10 −1 on J 2 and remove it from J 4 , -1+1 is the opposite. This is the irreducible error of the CMS method when neglecting the levels of the high atmosphere (<1 bar). An improvement in the precision of Juno cannot change this uncertainty. In that regard, using more than 1500 spheroids would be a second order correction on the method uncertainty.

Conclusion
In this paper, we have examined under which configuration the CMS method can be safely used in the context of the Juno mission. We have first derived analytical expressions, based on a polytropic eos, to evaluate the errors in the calculation of the gravitational moments with the CMS method in various cases. These analytical derivations relied on simplifying approximations and thus represent lower limits of the errors. When applied to calculations based on 512 spheroids, as done for instance in the works of Hubbard and Militzer (2016) these expressions show that it is not possible to match the Juno's constraints because of the discretization error. To verify the analytical calculations, we have developed a code based on the method exposed by Hubbard (2013). The numerical calculations confirmed the analytical results and also showed that the external layers of the planet have a huge impact of the determination of J 2 and need to be described with very high accuracy. At first, this appeared to be the biggest limitation of the CMS method: for a precise evaluation of the gravitational moments of Jupiter, this method requires layers in the high atmosphere smaller than 1 km. However, this implies densities so low that the very concept of an equilibrium equation of state to describe these layers becomes questionable, while density profiles fit from observations would not solve this issue (as explained at the end of §3.2).
We have shown, with a polytrope n = 1, that one needs a quadratic or exponential repartition of at least 1500 spheroids to reach precise enough values of the moments to match Juno's constraints, which implies densities of the order of a few 10 −2 kg m −3 for the outermost layers. With such specifications, the CMS method can in principle reach precisions better than Juno's uncertainties.
When these calculations are performed numerically with barotropes prescribed from an appropriate eos for the external layers, the above conclusions were confirmed: the CMS method requires a square, cubic or exponential spacing of at least 1500 spheroids. None of these calculations, however, was able to get a satisfying description of the first layers. Indeed, even though each of these individual layers contributes negligibly to J 2 , a small change in the density structure of these regions leads eventually to unacceptable changes in the total value of J 2 . Furthermore, when fitting J 2 and J 4 , we have shown that the description of these layers leads to a degeneracy of solutions which hampers the determination of physical quantities such as the heavy element mass fraction or the mass of the core to better than a few percents.
In order to overcome this issue, Wahl et al. (2017) propose to neglect the high atmosphere and impose the 1 bar pressure as an outer boundary condition. We have evaluated the errors coming from this assumption, and showed that they are of the order of half Juno error bars. In principle, the CMS method with this boundary condition can thus be applied to derive Jupiter models within the current precision of Juno, with an irreducible uncertainty arising from the neglected high atmosphere. This uncertainty is of the order of 0.1 on J 2 × 10 6 while the current precision of Juno is ±0.272. The remaining sources of errors on the physical quantities, however, stem from discretization errors, and 512 spheroids can give results with a precision of a few parts of a thousand. Considering the degeneracy of models that could fit with Jupiter's observed gravitational moments, this can be considered as an acceptable uncertainty. However, we show that, due to the aforementioned irreducible error of the method, increasing Juno precision will not enable us to derive more precise models. Quantitatively, the irreducible errors yield global uncertainties on the physical quantities derived in a model close to 1%.
The conclusion of this study is that the Concentric MacLaurin Spheroid method such as used in HM16 could not satisfy Juno's constraints and needed to be improved to derive accurate enough Jupiter model, invalidating any model up to this stage in the context of the Juno mission. First, a larger number of spheroids with a better repartition is mandatory. Second, there is no satisfying solution to safely take into account the impact of the high atmosphere region above 1 bar. We have shown that imposing the 1 bar radius as the outer boundary condition, as done in Wahl et al. (2017), is acceptable, although it leads to irreducible errors of a few parts of a thousand on the derived physical quantities, increased by an eventual discretization error.
Although such a precision can be considered as satisfactory, given all the other sources of uncertainty in the input physics of the model, this shows that, if the precision on Juno's data improves, the CMS method in its present form will not be able to exploit this improvement to refine further the models. Zharkov, V. N. and Trubitsyn, V. P. (1978). The physics of planetary interiors.

A Simplifications for the angular part
Here, we show in details how we got rid of the angular part in Eq.(4). The integral we wanted to simplify is: We approached this integral differently: instead of just considering an integral over r and µ, we rather followed the surfaces of constant density. Then, r (µ) was defined in order to have ρ(r (µ)) = cst when µ varies from 0 to 1. Another way to say it is that we took the value of ρ for a given equatorial radius r eq , and calculated the integral over constant ρ from equator to pole. In order to go further, we considered that the surfaces of constant potential were ellipsoids with the same ratio of polar to equatorial radius e. That gave, for every equipotential surface: where r eq is the equatorial radius of the spheroid. In reality, this is not exactly the case, the isobaric shperoids are not ellipsoids and the innermost ones are less flattened than the outermost ones, but it gives an idea of the integral over µ.
We just had to change variables : µ → µ and r → r eq / 1 + e 2 µ 2 . With µ varying from 0 to 1 and r eq from 0 to a J , we have not changed the domain of integration because outside of the outermost ellipsoid -in our ellipsoidal approximation -we had ρ(a J , µ > 0) = 0. The Jacobian is 1/ 1 + e 2 µ 2 and we obtained : 1 + e 2 µ 2 2k+3 dµ dr eq . (32) From Eq.(5), we see that the density only depends on l/l 0 , the mean radius of one isobaric layer above the external mean radius. With Eq.(31), we were then able to write l = β × r eq with the same β for all r eq . The β cancelled in the fraction, and we simply changed l/l 0 by r eq /a J .
We just had to calculate the other part of the integral, with the analytical expression of Eq.(31) : Since there is no simple way to do it, we kept it in our equations, but the interesting aspect is that it no longer depends on r : we have reduced our equation to a simple integral. We were able to check whether this is correct or not. So far, we have supposed that the equipotential surfaces are ellipsoids and that they have the same e value. Considering that we integrate along these surfaces by varying the equatorial radius from the center to a J , we showed that we can write: It was easy to calculate these integrals as we know their analytical expression. We just needed to choose the correct e and m values. From the observed equatorial and polar radii of Jupiter (see Table 1), we have : There is about 20% difference with the theoretical results, which is expected as we have chosen the outermost e value, which is the largest. With e = 0.36, we obtained less than 10% difference on these moments. Considering further approximations that are made in the next appendices, a 0.9 factor is acceptable to compare the theory with the numerical results and we retain this e = 0.36 value in our calculations.
Because there is a high number of layers, we were able to approximate |R i − r | R i everywhere. This approximation is less valid in the ∼ 10 deepest layers, but their impact on the gravitational moments is negligible (Figure 1).
we developped the sinus to order 2 : We know that α ∼ π. In the internal layers, a J R i , so γξ = ξα/a J ξ/R i . On the outside, a J ∼ R i so γξ ∼ ξ/R i . As mentioned above, ξ/R i 1, which means that everywhere γξ 1. Therefore: The zeroth order term cancelled in Eq.(38) so : We introduced This gave a simplified expression for the potential : Developping the power of r yielded: Then, the integral I is given by : We changed the variable of integration r → (r −R i ), and using Eq. (39) : (48) As a first calculation, we chose a constant ∆r with depth, as suggested in H13, ∆r = a J /N . If one wants another repartition of spheroids, they just have to change ∆r by ∆r i depending on depth: (49) This is easily calculated and we got: Puting this expression into Eq.(43), we obtained an intermediate formula : We considered two cases: in the internal layers, where a J R i : In the external layers, a J ≈ R i , there is no obvious approximation. We assumed that C 1 and C 2 were of the same order of magnitude as the coefficient in front of the sinusoidal functions and, remembering that R i /a J ∼ 1: Neglecting the factor α 3 /3 (which varies with height): which yielded the aproximated relation: In Figure 11, we plot the exact and approximated values of C2+ (2k + 2)C1 for k = 1. As expected, in the interior we recover the factor ten shift (α 3 /3) between the real and estimated value, whereas in the externalmost layers this term is smaller so we underestimate the error by a factor ∼ 2 − 5. Everything is thus consistent with our assumptions. With these assumptions, we were then able to write (redefining < P 2k > as its absolute value) : It is important to note that up to now, if ∆r was not constant we would still have the same result by using ∆r i instead of ∆r in the sum. This is done in Appendix C.3 In the case of a linear spacing of the spheroids, ∆r = a J /N , we went further: because, from Eq.(39), A ∼ π/3 if m 1, which is valid for Jupiter (and for the vast majority of celestial bodies). As r i+1 = r i − ∆r and r 0 = a J , r i = a J − i × ∆ r , we get: (60) Combining Eq.(58), Eq.(59) and Eq.(60) gave the final result, as written in Eq. (7) :

C Supplementary calculations for J2
The formula for J 2 is directly derived from Eq. (7) : C.1 Linear spacing First, we wanted to calculate P 2 from Eq.(33). A simple numerical integration yielded (in absolute value): We could have expanded the sum in Eq.(62) with the binomial theorem, but it is easy to show that the first order is equivalent to neglecting the factor 1/2 in the term (i + 1/2). So we had directly: (64) As Eq.(64) corresponds to the 1st order expansion in i in Eq.(62), we needed to make sure that these terms did not cancel in the sum for this approximation to be correct. This sum is If we wanted to keep every term up to this stage, we would find that the second order term is 1/2 N/5. Our approximations were thus justified. Finally, Eq.(62), Eq.(63) and Eq.(65) yielded:

C.2 Hubbard and Militzer spacing
In HM16, the spacing is not exactly linear. The planet is separated in two domains: first, from the outside boundary to r i = 0.5 a J , they use 341 spheroids, and 171 inside this limit. Then, outside ∆r ext = 3∆r/4 and inside ∆r in = 3∆r/2. Moreover, the first spheroid is a bit particular as it has a a zero density over half its size. We did not consider it here: since it is just one spheroid, the impact on the analytical derivation of the error is small. Using Eq.(58), we were able to write it with a general number N as : Which can be rewritten : For k = 1, using Eq.(59) and Eq.(65), we obtained : Compared to Eq.(9), this is approximately twice better in terms of errors.

C.3 Cubic spacing
In this appendix, we calculated the error for a cubic repartition of the spheroids, that is ∆r depending on i as r i = a J − a J × i 3 /N 3 . Using Eq.(58), we got: Straightforwardly, we obtained the new Eq. (7) : which gave for J 2 : (75) We developed again the various terms: In order to calculate the sum, we needed the leading term of but the terms in the sum cancelled two by two : With this formula, it is easy to show by recurrence that . Therefore, in the equation above only the k − 1 term can lead to N k . Then the result: From here, we obtained the leading order terms of the sum of Eq. (78) : The result is thus the same as in the linear case :
As we wanted the increment to be reasonable, e N α −1 could not be too large. In this paper, we chose where we varied γ between 5 and 10, depending on how sharp we wanted the exponential function to be (see Figure 2).
Then, using Eq.(58), we had to calculate, for J 2 : β = e α − 1 e N α − 1 , δ = 1 + 1 e N α − 1 , = 1 + e α 2(e N α − 1) . (93) Developing the power of 4 brackets and using the well known result for the sum of the terms of a geometric sequence: Here, we considered that 7α 1, e N α 1, which is a reasonable approximation in the range γ = 5 − 10 and N > 512. Then : We checked numerically that our approximations yield an error allways smaller than ≈ 10%.
Using the exterior e 0.38 and definingρ S as the density of Jupiter if it was a sphere of radius a 1bar , the neglected potential simply reads: This expression does not depend on r j nor on µ : the first order error is a constant neglected potential on each spheroid. This means that the hydrostatic condition is not perturbed by the neglected high atmosphere, because it is only sensitive to the gradient of the potential. As the 1 bar pressure is set at the observed radius, we do not expect a change in the pressure-density profile of the planet from the direct effect of this neglected potential.
On the other hand, this neglected potential affects the calculation of the shapes of the spheroids. Calling ∆U = φ neg × a 1bar /GM and U j the dimensionless potential of the j th spheroid, we derived the first order perturbation of equation (51) of H13 : For simplicity, we set x =ζ j (µ) the solution without the neglected potential and ∆x the variation to this solution. Both x and ∆x depend on µ. It was then straightforward to write Eq.(104) at the first order : (105) The second bracket could be simplified by the fact that x satisifies equation (51) of H13 : Now, we restricted ourselves to the outermost spheroid j = 0, λ j = 1. Only the first sum remained, and we recoginized Eq. (1) : With Juno data, we know the first few J ext 2k and the contribution of each J 2k to the sum is strongly decreasing with k so we were able to use only the first 4 even gravitational moments. From here, we approximated x given by the CMS program as a polynomial of order 15 (which gave errors of a few 10 −14 ), but allowed us to obtain as many evaluations of ∆x as we needed.
We noticed that around µ ∼ 0.5 the bracketed term tends to 0 so ∆x should not be defined. That means that at this point some neglected effect should be taken into account. However, with the polynomial approximation for x, we were able to use Riemann sum to evaluate the change on equation (40) of H13 : We found that the undefined region for ∆x has a negligible influence on the integral because it is almost ponctual. We also found that the change in the denominator which is, as expected, of the order of the neglected mass, is ten times smaller than the change in the numerator. Therefore, the change in the numerator dominates the uncertainties on J 2 . By comparing for different number of points in the Riemann sum, we obtained : This result has a variation of about a factor 2 (due to the bad handling of the divergence zone) and this evaluation is rather a lower bound for ∆J 2,0 . Because, in the linear limit, the 2k + 3 exponent cancel the 2k + 3 in the denominator, the absolute uncertainty is approximatively the same for all J k (decorrelated by the P 2k ) so the relative uncertainty is a rapidly growing function of k. Typically: ∆J 6,0 J 6,0 6 × 10 3 × ∆U ∆J 8,0 J 8,0 1 × 10 5 × ∆U