A complete study of the precision of the concentric MacLaurin spheroid method to calculate Jupiter’s gravitational moments

F. Debras; G. Chabrier

doi:10.1051/0004-6361/201731682

Home

All issues

Volume 609 (January 2018)

A&A, 609 (2018) A97

Full HTML

Free Access

Issue		A&A Volume 609, January 2018


Article Number		A97
Number of page(s)		19
Section		Planets and planetary systems
DOI		https://doi.org/10.1051/0004-6361/201731682
Published online		23 January 2018

A&A 609, A97 (2018)

A complete study of the precision of the concentric MacLaurin spheroid method to calculate Jupiter’s gravitational moments

F. Debras¹^,2 and G. Chabrier¹^,2

¹ École normale supérieure de Lyon, CRAL, UMR CNRS 5574, 69364 Lyon Cedex 07, France
e-mail: florian_debras@hotmail.com
² School of Physics, University of Exeter, Exeter, EX4 4QL, UK

Received: 31 July 2017
Accepted: 6 October 2017

Abstract

A few years ago, Hubbard (2012, ApJ, 756, L15; 2013, ApJ, 768, 43) 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.

Key words: planets and satellites: gaseous planets / planets and satellites: interiors / equation of state / methods: numerical

© ESO, 2018

1. Introduction

The concentric MacLaurin spheroid (CMS) method has been developed by Hubbard (2012, 2013, hereafter 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 = ω²a³/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ω²/ (4πρG), is also used in the theory of figures (see for example Zharkov & 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 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: $J_{k}^{ext} = \sum_{i = 0}^{N - 1} J_{k}^{i} \times {(λ_{i})}^{k},$ $\begin{equation} J_k^{\rm ext} = \sum_{i=0}^{N-1} J_k^i \times \left(\lambda_i\right)^k, \label{J} \end{equation}$ (1)where $J_{k}^{ext}$ $\hbox{$J_k^{\rm ext}$}$ is the moment of order k, $J_{k}^{i}$ $\hbox{$J_k^{i}$}$ is the moment of order k of the ith spheroid, λ_i is the ratio of the equatorial radius of the ith 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. Equation (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.

Table 1

Values of the planetary parameters of Jupiter.

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 Sect. 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 Sect. 3. We then apply our method with a realistic eos in Sect. 4. These sections show the need to improve the basic CMS method. In Sect. 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

2.1. Evaluation of the errors with 512 spheroids

In order to evaluate the errors analytically, we have considered a polytrope of index n = 1, P ∝ ρ². We call ρ₂ 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: $\begin{matrix} Δ φ & = & - \frac{4 πG}{r} \sum_{k = 0}^{\infty} P_{2 k} (μ) r^{- 2 k} \end{matrix}$ $\begin{eqnarray} \Delta \phi &= &-\dfrac{4\pi G}{r} \sum_{k=0}^{\infty}P_{2k}(\mu) r^{-2k} \nonumber \\ &&\quad\times \int_{0}^{1}\int_{0}^{a_J}\left[\rho (r') - \rho_2(r') \right] r'^{2k+2} P_{2k}(\mu) {\rm d}r' {\rm d}\mu, \label{potential_tot} \end{eqnarray}$ (2)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: $J_{2 k} = - \frac{4 π}{M a_{J}^{2 k}} \int_{0}^{1} \int_{0}^{a_{J}} ρ (r^{'}) r^{' 2 k + 2} P_{2 k} (μ) d r^{'} d μ .$ $\begin{equation} J_{2k} = - \dfrac{4\pi}{Ma_J^{2k}}\int_{0}^{1}\int_{0}^{a_J}\rho (r') r'^{2k+2} P_{2k}(\mu) {\rm d}r' \, {\rm d}\mu. \label{def_J2k} \end{equation}$ (3)By construction we can write ρ₂(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: $\begin{matrix} Δ J_{2 k} & = & - \frac{4 π}{M a_{J}^{2 k}} \sum_{i = 0}^{N - 1} \int_{0}^{1} \int_{r_{i + 1} (μ)}^{r_{i} (μ)} [ρ (r^{'}) - ρ_{i}] \end{matrix}$ $\begin{eqnarray} \Delta J_{2k} &=& - \dfrac{4\pi}{Ma_J^{2k}} \sum_{i=0}^{N-1} \int_{0}^{1}\int_{r_{i+1}(\mu)}^{r_i(\mu)} \left[\rho (r') - \rho_i \right] \nonumber \\ &&\quad \times r'^{2k+2} P_{2k}(\mu) {\rm d}r' \, {\rm d}\mu. \label{eq:def_deltaJ2k} \end{eqnarray}$ (4)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 $\hbox{$\bar{\rho}$}$ is given by (Zharkov & Trubitsyn 1978): $\begin{matrix} with A = \frac{π}{3} (1 - \frac{2}{3} m (1 - \frac{6}{π^{2}})) and α = π (1 + \frac{2 m}{π^{2}}), \end{matrix}$ $\begin{eqnarray} &&\dfrac{\rho(l)}{\bar{\rho}} = A \dfrac{\sin\left(\alpha \dfrac{l}{l_0}\right)}{\left(\dfrac{l}{l_0}\right)} + \dfrac{2}{3} m \label{density_ZT},\\ &&{\rm with}\,\, A = \dfrac{\pi}{3} \left( 1-\dfrac{2}{3}m \left( 1-\dfrac{6}{\pi^2} \right) \right) \text{ and } \alpha = \pi \left(1 + \dfrac{2m}{\pi^2} \right), \end{eqnarray}$ where l₀ is the outer mean radius.

At this stage, our method implied an important approximation, namely that ρ₂ 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$ $\hbox{$\rho_i = \rho_{\rm poly} \left(r_i+r_{i+1}\right)/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: $\begin{matrix} \end{matrix}$ $\begin{eqnarray} &&|\Delta J_{2k}| \sim (2k+3)\dfrac{\pi}{12}\dfrac{\langle P_{2k}\rangle}{N^3} \sum_{i=0}^{N-1} \left(1 - \dfrac{1}{N}\left(i+\dfrac{1}{2} \right)\right)^{2k+2}, \label{J_2k_final} \\ &&{\rm with} ~~\langle P_{2k}\rangle = \bigg |\int_{0}^{1} \left(\dfrac{1}{1+{\rm e}^2\mu^2}\right)^{k+1} P_{2k} (\mu) {\rm d}\mu \bigg |. \label{def_P2k} \end{eqnarray}$ 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 moments than the inner ones. Figure 1 displays the contribution to J₂ of the outer layers of the planet in the polytropic case (we have numerically integrated the expression for J₂ from Eqs. (3), (5) and (A.4)).

Fig. 1

Contribution of each part of the planet to the total value of J₂, normalized to its maximum value, for a polytrope of index n = 1.

To quantify the error, we focused on J₂, the case k = 1. As detailed in Appendix C, one gets: $| Δ J_{2} | \times 10^{6} ~ \frac{π}{12} \frac{⟨ P_{2} ⟩}{N^{2}} \times 10^{6} \approx \frac{9.4 \times 10^{3}}{N^{2}} \cdot$ $\begin{equation} |\Delta J_{2}| \times 10^{6} \sim \dfrac{\pi}{12} \dfrac{\langle P_{2}\rangle }{N^2}\times 10^{6} \approx \dfrac{9.4\times 10^{3}}{N^2}\cdot \label{J_2_fin} \end{equation}$ (9)With N = 512 (as used in H13), we get | ΔJ₂ | × 10⁶ ~ 3.6 × 10^-2. The first data from Juno after three orbits give an uncertainty of 2.72 × 10^-1 on J₂ × 10⁶. 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 shape and the density of each spheroid, and we have neglected the factor α³/ 3 ~ 2−10 in Eq. (B.20). 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.

Hubbard & Militzer (2016, hereafter 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: $\begin{matrix} | Δ J_{2} | \times 10^{6} = \frac{9}{256} \frac{π ⟨ P_{2} ⟩}{N^{2}} \times 10^{6} \approx 1.51 \times 10^{-2} . \end{matrix}$ $\begin{eqnarray} |\Delta J_{2}|\times 10^{6} = \dfrac{9}{256} \dfrac{\pi \langle P_{2}\rangle }{N^2}\times 10^{6} \approx 1.51 \times 10^{-2}. \label{J2_hub} \end{eqnarray}$ (10)This is about 2.5 times better than with the linear spacing. But once again, when considering the factor of ten neglected in Eq. (B.20) 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 Sect. 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₂ 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 Sect. 4.3.

It should be mentioned that Wisdom & Hubbard (2016) calculated similar errors by comparing their respective methods, finding (their Eq. (15)): $\begin{matrix} \log_{10} | \frac{Δ J_{2}}{J_{2}} | ≃ 0.7 - 1.81 \log_{10} (N), \\ whileweget \\ \log_{10} | \frac{Δ J_{2}}{J_{2}} | ≃ - 0.6 - 2 \log_{10} (N) . \end{matrix}$ $\begin{eqnarray} &&\log_{10} \left|\frac{\Delta J_2}{J_2}\right| \simeq 0.7 {-} 1.81 \log_{10}(N),\\ &&\text{while we get} \\ &&\log_{10} \left|\frac{\Delta J_2}{J_2}\right| \simeq -0.6 {-} 2 \log_{10}(N). \end{eqnarray}$ 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 & 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 Sects. 4 and 5.

Another option would be to use Wisdom’s method (Wisdom 1996), which has virtually an unlimited precision. Table 3 of Wisdom & 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₁₀) 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.

2.2. 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₂: $| Δ J_{2}^{cubic} | ~ \frac{π}{12} \frac{⟨ P_{2} ⟩}{N^{2}} \cdot$ $\begin{equation} |\Delta J_{2}^{\rm cubic}| \sim \dfrac{\pi}{12} \dfrac{\langle P_{2}\rangle }{N^2}\cdot \label{J_2k_cubic} \end{equation}$ (14)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. (B.20) 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. (C.11), we got: $\begin{matrix} Δ r_{i}^{cubic} = Δ r \Leftarrow \Rightarrow \frac{Δ i^{3}}{N^{2}} = 1 \\ \Leftarrow \Rightarrow 3 i^{2} + 3 i + 1 = N^{2} \\ = \Rightarrow i ~ \frac{N}{\sqrt{3}} \cdot \end{matrix}$ $\begin{eqnarray} &&\Delta r_i^{{\rm cubic}} = \Delta r \Longleftrightarrow \dfrac{\Delta i^3}{N^2} = 1 \nonumber \\ &&\Longleftrightarrow 3i^2+3i +1 = N^2 \nonumber \\ &&\Longrightarrow i \sim \dfrac{N}{ \sqrt3}\cdot \end{eqnarray}$ (15)More generally, with a spacing of power k, $Δ r_{i}^{k}$ $\hbox{$\Delta r_i^k$}$ is smaller than the linear Δr for $i \leq N / k^{\frac{1}{k - 1}}$ $\hbox{$i \le {N}/{k^{\frac{1}{k-1}}}$}$ , and larger for larger values of i. That means that above $r_{i} = a_{J}^{(} 1 - 1 / k^{\frac{1}{k - 1}}^{)}$ $\hbox{$r_i = a_J \left(1 - 1/k^{\frac{1}{k-1}} \right)$}$ 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 Sect. 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 Sect. 5 is to impose a non zero density as outer boundary condition.

2.3. 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 Fig. 2.

Fig. 2

Equatorial radius as a function of the spheroid number for various functional forms for the repartition of the spheroids. a) from 0 to a_J; b) zoom of the external layers.

For the exponential repartition, the value of the error is: $| Δ J_{2}^{\exp} | \times 10^{6} ≃ \frac{π ⟨ P_{2} ⟩}{252} \times \frac{γ^{2}}{N^{2}} \times 10^{6} \approx 449 {(\frac{γ}{N})}^{2},$ $\begin{equation} |\Delta J_{2}^{\exp}| \times 10^{6} \simeq \dfrac{\pi \langle P_{2}\rangle }{252}\times \dfrac{\gamma^2}{ N^2}\times 10^{6} \approx 449 \left(\dfrac{\gamma}{N}\right)^2, \label{J2_exp} \end{equation}$ (16)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_{2}^{\exp} | \times 10^{6} \in [1.54 \times 10^{-2}, 1.6 \times 10^{1}] .$ $\begin{equation} |\Delta J_{2}^{\exp}|\times 10^{6} \in [1.54\times 10^{-2}, 1.6 \times 10^1]. \end{equation}$ (17)As seen, the smaller γ the smaller the error, but also the larger the range of outer layers which is a problem, as examined in Sect. 3. Note also that the neglected factor of Eq. (B.20) in the analytical calculations is much smaller because of the very large number of spheroids in the outermost part of the envelope (see Fig. 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 Sect. 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 Sect. 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.

3. Numerical calculations

3.1. How to match Juno’s error bars

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.

Table 2

Difference between the values calculated with our code and the results from Hubbard (2013).

Figure 3 displays the difference between the expected and numerical values of J₂ 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.

Fig. 3

Absolute value of the difference between the expected and numerical results on J₂ × 10⁶ 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.

The second and quite important problem, as briefly aluded to at the end of Sect. 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₂ 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₂ × 10⁶ is ~ 5 × 10¹. We discuss that in more details in Sect. 3.2.

Fig. 4

Numerical error as a function of the number of spheroids for the cubic and square repartitions. a) Error between the analytical and measured J₂ × 10⁶; b) difference on J₂ × 10⁶ when the first layer has 0 or non-0 density.

Fig. 5

Numerical error as a function of the number of spheroids for exponential repartitions. a) Error between the analytical and measured J₂ × 10⁶; b) difference on the J₂ × 10⁶ when the first layer has 0 or non-0 density.

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 Fig. 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 (Fig. 4b). 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 (Fig. 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 Fig. 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 Sect. 5.

Fig. 6

Normalized value of $δ ρ_{i} λ_{i}^{5}$ $\hbox{$\delta \rho_i \lambda_i^5$}$ as a function of altitude for linear, cubic and exponential (γ = 8) repartitions of spheroids. a) From center to exterior; b) zoom on the external layers.

As seen in Fig. 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 Fig. 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 Sect. 3.2). Considering Figs. 5a and b, 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 Sect. 5, however, we show that the cubic repartition can also be used, with a different boundary condition.

3.2. 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. Although this could seem quite arbitrary, it has an analytical justification emerging from the linear dependence of the density with radius (Hubbard, priv. comm.). 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₂ × 10⁶ = 13961.23 with an error of 2.7 × 10¹ 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 Fig. 1, we see that the region outside 90% in radius contributes as much to J₂ 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 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 H13 gives the explicit expression of J₂ for each layer. With the notation of the paper: $\begin{matrix} J_{i}^{2} = - (\frac{3}{5}) \frac{δ ρ_{i} \int_{0}^{1} d μ P_{2 k} (μ) ξ_{i} (μ)^{5}}{\sum_{j = 0}^{N - 1} δ ρ_{j} \int_{0}^{1} d μ ξ_{i} (μ)^{3}}, \\ δ ρ_{i} = ρ_{i} - ρ_{i + 1} if i > 0, \\ δ ρ_{0} = ρ_{0}, \\ ξ_{i} (μ) = r_{i} (μ) / a_{J} . \end{matrix}$ $\begin{eqnarray} &&J_i^2 = - \left(\dfrac{3}{5} \right) \dfrac{\delta \rho_i \int _0^1 {\rm d}\mu P_{2k}(\mu) \xi_i(\mu)^5}{\sum_{j=0}^{N-1} \delta \rho_j \int_0^1 {\rm d}\mu \xi_i(\mu)^3}, \\ &&\delta \rho_i = \rho_i - \rho_{i+1} \text{ if } i > 0, \nonumber \\ &&\delta \rho_0 = \rho_0, \\ &&\xi_i(\mu) = r_i(\mu)/a_J. \end{eqnarray}$ To illustrate our argument, we said that every spheroid has the same shape and that δρ_i is the same as in the rest case, given by a slight modification of Eq. (63) of H13: $\begin{matrix} δ ρ_{i} = ρ_{c} \times δ λ_{i} (\frac{\cos (π λ_{i})}{λ_{i}} - \frac{\sin π λ_{i}}{λ_{i}^{2}}), \\ λ_{i} = ξ_{i} (0) and δ λ_{i} = λ_{i} - λ_{i + 1} . \end{matrix}$ $\begin{eqnarray} &&\delta \rho_i = \rho_{\rm c} \times \delta \lambda_i \left( \dfrac{\cos(\pi \lambda_i)}{\lambda_i} - \dfrac{\sin{\pi \lambda_i}}{\lambda_i^2} \right), \\ &&\lambda_i = \xi_i(0) \text{ and } \delta \lambda_i = \lambda_i - \lambda_{i+1}. \end{eqnarray}$ Then, the only different terms for each layer stemmed from the numerator: $δ ρ_{i} \times λ_{i}^{5}$ $\hbox{$\delta \rho_i \times \lambda_i^5$}$ . We plot this value for different spacings of the spheroids in Fig. 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 Fig. 6), but provides a nice cut-off between linear and cubic spacings in the uppermost layer region.

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 Figs. 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 Sect. 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, MacLaurin spheroids imply constant density layers, there is an additional 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 Sect. 5, we examine the solution used in Wahl et al. (2017) to overcome this issue.

3.3. One 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 = 71 492 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 = 71 505 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_{< 1 bar}} × 10⁶ ≃ 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 Sect. 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 Sect. 3 still hold when using a realistic eos.

4. 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).

4.1. 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₂ 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.

Fig. 7

Contribution of each spheroid to the total value of J₂ × 10⁶ and to the truncated J₂ × 10⁶, 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₂ × 10⁶ 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₂.

Table 3

Value of J₂ × 10⁶ for an exponential γ = 7 spacing with 1500 spheroids and the SCvH eos with no core.

We have also tried a solution where we converged the 1 bar radius to the observed R_{1 bar} value, in order to avoid this modification in the structure, but the result still holds: it is not possible 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 Sect. 5.

Table 4

Difference between square, exponential γ = 7 and HM16 spacing for 3 different tests with 1500 spheroids.

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₂ value larger than Juno error bars. A quadratic repartition yields similar results with differences on J₂ × 10⁶ 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 Sect. 3.2: a very accurate evaluation of the first layers is mandatory in order to properly converge the CMS method.

4.2. Errors arising from a (quasi) linear repartition

In order to assess quantitatively the errors arising from an ill-adapted 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₂ 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 = 71 492 km. The results are given in Table 4.

The difference on J₂ × 10⁶ between the square and exponential spacings with 1500 spheroids are < 2 × 10^-1, as expected from Figs. 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₂ × 10⁶ 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₂ × 10⁶ we obtained for the HM16 spacing with 512 spheroids is 15 075.52. When we compare to Table 4, the difference is 5 × 10¹: 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 Fig. 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₂ from the bad handling of the high atmosphere reach almost 100 × Juno’s error bars, as shown in Fig. 8a. On the other hand, when fitting J₂, Fig. 8b shows that the errors on J₄ and J₆ are comparable or inferior to Juno’s. Nonetheless, we show in Sect. 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.

Fig. 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₄ × 10⁶ vs. J₂ × 10⁶ from test 1; b) J₆ × 10⁶ vs. −J₄ × 10⁶ from test 3, when J₂ is fit.

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₂. 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.

4.3. Intrinsic uncertainties on the Jupiter models

In this section, we derived interior models that match both J₂ and J₄ from the first Juno data (Bolton et al. 2017) and determined how they were affected by the number of spheroids and 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 = 20 000 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}^{'}$ $\hbox{$Y^\prime_2$}$ . All the models were converged to the appropriate J₂ and J₄ values with a precision of 10^-8, within Juno’s error bars.

Table 5 gives the values of Y′, $Y_{2}^{'}$ $\hbox{$Y^\prime_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₂ × 10⁶ became 14 647.9, which leads to an error of 5 × 10¹, 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 Sect. 5.

Table 5

Effective helium mass fraction above (Y′) and underneath ( $Y_{2}^{'}$ $\hbox{$Y^\prime_2$}$ ) the entropy jump, and mass of the core (M_c in units of 10²⁴ kg) for different repartitions and numbers of spheroids.

With the quadratic repartition, using the values obtained with 512 spheroids in the 1500 spheroid calculations yields an error on J₂ × 10⁶ of 4 × 10⁰, 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₂ and J₄ 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₂ and J₄. 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 zerodensity first layer, as the values of the moments would be changedby 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.

5. Taking the 1 bar level as the outer boundary condition

5.1. Irreducible errors due to the high atmosphere region (<1 bar)

According to their second Appendix, Wahl et al. (2017) impose the 1 bar level as the outer boundary condition in their CMS calculations. As shown with a polytropic eos in Sect. 3.2, the atmosphere above this level, called here the high atmosphere, has a contribution to J₂ × 10⁶ ≃ 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₂ 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: $\begin{matrix} φ_{neg} ≃ 2.87 \times \frac{GM}{a_{1 bar}} (\frac{ρ_{ext}}{\bar{ρ_{S}}}) (\frac{Δ a}{a_{1 bar}}), \end{matrix}$ $\begin{eqnarray} \phi_{\rm neg} \simeq 2.87\times \dfrac{G M}{a_{\rm 1\,bar}} \left(\dfrac{\rho_{\rm ext}}{\bar{\rho_{\rm S}}} \right) \left(\dfrac{\Delta a}{a_{\rm 1\,bar}} \right), \label{eq:phi_neg} \end{eqnarray}$ (23)where ρ_ext is the (constant) density of the high atmosphere, $\bar{ρ_{S}}$ $\hbox{$\bar{\rho_{\rm S}}$}$ is the mean density of Jupiter if it was a sphere of radius a_{1 bar} = 71 492 km, and Δa is the depth of the high atmosphere.

As mentioned above, Δa ≃ 70 km and we know that $\bar{ρ_{S}} ≃ 10^{3}$ $\hbox{$\bar{\rho_{\rm S}}\simeq 10^3$}$ kg m^-3. Taking the outer density in a range (0.01−1) × ρ_{1 bar}, and the conditions of Jupiter, ρ_{1 bar} ≃ 0.17 kg m^-3, yields: $\begin{matrix} φ_{neg} \in [5 \times 10^{-9}; 5 \times 10^{-7}] \times \frac{GM}{a_{1 bar}} \cdot \end{matrix}$ $\begin{eqnarray} \phi_{\rm neg} \in \left[5 \times 10^{-9} ; 5 \times 10^{-7}\right] \times \dfrac{G M}{a_{\rm 1\,bar}}\cdot \end{eqnarray}$ (24)This is clearly a first order correction since the potential on each spheroid, given by the CMS method, ranges from GM/a_{1 bar} in the outermost layers to 2GM/a_{1 bar} 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 Eq. (40) of H13 and obtained the following result: $\frac{Δ J_{2, 0}}{J_{2, 0}} ≃ 287 \times (\frac{ρ_{ext}}{\bar{ρ_{S}}}) (\frac{Δ a}{a_{1 bar}}) \cdot$ $\begin{equation} \dfrac{\Delta J_{2,0}}{J_{2,0}} \simeq 287 \times \left(\dfrac{\rho_{\rm ext}}{\bar{\rho_{\rm S}}} \right) \left(\dfrac{\Delta a}{a_{\rm 1\,bar}} \right)\cdot \end{equation}$ (25)Within the above range of outer densities we get: $\frac{Δ J_{2, 0}}{J_{2, 0}} \in [5 \times 10^{-7}; 5 \times 10^{-5}] .$ $\begin{equation} \dfrac{\Delta J_{2,0}}{J_{2,0}} \in \left[5 \times 10^{-7};5 \times 10^{-5}\right]. \label{eq:delta_J2,0} \end{equation}$ (26)Since the potential is smallest in the outside layers, we expect the relative change in J_2,i of the ith 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₂, 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_{2}^{ext} \times 10^{6} ≃ 8 \times 10^{-6} (J_{2}^{ext} \times 10^{6}) ≃ 1 \times 10^{-1} .$ $\begin{equation} \Delta J_{2}^{\rm ext}\times 10^6 \simeq 8 \times 10^{-6} (J_{2}^{\rm ext} \times 10^6) \simeq 1 \times 10^{-1}. \label{eq:delta_J2_highatm} \end{equation}$ (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₂. From the evaluation of J₄ and J₆ in Appendix D, we calculated that: $\begin{matrix} Δ J_{4}^{ext} \times 10^{6} ≃ 6 \times 10^{-2}, \end{matrix}$ $\begin{eqnarray} &&\Delta J_{4}^{\rm ext}\times 10^6 \simeq 6 \times 10^{-2}, \\ &&\Delta J_{6}^{\rm ext}\times 10^6 \simeq 2 \times 10^{-2}. \label{eq:delta_Jn_highatm} \end{eqnarray}$ We see that the agreement with Juno’s error bars improves with higher order moments, although it stabilizes around 10^-2 from J₆ 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 Fig. 8, we show graphically these uncertainties in Fig. 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 Fig. 9.

Fig. 9

Comparison of the errors of the CMS arising from neglecting the high atmosphere with Juno’s error bar, shifted to be centered on the observed J_n. a) −J₄ × 10⁶ vs. J₂ × 10⁶; b) J₆ × 10⁶ vs. −J₄ × 10⁶.

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₂ × 10⁶ 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. (D.11) show that, globally, this error is of the order of 10^-2 for higher order moments. Therefore, the J_k × 10⁶ cannot be determined within better than a few 10^-2. On J₁₀ 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).

5.2. 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 Sect. 3 are affected or not by the change of outer boundary condition. Figure 10 displays the dependence of the error upon the number of spheroids for different spheroid repartitions in the present case.

Fig. 10

Numerical error as a function of the number of spheroids when neglecting the high atmosphere for a polytropic eos.

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, however, this needs to be clearly specified when deriving Jupiter models in the context of Juno to verify their validity.

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 ρ_{1 bar} ≃ 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 Sect. 4.3, we give in Table 6 the values of Y′, $Y_{2}^{'}$ $\hbox{$Y_2'$}$ and M_c obtained when matching J₂ and J₄ 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 71 492 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 Sect. 4.

Table 6

Same as Table 5 when neglecting the high atmosphere.

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₂ × 10⁶ 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₂ × 10⁶ and removed it from J₄ × 10⁶. We stress again that this must be done because the value of both J₂ and J₄ 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 Tables 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 7

Same as Table 6 except that J₂ × 10⁶ and J₄ × 10⁶ have been changed by ± 1 × 10^-1.

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.

6. 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 & 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 on the determination of J₂ 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 Sect. 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₂, a small change in the density structure of these regions leads eventually to unacceptable changes in the total value of J₂. Furthermore, when fitting J₂ and J₄, 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₂ × 10⁶ 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.

Acknowledgments

The authors are deeply indebted to W. Hubbard for kindly providing a version of his CMS code, which enabled us to carefully compare our results with the published ones, and for useful conversations.

References

Archinal, B. A., A’Hearn, M. F., Bowell, E., et al. 2011, Celes. Mech. Dyn. Astron., 109, 101 [Google Scholar]
Bolton, S. J., Adriani, A., Adumitroaie, V., et al. 2017, Science, 356, 821 [NASA ADS] [CrossRef] [Google Scholar]
Chabrier, G., Saumon, D., Hubbard, W. B., & Lunine, J. I. 1992, ApJ, 391, 817 [NASA ADS] [CrossRef] [Google Scholar]
Cohen, E. R., & Taylor, B. N. 1987, Rev. Mod. Phys., 59, 1121 [NASA ADS] [CrossRef] [Google Scholar]
Folkner, W. M., Iess, L., Anderson, J. D., et al. 2017, Geophys. Res. Lett., 44, 4694 [NASA ADS] [CrossRef] [Google Scholar]
Hubbard, W. B. 1975, Sov. Astron., 18, 621 [NASA ADS] [Google Scholar]
Hubbard, W. B. 2012, ApJ, 756, L15 [NASA ADS] [CrossRef] [Google Scholar]
Hubbard, W. B. 2013, ApJ, 768, 43 [NASA ADS] [CrossRef] [Google Scholar]
Hubbard, W. B., & Militzer, B. 2016, ApJ, 820, 80 [NASA ADS] [CrossRef] [Google Scholar]
Hubbard, W. B., Schubert, G., Kong, D., & Zhang, K. 2014, Icarus, 242, 138 [NASA ADS] [CrossRef] [Google Scholar]
Nettelmann, N. 2017, A&A, 606, A139 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Riddle, A. C., & Warwick, J. W. 1976, Icarus, 27, 457 [NASA ADS] [CrossRef] [Google Scholar]
Saumon, D., Chabrier, G., & van Horn, H. M. 1995, ApJS, 99, 713 [NASA ADS] [CrossRef] [Google Scholar]
Wahl, S. M., Hubbard, W. B., Militzer, B., et al. 2017, Geophys. Res. Lett., 44, 4649 [NASA ADS] [CrossRef] [Google Scholar]
Wisdom, J. 1996, http://web.mit.edu/wisdom/www/interior.pdf [Google Scholar]
Wisdom, J., & Hubbard, W. B. 2016, Icarus, 267, 315 [NASA ADS] [CrossRef] [Google Scholar]
Zharkov, V. N., & Trubitsyn, V. P. 1978, The physics of planetary interiors (Tuscon: Pachart Publishing House) [Google Scholar]

Appendix 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: $II = \int_{0}^{1} \int_{0}^{a_{J}} ρ (r^{'}) r^{' 2 k + 2} P_{2 k} (μ) d r^{'} d μ .$ $\appendix \setcounter{section}{1} \begin{equation} II = \int_{0}^{1}\int_{0}^{a_J} \rho (r') r'^{2k+2} P_{2k}(\mu) {\rm d}r' \, {\rm d}\mu . \end{equation}$ (A.1)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′(μ)) = const. 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: $r^{'} (μ) = \frac{r_{eq}}{\sqrt{1 + e^{2} μ^{2}}},$ $\appendix \setcounter{section}{1} \begin{equation} r'(\mu) = \dfrac{r_{\rm eq}}{\sqrt{1+{\rm e}^2\mu^2}}, \label{r_of_mu} \end{equation}$ (A.2)where r_eq is the equatorial radius of the spheroid. In reality, this is not exactly the case, the isobaric spheroids 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^{'} \to r_{eq} / \sqrt{1 + e^{2} μ^{2}}$ $\hbox{$r' \rightarrow r_{\rm eq}/\!\sqrt{1+{\rm e}^2\mu^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 / \sqrt{1 + e^{2} μ^{2}}$ $\hbox{$1/\sqrt{1+{\rm e}^2 \mu^2}$}$ and we obtained: $II = \int_{0}^{a_{J}} ρ (r_{eq}) r_{eq}^{2 k + 2} \int_{0}^{1} \frac{P_{2 k} (μ)}{{(\sqrt{1 + e^{2} μ^{2}})}^{2 k + 3}} d μ d r_{eq} .$ $\appendix \setcounter{section}{1} \begin{equation} II = \int_{0}^{a_J} \rho (r_{\rm eq})r_{\rm eq}^{2k+2} \int_{0}^{1} \dfrac{P_{2k}(\mu)}{\left(\sqrt{1+{\rm e}^2\mu^2}\right)^{2k+3}} {\rm d}\mu\, {\rm d}r_{\rm eq}. \end{equation}$ (A.3)From Eq. (5), we see that the density only depends on l/l₀, the mean radius of one isobaric layer above the external mean radius. With Eq. (A.2), 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₀ by r_eq/a_J.

We just had to calculate the other part of the integral, with the analytical expression of Eq. (A.2): $⟨ P_{2 k} ⟩ = \int_{0}^{1} {(\frac{1}{1 + e^{2} μ^{2}})}^{\frac{2 k + 3}{2}} P_{2 k} (μ) d μ .$ $\appendix \setcounter{section}{1} \begin{equation} \langle P_{2k}\rangle = \int_{0}^{1} \left(\dfrac{1}{1+{\rm e}^2\mu^2}\right)^{\frac{2k+3}{2}} P_{2k}(\mu) {\rm d}\mu. \label{def_<P>} \end{equation}$ (A.4)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: $\begin{matrix} J_{2 k} ~ & - & \frac{4 π}{M a_{J}^{2 k}} (\int_{0}^{1} {(\frac{1}{1 + e^{2} μ^{2}})}^{\frac{2 k + 3}{2}} P_{2 k} (μ) d μ) \\ \times \end{matrix}$ $\appendix \setcounter{section}{1} \begin{eqnarray} J_{2k} \sim &\!\!-\!\!& \dfrac{4\pi}{Ma_J^{2k}} \left( \int_{0}^{1} \left(\dfrac{1}{1+{\rm e}^2\mu^2}\right)^{\frac{2k+3}{2}} P_{2k}(\mu) {\rm d}\mu \right) \nonumber \\ &\!\!\times\!\!& \int_0^{a_J} \rho(r_{\rm eq}) r_{\rm eq}^{2k+2} {\rm d}r_{\rm eq}. \text{ } \label{eval_J2} \end{eqnarray}$ (A.5)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: $e = \sqrt{71 492^{2} / 66854^{2} - 1} \approx 0.378897.$ $\appendix \setcounter{section}{1} \begin{equation} e = \sqrt{71\,492^2/66854^2-1} \approx 0.378897. \end{equation}$ (A.6)These radii also give us the value for the mean density: $ρ̅ = \frac{M_{J}}{\frac{4}{3} π a_{J}^{2} (a_{J})_{polaire}} \approx 1326.5$ $\hbox{$\bar{\rho} = \dfrac{M_J}{\dfrac{4}{3}\pi a_J^2 (a_J)_{\rm polaire}} \approx 1326.5$}$ kg m^-3, which allowed us to calculate $m = \frac{3 ω^{2}}{4 πG ρ̅} \approx 0.083408.$ $\appendix \setcounter{section}{1} \begin{equation} m = \dfrac{3 \omega^2}{4 \pi G \bar{\rho}} \approx 0.083408. \end{equation}$ (A.7)Implemented in Eq. (A.5) with the expression of density from Eq. (5), a numerical integration gave: $\begin{matrix} J_{2} \times 10^{6} = 16348.804 with J_{2}^{theory} \times 10^{6} = 13988.511 \\ - J_{4} \times 10^{6} = 643.983 with - J_{4}^{theory} \times 10^{6} = 531.828 \\ J_{6} \times 10^{6} = 35.16 with J_{6}^{theory} \times 10^{6} = 30.12. \end{matrix}$ $\appendix \setcounter{section}{1} \begin{eqnarray} &&J_2 \times 10^6 = 16348.804 \text{ with } J_2^{\rm theory} \times 10^6 = 13988.511 \nonumber \\ &&-J_4\times 10^6 = 643.983 \text{ with } -J_4^{\rm theory}\times 10^6 = 531.828 \nonumber \\ &&J_6\times 10^6 = 35.16 \text{ with } J_6^{\rm theory}\times 10^6 = 30.12. \end{eqnarray}$ (A.8)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.

Appendix B: Supplementary calculations for the general case

Going back to Eqs. (4), (5) and (A.4) we have: $\begin{matrix} | Δ J_{2 k} | ~ | \frac{4 π}{M a_{J}^{2 k}} \sum_{i = 0}^{N - 1} ⟨ P_{2 k} ⟩ \end{matrix}$ $\appendix \setcounter{section}{2} \begin{eqnarray} &&\left| \Delta J_{2k} \right| \sim \bigg| \dfrac{4\pi}{Ma_J^{2k}}\sum_{i=0}^{N-1} \langle P_{2k}\rangle \nonumber \\ &&\qquad \quad\times \int_{r_{i+1}}^{r_i} \bar{A}\left[\dfrac{\sin(\alpha \dfrac{r'}{a_J})}{\dfrac{r'}{a_J}} - \dfrac{\sin(\alpha \dfrac{R_i}{a_J})}{\dfrac{R_i}{a_J}} \right] r'^{2k+2} {\rm d}r' \label{J2k_maths} \bigg|,\\ &&{\rm with}\,\,R_i = \dfrac{r_{i+1}+r_{i}}{2}, \label{def_RA}\\ &&{\rm and}\,\,\bar{A} = \bar{\rho} A \label{def_Abar}. \end{eqnarray}$ 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 (Fig. 1).

With: $γ = \frac{α}{a_{J}},$ $\appendix \setcounter{section}{2} \begin{equation} \gamma = \dfrac{\alpha}{a_J}, \end{equation}$ (B.4)we developped the sinus to order 2: $\begin{matrix} \frac{\sin (γ r^{'})}{r^{'}} = \frac{\sin (γ R_{i} (1 + \frac{r^{'} - R_{i}}{R_{i}}))}{R_{i} (1 + \frac{r^{'} - R_{i}}{R_{i}})} \\ \frac{\sin (γ r^{'})}{r^{'}} ≃ \frac{1}{R_{i}} (1 - \frac{ξ}{R_{i}} + {(\frac{ξ}{R_{i}})}^{2}) (\sin (γ R_{i}) \cos (γξ) \\ + \cos (γ R_{i}) \sin (γξ)), \\ where ξ = r^{'} - R_{i} ≪ R_{i} . \end{matrix}$ $\appendix \setcounter{section}{2} \begin{eqnarray} &&\dfrac{\sin\left( \gamma r' \right)}{r'} = \dfrac{\sin\left( \gamma R_i\left(1 + \dfrac{r'-R_i}{R_i} \right) \right)} {R_i \left(1+\dfrac{r'-R_i}{R_i} \right)} \nonumber \\ &&\dfrac{\sin\left( \gamma r' \right)}{r'} \simeq \dfrac{1}{R_i} \left(1-\dfrac{\xi}{R_i} +\left(\dfrac{\xi}{R_i}\right)^2 \right) \bigg(\sin(\gamma R_i)\cos(\gamma \xi) \nonumber \\ &&\qquad \qquad \,+ \cos(\gamma R_i) \sin(\gamma \xi) \bigg), \nonumber \\ &&\text{ where } \xi = r'-R_i \ll R_i. \nonumber \end{eqnarray}$ 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: $\begin{matrix} \frac{\sin (γ r^{'})}{r^{'}} & ≃ & \frac{1}{R_{i}} (1 - \frac{ξ}{R_{i}} + {(\frac{ξ}{R_{i}})}^{2}) (\sin (γ R_{i}) (1 - \frac{(γξ)^{2}}{2}) \\ + \cos (γ R_{i}) (γξ)), \\ \frac{\sin (γ r^{'})}{r^{'}} & ≃ & \frac{\sin (γ R_{i})}{R_{i}} + ξ (\frac{γ}{R_{i}} \cos (γ R_{i}) - \frac{1}{R_{i}^{2}} \sin (γ R_{i})) \\ + ξ^{2} (- \frac{γ^{2}}{2 R_{i}} \sin (γ R_{i}) - \frac{γ}{R_{i}^{2}} \cos (γ R_{i}) + \frac{1}{R_{i}^{3}} \sin (γ R_{i})) . \end{matrix}$ $\appendix \setcounter{section}{2} \begin{eqnarray} \dfrac{\sin\left( \gamma r' \right)}{r'} &\simeq& \dfrac{1}{R_i} \left(1-\dfrac{\xi}{R_i} +\left(\dfrac{\xi}{R_i}\right)^2 \right) \bigg(\sin(\gamma R_i)\left(1 - \dfrac{(\gamma \xi)^2}{2} \right) \nonumber \\ &&\quad+ \cos(\gamma R_i) \left( \gamma\xi \right)\bigg), \nonumber \\ \dfrac{\sin\left( \gamma r' \right)}{r'} &\simeq &\dfrac{\sin(\gamma R_i)}{R_i} + \xi \left( \dfrac{\gamma}{R_i} \cos(\gamma R_i) - \dfrac{1}{R_i^2}\sin(\gamma R_i) \right) \nonumber \\ &&\quad +\xi^2 \left( - \dfrac{\gamma^2}{2R_i}\sin(\gamma R_i) - \dfrac{\gamma}{R_i^2} \cos(\gamma R_i) + \dfrac{1}{R_i^3} \sin(\gamma R_i) \right). \end{eqnarray}$ (B.5)The zeroth order term cancelled in Eq. (B.1) so: $\begin{matrix} | Δ J_{2 k} | & ~ & | \frac{4 π}{M a_{J}^{2 k}} \sum_{i = 0}^{N - 1} ⟨ P_{2 k} ⟩ \frac{A̅ a_{J}}{R_{i}} \\ \times \int_{r_{i + 1}}^{r_{i}} [(r^{'} - R_{i}) (γ \cos (γ R_{i}) - \frac{1}{R_{i}} \sin (γ R_{i})) \\ + (r^{'} - R_{i})^{2} (- \frac{γ^{2}}{2} \sin (γ R_{i}) - \frac{γ}{R_{i}} \cos (γ R_{i}) \end{matrix}$ $\appendix \setcounter{section}{2} \begin{eqnarray} |\Delta J_{2k}| &\sim &\bigg| \dfrac{4\pi}{Ma_J^{2k}} \displaystyle{\sum_{i=0}^{N-1}} \langle P_{2k}\rangle\dfrac{\bar{A}a_J}{R_i} \nonumber \\ &&\quad\times \displaystyle{\int_{r_{i+1}}^{r_i}} \Bigg[ (r' - R_i) \left( \gamma \cos(\gamma R_i) - \dfrac{1}{R_i}\sin(\gamma R_i) \right) \nonumber \\ &&\quad+ (r' - R_i)^2 \Bigg( - \dfrac{\gamma^2}{2}\sin(\gamma R_i) - \dfrac{\gamma}{R_i} \cos(\gamma R_i) \nonumber \\ &&\quad + \dfrac{1}{R_i^2} \sin(\gamma R_i) \Bigg) \Bigg]r'^{2k+2} {\rm d}r'\bigg|. \label{potential} \end{eqnarray}$ (B.6)We introduced $\begin{matrix} C_{1} = (γ \cos (γ R_{i}) - \frac{1}{R_{i}} \sin (γ R_{i})) \\ C_{2} = (- \frac{γ^{2}}{2} \sin (γ R_{i}) - \frac{γ}{R_{i}} \cos (γ R_{i}) + \frac{1}{R_{i}^{2}} \sin (γ R_{i})) . \end{matrix}$ $\appendix \setcounter{section}{2} \begin{eqnarray} &&C_1 = \left( \gamma \cos(\gamma R_i) - \dfrac{1}{R_i}\sin(\gamma R_i) \right) \nonumber \\ &&C_2 = \left( - \dfrac{\gamma^2}{2}\sin(\gamma R_i) - \dfrac{\gamma}{R_i} \cos(\gamma R_i) + \dfrac{1}{R_i^2} \sin(\gamma R_i) \right). \end{eqnarray}$ (B.7)This gave a simplified expression for the potential: $\begin{matrix} | Δ J_{2 k} | & ~ & | \frac{4 π}{M a_{J}^{2 k}} \sum_{i = 0}^{N - 1} ⟨ P_{2 k} ⟩ \frac{A̅ a_{J}}{R_{i}} \\ \times \int_{r_{i + 1}}^{r_{i}} (C_{1} (r^{'} - R_{i}) + C_{2} (r^{'} - R_{i})^{2}) r^{' 2 k + 2} d r^{'} | . \end{matrix}$ $\appendix \setcounter{section}{2} \begin{eqnarray} |\Delta J_{2k}| &\sim& \bigg| \dfrac{4\pi}{Ma_J^{2k}} \sum_{i=0}^{N-1} \langle P_{2k}\rangle\dfrac{\bar{A}a_J}{R_i} \nonumber \\ &&\quad \times \displaystyle{\int_{r_{i+1}}^{r_i}} \left (C_1 (r' - R_i) + C_2 (r' - R_i)^2 \right)r'^{2k+2} {\rm d}r' \bigg|. \end{eqnarray}$ (B.8)Developping the power of r′ yielded: $\begin{matrix} r^{' 2 k + 2} = R_{i}^{2 k + 2} {(1 + \frac{r^{'} - R_{i}}{R_{i}})}^{2 k + 2}, \\ r^{' 2 k + 2} ≃ R_{i}^{2 k + 2} (1 + (2 k + 2) \frac{r^{'} - R_{i}}{R_{i}}) \cdot \end{matrix}$ $\appendix \setcounter{section}{2} \begin{eqnarray} &&r'^{2k+2} = R_i^{2k+2} \left(1+\dfrac{r'-R_i}{R_i} \right)^{2k+2}, \nonumber \\ &&r'^{2k+2} \simeq R_i^{2k+2} \left(1 + (2k+2)\dfrac{r'-R_i}{R_i} \right)\cdot \end{eqnarray}$ (B.9)Then, the integral I is given by: $\begin{matrix} I & ≃ & R_{i}^{2 k + 2} \int_{r_{i + 1}}^{r_{i}} [C_{1} (r^{'} - R_{i}) \\ + (C_{2} + (2 k + 2) \frac{C_{1}}{R_{i}}) (r^{'} - R_{i})^{2}] d r^{'} . \end{matrix}$ $\appendix \setcounter{section}{2} \begin{eqnarray} I &\simeq& R_i^{2k+2} \displaystyle{\int_{r_{i+1}}^{r_i}} \bigg[C_1 (r' - R_i) \nonumber \\ &&\quad+\left( C_2 + (2k+2) \dfrac{C_1}{R_i} \right) (r' - R_i)^2 \bigg] {\rm d}r'. \end{eqnarray}$ (B.10)We changed the variable of integration r′ → (r′−R_i), and using Eq. (B.2): $I ≃ R_{i}^{2 k + 2} \int_{\frac{r_{i + 1} - r_{i}}{2}}^{\frac{r_{i} - r_{i + 1}}{2}} [C_{1} r^{'} + (C_{2} + (2 k + 2) \frac{C_{1}}{R_{i}}) r^{' 2}] d r^{'} .$ $\appendix \setcounter{section}{2} \begin{equation} I \simeq R_i^{2k+2} \displaystyle{\int_{\frac{r_{i+1}-r_i}{2}}^{\frac{r_i-r_{i+1}}{2}}} \left [C_1 r' + \left( C_2 + (2k+2) \dfrac{C_1}{R_i} \right) r'^2 \right] {\rm d}r'. \end{equation}$ (B.11)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: $I ≃ R_{i}^{2 k + 2} \int_{- \frac{Δ r}{2}}^{\frac{Δ r}{2}} [C_{1} r^{'} + (C_{2} + (2 k + 2) \frac{C_{1}}{R_{i}}) r^{' 2}] d r^{'} .$ $\appendix \setcounter{section}{2} \begin{equation} I \simeq R_i^{2k+2} \displaystyle{\int_{-\frac{\Delta r}{2}}^{\frac{\Delta r}{2}}} \left [C_1 r' + \left( C_2 + (2k+2) \dfrac{C_1}{R_i} \right) r'^2 \right] {\rm d}r'. \end{equation}$ (B.12)This is easily calculated and we got: $\begin{matrix} I ≃ R_{i}^{2 k + 2} (C_{2} + (2 k + 2) \frac{C_{1}}{R_{i}}) \frac{Δ r^{3}}{12} \cdot \end{matrix}$ $\appendix \setcounter{section}{2} \begin{eqnarray} I \simeq R_i^{2k+2} \, \left( C_2 + (2k+2) \dfrac{C_1}{R_i} \right) \dfrac{\Delta r^3}{12}\cdot \end{eqnarray}$ (B.13)Puting this expression into Eq. (B.6), we obtained an intermediate formula: $\begin{matrix} | Δ J_{2 k} | & ~ & | \frac{4 π}{M a_{J}^{2 k}} \sum_{i = 0}^{N - 1} ⟨ P_{2 k} ⟩ A̅ a_{J} R_{i}^{2 k + 1} \end{matrix}$ $\appendix \setcounter{section}{2} \begin{eqnarray} |\Delta J_{2k}| &\sim& \bigg| \dfrac{4\pi}{Ma_J^{2k}} \sum_{i=0}^{N-1} \langle P_{2k}\rangle\bar{A}a_J R_i^{2k+1} \nonumber \\ &&\quad \times \left( C_2 + (2k+2) \dfrac{C_1}{R_i} \right) \dfrac{\Delta r^3}{12} \bigg|\cdot \label{moments_maths} \end{eqnarray}$ (B.14)We considered two cases: in the internal layers, where a_J ≫ R_i: $\begin{matrix} C_{1} ~ γ - \frac{γ^{3} R_{i}^{2}}{2} - γ + \frac{γ^{3} R_{i}^{2}}{6} = - \frac{α^{3}}{3} \frac{R_{i}^{2}}{a_{j}^{3}}, \\ C_{2} ~ - \frac{γ^{3} R_{i}}{2} - \frac{γ}{R_{i}} + \frac{γ^{3} R_{i}}{2} + \frac{γ}{R_{i}} - \frac{γ^{3} R_{i}}{6} = - \frac{α^{3}}{6} \frac{R_{i}}{a_{j}^{3}} \cdot \end{matrix}$ $\appendix \setcounter{section}{2} \begin{eqnarray} &&C_1 \sim \gamma -\dfrac{\gamma ^3 R_i^2}{2} - \gamma + \dfrac{\gamma^3 R_i^2}{6} = -\dfrac{\alpha^3}{3} \dfrac{R_i^2}{a_j^3}, \\ &&C_2 \sim -\dfrac{\gamma^3 R_i}{2} - \dfrac{\gamma}{R_i} + \dfrac{\gamma^3 R_i}{2} + \dfrac{\gamma}{R_i} - \dfrac{\gamma^3 R_i}{6} = - \dfrac{\alpha^3}{6}\dfrac{R_i}{a_j^3}\cdot \end{eqnarray}$ In the external layers, a_J ≈ R_i, there is no obvious approximation. We assumed that C₁ and C₂ were of the same order of magnitude as the coefficient in front of the sinusoidal functions and, remembering that R_i/a_J ~ 1: $\begin{matrix} C_{1} ~ γ ~ \frac{1}{R_{i}} = {(\frac{R_{i}}{a_{j}})}^{2} \frac{1}{a_{J}} = \frac{R_{i}^{2}}{a_{j}^{3}}, \\ C_{2} ~ γ^{2} ~ \frac{1}{R_{i}^{2}} = \frac{R_{i}}{a_{J}} \frac{1}{a_{J}^{2}} = \frac{R_{i}}{a_{j}^{3}} \cdot \end{matrix}$ $\appendix \setcounter{section}{2} \begin{eqnarray} &&C_1 \sim \gamma \sim \dfrac{1}{R_i} = \left(\dfrac{R_i}{a_j}\right)^2 \dfrac{1}{a_J} = \dfrac{R_i^2}{a_j^3}, \\ &&C_2 \sim \gamma^2 \sim \dfrac{1}{R_i^2} = \dfrac{R_i}{a_J}\dfrac{1}{a_J^2} = \dfrac{R_i}{a_j^3}\cdot \end{eqnarray}$ Neglecting the factor α³/ 3 (which varies with height): $\begin{matrix} C_{1} ~ \frac{R_{i}^{2}}{a_{j}^{3}} \\ C_{2} ~ \frac{R_{i}}{a_{j}^{3}}, \end{matrix}$ $\appendix \setcounter{section}{2} \begin{eqnarray} &&C_1 \sim \dfrac{R_i^2}{a_j^3} \nonumber \\ &&C_2 \sim \dfrac{R_i}{a_j^3}, \end{eqnarray}$ (B.19)which yielded the aproximated relation: $\begin{matrix} C_{2} + (2 k + 2) \frac{C_{1}}{R_{i}} ~ (2 k + 3) \frac{R_{i}}{a_{J}^{3}} \cdot \end{matrix}$ $\appendix \setcounter{section}{2} \begin{eqnarray} C_2 + (2k+2) \dfrac{C_1}{R_i} \sim (2k+3)\dfrac{R_i}{a_J^3}\cdot \label{C1_C2} \end{eqnarray}$ (B.20)In Fig. B.1, 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) 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.

Fig. B.1

Normalized value of C2 + 4C1 /R_i and of its approximation $5 R_{i} / a_{J}^{3}$ $\hbox{$5R_i/a_J^3$}$ (Eq. (B.20)) as a function of altitude R_i. a) From center to exterior; b) zoom on the external layers

With these assumptions, we were then able to write (redefining ⟨ P_2k ⟩ as its absolute value): $| Δ J_{2 k} | ~ 4 π (2 k + 3) \frac{Δ r^{3}}{12} \frac{A̅ ⟨ P_{2 k} ⟩}{M} \sum_{i = 0}^{N - 1} {(\frac{R_{i}}{a_{J}})}^{2 k + 2} \cdot$ $\appendix \setcounter{section}{2} \begin{equation} |\Delta J_{2k}| \sim 4\pi(2k+3)\dfrac{\Delta r^3}{12}\dfrac{\bar{A}\langle P_{2k}\rangle}{M} \sum_{i=0}^{N-1} \left( \dfrac{R_i}{a_J}\right)^{2k+2}\cdot \label{J_2k_deltar} \end{equation}$ (B.21)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: $\frac{Δ r^{3} A̅}{M} = \frac{a_{J}^{3}}{N^{3}} \frac{A ρ̅}{M} = (\frac{a_{J}^{3} ρ̅}{M}) \frac{A}{N^{3}} = \frac{3}{4 π} \frac{A}{N^{3}} ~ \frac{1}{4 N^{3}},$ $\appendix \setcounter{section}{2} \begin{equation} \dfrac{\Delta r^3\bar{A}}{M} = \dfrac{a_J^3}{N^3}\dfrac{A \bar{\rho}}{M} = \left(\dfrac{a_J^3 \bar{\rho}}{M} \right)\dfrac{A}{N^3} =\dfrac{3}{4 \pi}\dfrac{A}{N^3} \sim \dfrac{1}{4N^3}, \label{simp_deltar} \end{equation}$ (B.22)because, from Eq. (B.2), 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₀ = a_J, r_i = a_J−i × Δ_r, we get: $R_{i} = \frac{r_{i} + r_{i + 1}}{2} = a_{J} - (i + \frac{1}{2}) Δ r = a_{J} (1 - (i + \frac{1}{2}) / N) .$ $\appendix \setcounter{section}{2} \begin{equation} R_i = \dfrac{r_i+r_{i+1}}{2} = a_J - \left(i+\dfrac{1}{2} \right) \Delta r = a_J \left(1 - \left(i+\dfrac{1}{2} \right)\bigg/N\right). \label{simp_Ri} \end{equation}$ (B.23)Combining Eqs. (B.21), (B.22) and (B.23) gave the final result, as written in Eq. (7): $| Δ J_{2 k} | ~ (2 k + 3) \frac{π}{12} \frac{⟨ P_{2 k} ⟩}{N^{3}} \sum_{i = 0}^{N - 1} {(1 - (i + \frac{1}{2}) / N)}^{2 k + 2} .$ $\appendix \setcounter{section}{2} \begin{equation} |\Delta J_{2k}| \sim (2k+3)\dfrac{\pi}{12}\dfrac{\langle P_{2k}\rangle}{N^3} \sum_{i=0}^{N-1} \left(1 - \left(i+\dfrac{1}{2} \right)\bigg/N\right)^{2k+2} . \end{equation}$ (B.24)

Appendix C: Supplementary calculations for J₂

The formula for J₂ is directly derived from Eq. (7): $| Δ J_{2} | ~ \frac{5 π}{12} \frac{⟨ P_{2} ⟩}{N^{3}} \sum_{i = 0}^{N - 1} {(1 - (i + \frac{1}{2}) / N)}^{4} .$ $\appendix \setcounter{section}{3} \begin{equation} |\Delta J_{2}| \sim \dfrac{5\pi}{12}\dfrac{\langle P_{2}\rangle}{N^3} \sum_{i=0}^{N-1} \left(1 - \left(i+\dfrac{1}{2} \right)\bigg/N\right)^{4}. \label{J_2} \end{equation}$ (C.1)

Appendix C.1: Linear spacing

First, we wanted to calculate ⟨ P₂ ⟩ from Eq. (A.4). A simple numerical integration yielded (in absolute value): $\begin{matrix} ⟨ P_{2} ⟩ \approx 0.035982. \end{matrix}$ $\appendix \setcounter{section}{3} \begin{eqnarray} \langle P_{2}\rangle \approx 0.035982. \label{P_2_moyen} \end{eqnarray}$ (C.2)We could have expanded the sum in Eq. (C.1) 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: $P = {(1 - (i + \frac{1}{2}) / N)}^{4} ~ 1 - \frac{4}{N} i + \frac{6}{N^{2}} i^{2} - \frac{4}{N^{3}} i^{3} + \frac{i^{4}}{N^{4}} \cdot$ $\appendix \setcounter{section}{3} \begin{equation} P = \left(1 - \left(i+\dfrac{1}{2} \right)/N\right)^{4} \sim 1 - \dfrac{4}{N}i + \dfrac{6}{N^2}i^2 - \dfrac{4}{N^3}i^3 + \dfrac{i^4}{N^4} \cdot \label{P_dominant} \end{equation}$ (C.3)As Eq. (C.3) corresponds to the 1st order expansion in i in Eq. (C.1), we needed to make sure that these terms did not cancel in the sum for this approximation to be correct. This sum is $\begin{matrix} \sum_{i = 0}^{N - 1} P & \approx & N - \frac{N (N - 1)}{2} \frac{4}{N} + \frac{(N - 1) N (2 N - 1)}{6} \frac{6}{N^{2}} \\ - \frac{(N - 1)^{2} N^{2}}{4} \frac{4}{N^{3}} i^{3} \\ + \frac{(N - 1) (N) (2 N - 1) (3 (N - 1)^{2} + 3 (N - 1) - 1)}{30 N^{4}} \cdot \\ \sum_{i = 0}^{N - 1} P & \approx & N - 2 N + 2 N - N + \frac{N}{5}, \\ so \end{matrix}$ $\appendix \setcounter{section}{3} \begin{eqnarray} \sum_{i=0}^{N-1} P &\approx& N - \dfrac{N(N-1)}{2} \dfrac{4}{N} + \dfrac{(N-1)N(2N-1)}{6}\dfrac{6}{N^2} \nonumber \\ &&\quad -\dfrac{(N-1)^2N^2}{4}\dfrac{4}{N^3}i^3 \nonumber \\ &&\quad+\dfrac{(N-1)(N)(2N-1)(3(N-1)^2+3(N-1)-1)}{30N^4}\cdot \nonumber \\ \sum_{i=0}^{N-1} P &\approx& N - 2N +2N - N + \dfrac{N}{5}, \nonumber \\ \!\!\!\!\text{ so }\; ~~~ &&\sum_{i=0}^{N-1}\left(1 - \left(i+\dfrac{1}{2} \right)/N\right)^{4} \sim \dfrac{N}{5}\cdot \label{P_fin} \end{eqnarray}$ (C.4)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, Eqs. (C.1), (C.2) and (C.4) yielded: $| Δ J_{2} | ~ \frac{π}{12} \frac{⟨ P_{2} ⟩}{N^{2}} \approx \frac{9.42 \times 10^{-3}}{N^{2}} \cdot$ $\appendix \setcounter{section}{3} \begin{equation} |\Delta J_{2}| \sim \dfrac{\pi}{12} \dfrac{\langle P_{2}\rangle }{N^2} \approx \dfrac{9.42\times 10^{-3}}{N^2}\cdot \end{equation}$ (C.5)

Appendix 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ⁱⁿ = 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. (B.21), we were able to write it with a general number N as: $\begin{matrix} | Δ J_{2 k} | & ~ & \frac{4 π (2 k + 3)}{12} \frac{A̅ ⟨ P_{2 k} ⟩}{M} \\ \times [\sum_{i = 0}^{\frac{2 N}{3}} {(\frac{3 Δ r}{4})}^{3} {(\frac{R_{i}^{ext}}{a_{J}})}^{2 k + 2} + \sum_{i = \frac{2 N}{3} + 1}^{N - 1} {(\frac{3 Δ r}{2})}^{3} {(\frac{R_{i}^{int}}{a_{J}})}^{2 k + 2}] \cdot \\ R_{i}^{ext} & = & a_{J} (1 - \frac{3}{4 N} (i + \frac{1}{2})), \\ R_{i}^{in} & = & a_{J} (0.5 - \frac{3}{2 N} ((i - \frac{2 N}{3}) + \frac{1}{2})) \cdot \end{matrix}$ $\appendix \setcounter{section}{3} \begin{eqnarray} |\Delta J_{2k}| &\sim& \dfrac{4\pi(2k+3)}{12}\dfrac{\bar{A}\langle P_{2k}\rangle}{M} \nonumber \\ &&\quad \times \left[\sum_{i=0}^{\frac{2N}{3}} \left(\dfrac{3 \Delta r}{4}\right)^3 \left( \dfrac{R_i^{\rm ext}}{a_J}\right)^{2k+2} + \sum_{i=\frac{2N}{3}+1}^{N-1} \left(\dfrac{3 \Delta r}{2}\right)^3 \left( \dfrac{R_i^{\rm int}}{a_J}\right)^{2k+2}\right] \cdot \\ R_i^{\rm ext} &=& a_J \left(1- \dfrac{3}{4N}\left(i+\dfrac{1}{2}\right)\right), ~~~\\ R_i^{\rm in} &= &a_J\left(0.5- \dfrac{3}{2N}\left(\left(i-\dfrac{2N}{3}\right)+\dfrac{1}{2}\right)\right)\cdot \end{eqnarray}$ Which can be rewritten: $\begin{matrix} | Δ J_{2 k} | & ~ & \frac{4 π (2 k + 3)}{12} \frac{A̅ ⟨ P_{2 k} ⟩}{M} \frac{27 Δ r^{3}}{8} \\ \times [\frac{1}{8} \sum_{i = 0}^{\frac{2 N}{3}} {(\frac{R_{i}^{ext}}{a_{J}})}^{2 k + 2} + \sum_{i = 0}^{\frac{N}{3}} {(\frac{R_{i - \frac{2 N}{3}}^{int}}{a_{J}})}^{2 k + 2}] \cdot \end{matrix}$ $\appendix \setcounter{section}{3} \begin{eqnarray} |\Delta J_{2k}| &\sim& \dfrac{4\pi(2k+3)}{12}\dfrac{\bar{A}\langle P_{2k}\rangle}{M} \dfrac{27 \Delta r^3}{8} \nonumber \\ &&\quad \times \left[\dfrac{1}{8}\sum_{i=0}^{\frac{2N}{3}} \left( \dfrac{R_i^{\rm ext}}{a_J}\right)^{2k+2} + \sum_{i=0}^{\frac{N}{3}} \left( \dfrac{R_{i-\frac{2N}{3}}^{\rm int}}{a_J}\right)^{2k+2}\right]\cdot \end{eqnarray}$ (C.9)For k = 1, using Eqs. (B.22) and(C.4), we obtained: $\begin{matrix} | Δ J_{2} | ~ \frac{5 \times 27 \times π ⟨ P_{2} ⟩}{12 \times 8 \times N^{3}} [\frac{1}{8} (\frac{4}{3} \frac{2 \times N}{3}) \frac{1}{5} + \frac{1}{2^{4}} (\frac{1}{3} \frac{N}{3}) \frac{1}{5}] \end{matrix}$ $\appendix \setcounter{section}{3} \begin{eqnarray} &&|\Delta J_{2}| \sim \dfrac{5\times 27\times \pi\langle P_{2}\rangle}{12\times 8\times N^3}\left[\dfrac{1}{8}\left(\dfrac{4}{3}\dfrac{2\times N}{3}\right) \dfrac{1}{5}+\dfrac{1}{2^4}\left(\dfrac{1}{3}\dfrac{N}{3}\right)\dfrac{1}{5}\right] \nonumber \\ &&|\Delta J_{2}| \sim \dfrac{9}{256} \dfrac{\pi \langle P_{2}\rangle }{N^2}\cdot \label{J2_hub_calc} \end{eqnarray}$ (C.10)Compared to Eq. (9), this is approximately twice better in terms of errors.

Appendix 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³/N³.

Using Eq. (B.21), we got: $\begin{matrix} R_{i} = a_{J} (1 - \frac{i^{3} + (i + 1)^{3}}{2 N^{3}}) \cdot \end{matrix}$ $\appendix \setcounter{section}{3} \begin{eqnarray} &&\Delta r_i = ((i+1)^3 -i^3) \dfrac{a_J}{N^3} = \dfrac{\Delta (i^3)}{N^2} \Delta r, \label{delta_ri} \\ &&R_i = a_J \left(1- \dfrac{i^3+(i+1)^3}{2N^3} \right)\cdot \end{eqnarray}$ Straightforwardly, we obtained the new Eq. (7): $\begin{matrix} | Δ J_{2 k}^{c} | & ≃ & (2 k + 3) \frac{π}{12} \frac{⟨ | P_{2 k} | ⟩}{N^{5}} \end{matrix}$ $\appendix \setcounter{section}{3} \begin{eqnarray} |\Delta J_{2k}^c| &\simeq& (2k+3)\dfrac{\pi}{12}\dfrac{\langle |P_{2k}|\rangle}{N^5} \nonumber \\ &&\quad \times \sum_{i=0}^{N-1} \Delta (i^3) \left(1- \dfrac{i^3+(i+1)^3}{2N^3} \right)^{2k+2}, \label{J_2k_c} \end{eqnarray}$ (C.13)which gave for J₂: $| Δ J_{2}^{c} | ≃ \frac{5 π ⟨ P_{2} ⟩}{12 N^{5}} \sum_{i = 0}^{N - 1} Δ (i^{3}) {(1 - \frac{i^{3} + (i + 1)^{3}}{2 N^{3}})}^{4} \cdot$ $\appendix \setcounter{section}{3} \begin{equation} |\Delta J_{2}^c| \simeq \dfrac{5 \pi \langle P_{2}\rangle }{12N^5} \sum_{i=0}^{N-1} \Delta (i^3) \left(1- \dfrac{i^3+(i+1)^3}{2N^3} \right)^{4}\cdot \end{equation}$ (C.14)We developed again the various terms: $\begin{matrix} Δ (i^{3}) = 3 i^{2} + 3 i + 1 ~ 3 i^{2}, \\ P^{c} = {(1 - \frac{i^{3} + (i + 1)^{3}}{2 N^{3}})}^{4}, \\ P^{c} = {(1 - \frac{1}{N^{3}} i^{3} - \frac{3}{2 N^{3}} i^{2} - \frac{3}{2 N^{3}} i - \frac{1}{2 N^{3}})}^{4} ≃ {(1 - \frac{i^{3}}{N^{3}})}^{4}, \\ P^{c} ≃ 1 - 4 \frac{i^{3}}{N^{3}} + 6 \frac{i^{6}}{N^{6}} - 4 \frac{i^{9}}{N^{9}} + \frac{i^{12}}{N^{12}}, \end{matrix}$ $\appendix \setcounter{section}{3} \begin{eqnarray} &&\Delta (i^3) = 3i^2 + 3i + 1 \sim 3i^2,\\ &&P^c = \left(1- \dfrac{i^3+(i+1)^3}{2N^3} \right)^{4}, \nonumber \\ &&P^c = \left(1- \dfrac{1}{N^3}i^3 - \dfrac{3}{2N^3}i^2 - \dfrac{3}{2N^3}i - \dfrac{1}{2N^3} \right)^4 \simeq \left(1-\dfrac{i^3}{N^3} \right)^4, \nonumber \\ && P^c \simeq 1 -4\dfrac{i^3}{N^3} + 6 \dfrac{i^6}{N^6} - 4 \dfrac{i^9}{N^9} + \dfrac{i^{12}}{N^{12}}, \\ &&\Delta (i^3) P^c \simeq 3 \left(i^2 -4\dfrac{i^5}{N^3} + 6\dfrac{i^8}{N^6} - 4\dfrac{i^{11}}{N^9} + \dfrac{i^{14}}{N^{12}} \right)\cdot \label{P_cubic} \end{eqnarray}$ In order to calculate the sum, we needed the leading term of $\sum_{i = 0}^{N - 1} i^{k}$ $\hbox{$\sum_{i=0}^{N-1} i^k$}$ . With: $\begin{matrix} (a + 1)^{k} - a^{k} = \sum_{i = 0}^{k - 1} (\begin{matrix} k \\ i \end{matrix}) a^{i}, \\ \sum_{a = 1}^{N - 1} ((a + 1)^{k} - a^{k}) = \sum_{i = 0}^{k - 1} (\begin{matrix} k \\ i \end{matrix}) \sum_{a = 1}^{N - 1} a^{i}, \end{matrix}$ $\appendix \setcounter{section}{3} \begin{eqnarray} &&(a+1)^k - a^k = \sum_{i=0}^{k-1} \dbinom{k}{i} a^i, \\ &&\sum_{a=1}^{N-1} \left((a+1)^k - a^k \right) = \sum_{i=0}^{k-1} \dbinom{k}{i} \sum_{a=1}^{N-1}a^i, \end{eqnarray}$ but the terms in the sum cancelled two by two: $\begin{matrix} \sum_{a = 1}^{N - 1} ((a + 1)^{k} - a^{k}) = N^{k} - 1, \\ so: N^{k} - 1 = \sum_{i = 0}^{k - 1} (\begin{matrix} k \\ i \end{matrix}) \sum_{a = 1}^{N - 1} a^{i} . \end{matrix}$ $\appendix \setcounter{section}{3} \begin{eqnarray} &&\sum_{a=1}^{N-1} \left((a+1)^k - a^k \right) = N^k - 1, \\ &&\text{ so: }\,\,\,\,N^k-1 = \sum_{i=0}^{k-1} \dbinom{k}{i} \sum_{a=1}^{N-1}a^i. \end{eqnarray}$ With this formula, it is easy to show by recurrence that $\sum_{a = 1}^{N - 1} a^{i} = O (N^{i + 1})$ $\hbox{$\sum_{a=1}^{N-1}a^i = O(N^{i+1})$}$ . Therefore, in the equation above only the k−1 term can lead to N^k. Then the result: $\sum_{a = 1}^{N - 1} a^{k - 1} ~ \frac{N^{k}}{(\begin{matrix} k \\ k - 1 \end{matrix})} = \frac{N^{k}}{k} \cdot$ $\appendix \setcounter{section}{3} \begin{equation} \sum_{a=1}^{N-1}a^{k-1} \sim \dfrac{N^k}{\dbinom{k}{k-1}} = \dfrac{N^k}{k}\cdot \end{equation}$ (C.22)From here, we obtained the leading order terms of the sum of Eq. (C.17): $\begin{matrix} \sum_{i = 0}^{N - 1} Δ (i^{3}) P^{c} ~ 3 (\frac{N^{3}}{3} - \frac{4 N^{3}}{6} + \frac{6 N^{3}}{9} - \frac{4 N^{3}}{12} + \frac{N^{3}}{15}) = \frac{N^{3}}{5} \cdot \end{matrix}$ $\appendix \setcounter{section}{3} \begin{eqnarray} \sum_{i=0}^{N-1}\Delta (i^3) P^c \sim 3 \left(\dfrac{N^3}{3} - \dfrac{4N^3}{6} + \dfrac{6N^3}{9} - \dfrac{4N^3}{12} + \dfrac{N^3}{15} \right) = \dfrac{N^3}{5}\cdot \end{eqnarray}$ (C.23)The result is thus the same as in the linear case: $| Δ J_{2}^{c} | ~ \frac{π}{12} \frac{⟨ P_{2} ⟩}{N^{2}} \cdot$ $\appendix \setcounter{section}{3} \begin{equation} |\Delta J_{2}^c| \sim \dfrac{\pi}{12} \dfrac{\langle P_{2}\rangle }{N^2} \cdot \end{equation}$ (C.24)

Appendix C.4: Exponential spacing

For an exponential repartition of the spheroids, we imposed: $λ_{i + 1} = λ_{i} - β e^{iα} .$ $\appendix \setcounter{section}{3} \begin{equation} \lambda_{i+1} = \lambda_{i} - \beta {\rm e}^{i\alpha}. \end{equation}$ (C.25)Then, considering the difference λ_{i + 1}−λ_i as a geometric sequence, we obtained: $λ_{i} = 1 - β \frac{e^{iα} - 1}{e^{α} - 1} \cdot$ $\appendix \setcounter{section}{3} \begin{equation} \lambda_{i} = 1 - \beta \dfrac{{\rm e}^{i\alpha}-1}{{\rm e}^{\alpha}-1}\cdot \end{equation}$ (C.26)Now, with the condition that the first layer (normalized to the radius of Jupiter) is 1 and the last one is 0: $λ_{i} = 1 - \frac{e^{iα} - 1}{e^{Nα} - 1} \cdot$ $\appendix \setcounter{section}{3} \begin{equation} \lambda_{i} = 1 - \dfrac{{\rm e}^{i\alpha}-1}{{\rm e}^{N\alpha}-1}\cdot \label{lambda_exp} \end{equation}$ (C.27)To determine the error we needed: $\begin{matrix} R_{i} & = & 1 + \frac{1}{e^{Nα} - 1} - \frac{e^{iα}}{2 (e^{Nα} - 1)} - \frac{e^{(i + 1) α}}{2 (e^{Nα} - 1)} \\ = & 1 + \frac{1}{e^{Nα} - 1} - e^{iα} \frac{1 + e^{α}}{2 (e^{Nα} - 1)}, \\ Δ r & = & \frac{e^{(i + 1) α} - e^{iα}}{e^{Nα} - 1} = e^{iα} \frac{e^{α} - 1}{e^{Nα} - 1} = β e^{iα} . \end{matrix}$ $\appendix \setcounter{section}{3} \begin{eqnarray} R_{i} &=& 1 + \dfrac{1}{{\rm e}^{N\alpha}-1} - \dfrac{{\rm e}^{i\alpha}}{2({\rm e}^{N\alpha}-1)} -\dfrac{{\rm e}^{(i+1)\alpha}}{2({\rm e}^{N\alpha}-1)} \nonumber \\ &=&1 + \dfrac{1}{{\rm e}^{N\alpha}-1} - {\rm e}^{i\alpha}\dfrac{1+{\rm e}^{\alpha}}{2({\rm e}^{N\alpha}-1)}, \\ \Delta r &=& \dfrac{{\rm e}^{(i+1)\alpha}-{\rm e}^{i\alpha}}{{\rm e}^{N\alpha}-1} = {\rm e}^{i\alpha}\dfrac{{\rm e}^{\alpha}-1}{{\rm e}^{N\alpha}-1} = \beta {\rm e}^{i\alpha}. \end{eqnarray}$ As we wanted the increment to be reasonable, e^Nα−1 could not be too large. In this paper, we chose $α = \frac{γ}{N},$ $\appendix \setcounter{section}{3} \begin{equation} \alpha = \dfrac{\gamma}{N}, \end{equation}$ (C.30)where we varied γ between 5 and 10, depending on how sharp we wanted the exponential function to be (see Fig. 2).

Then, using Eq. (B.21), we had to calculate, for J₂: $\begin{matrix} \sum_{i = 0}^{N - 1} (Δ r_{i})^{3} R_{i}^{4} = β^{3} \sum_{i = 0}^{N - 1} e^{3 iα} {([1 + \frac{1}{e^{Nα} - 1}] - e^{iα} \frac{1 + e^{α}}{2 (e^{Nα} - 1)})}^{4}, \\ \sum_{i = 0}^{N - 1} (Δ r_{i})^{3} R_{i}^{4} = β^{3} \sum_{i = 0}^{N - 1} e^{3 iα} (δ - ϵ e^{iα})^{4}, \\ β = \frac{e^{α} - 1}{e^{Nα} - 1}, δ = 1 + \frac{1}{e^{Nα} - 1}, ϵ = \frac{1 + e^{α}}{2 (e^{Nα} - 1)} \cdot \end{matrix}$ $\appendix \setcounter{section}{3} \begin{eqnarray} &&\sum_{i=0}^{N-1} (\Delta r_i)^3 R_i^4 = \beta^3 \sum_{i=0}^{N-1} {\rm e}^{3i\alpha} \left( \left[1 + \dfrac{1}{{\rm e}^{N\alpha}-1} \right] - {\rm e}^{i\alpha}\dfrac{1+{\rm e}^{\alpha}}{2({\rm e}^{N\alpha}-1)} \right)^4, \nonumber \\ &&\sum_{i=0}^{N-1} (\Delta r_i)^3 R_i^4 = \beta^3 \sum_{i=0}^{N-1} {\rm e}^{3i\alpha}(\delta - \epsilon {\rm e}^{i\alpha})^4, \\ &&\beta = \dfrac{{\rm e}^{\alpha}-1}{{\rm e}^{N\alpha}-1} \text{ , } \delta = 1 + \dfrac{1}{{\rm e}^{N\alpha}-1} \text{ , } \epsilon = \dfrac{1+{\rm e}^{\alpha}}{2({\rm e}^{N\alpha}-1)}\cdot \end{eqnarray}$ Developing the power of 4 brackets and using the well known result for the sum of the terms of a geometric sequence: $\begin{matrix} \sum_{i = 0}^{N - 1} (Δ r_{i})^{3} R_{i}^{4} & = & β^{3} δ^{4} (\frac{e^{3 Nα} - 1}{e^{3 α} - 1} - 4 \frac{ϵ}{δ} \frac{e^{4 Nα} - 1}{e^{4 α} - 1} \\ + 6 {(\frac{ϵ}{δ})}^{2} \frac{e^{5 Nα} - 1}{e^{5 α} - 1} - 4 {(\frac{ϵ}{δ})}^{3} {\frac{e^{6 Nα} - 1}{e^{6 α} - 1}}^{} \\ + {(\frac{ϵ}{δ})}^{4} \frac{e^{7 Nα} - 1}{e^{7 α} - 1}) \cdot \end{matrix}$ $\appendix \setcounter{section}{3} \begin{eqnarray} \sum_{i=0}^{N-1} (\Delta r_i)^3 R_i^4 &=& \beta^3 \delta^4 \left(\dfrac{{\rm e}^{3N\alpha} -1}{{\rm e}^{3\alpha}-1} - 4 \dfrac{\epsilon}{\delta} \dfrac{{\rm e}^{4N\alpha} -1}{{\rm e}^{4\alpha}-1} \nonumber \right. \\&&\quad \left. + \,6 \left(\dfrac{\epsilon}{\delta}\right)^2\dfrac{{\rm e}^{5N\alpha} -1} {{\rm e}^{5\alpha}-1} - 4 \left(\dfrac{\epsilon}{\delta}\right)^3\dfrac{{\rm e}^{6N\alpha} -1} {{\rm e}^{6\alpha}-1} \right. \nonumber \\ &&\quad +\left. \left(\dfrac{\epsilon}{\delta}\right)^4\dfrac{{\rm e}^{7N\alpha} -1} {{\rm e}^{7\alpha}-1} \right)\cdot \end{eqnarray}$ (C.33)Here, we considered that 7α ≪ 1, e^Nα ≫ 1, which is a reasonable approximation in the range γ = 5−10 and N> 512. Then: $\begin{matrix} β ~ \frac{α}{e^{Nα}}, δ ~ 1, ϵ ~ \frac{1}{e^{Nα}}, \frac{e^{3 Nα} - 1}{e^{3 α} - 1} ~ \frac{e^{3 Nα}}{3 α}, \\ \sum_{i = 0}^{N - 1} (Δ r_{i})^{3} R_{i}^{4} ~ {(\frac{α}{e^{Nα}})}^{3} (\frac{e^{3 Nα}}{3 α} - 4 \frac{e^{3 Nα}}{4 α} \\ + 6 \frac{e^{3 Nα}}{5 α} - 4 \frac{e^{3 Nα}}{6 α} + \frac{e^{3 Nα}}{7 α}), \\ \sum_{i = 0}^{N - 1} (Δ r_{i})^{3} R_{i}^{4} ~ - \frac{α^{2}}{105} = - \frac{γ^{2}}{105 N^{2}} \cdot \end{matrix}$ $\appendix \setcounter{section}{3} \begin{eqnarray} &&\beta \sim \dfrac{\alpha}{{\rm e}^{N\alpha}} \text{ , } \delta \sim 1 \text{ , } \epsilon \sim \dfrac{1}{{\rm e}^{N\alpha}} \text{ , } \dfrac{{\rm e}^{3N\alpha} -1}{{\rm e}^{3\alpha}-1} \sim \dfrac{{\rm e}^{3N\alpha}}{3\alpha}, \\ &&\sum_{i=0}^{N-1} (\Delta r_i)^3 R_i^4 \sim \left(\dfrac{\alpha}{{\rm e}^{N\alpha}} \right)^3 \left(\dfrac{{\rm e}^{3N\alpha}}{3\alpha} - 4 \dfrac{{\rm e}^{3N\alpha}}{4\alpha} \right. \nonumber \\ &&\phantom{\sum_{i=0}^{N-1} (\Delta r_i)^3 R_i^4 \sim } +\left. 6\dfrac{{\rm e}^{3N\alpha}}{5\alpha} - 4 \dfrac{{\rm e}^{3N\alpha}}{6\alpha} + \dfrac{{\rm e}^{3N\alpha}}{7\alpha} \right), \nonumber \\ &&\sum_{i=0}^{N-1} (\Delta r_i)^3 R_i^4 \sim -\dfrac{\alpha^2}{105} = -\dfrac{\gamma^2}{105 \, N^2}\cdot \end{eqnarray}$ We checked numerically that our approximations yield an error always smaller than ≈10%.

Implementing these results into Eq. (B.21) yields: $| Δ J_{2}^{\exp} | ~ \frac{π ⟨ P_{2} ⟩}{252} * \frac{γ^{2}}{N^{2}} \approx 0.000449 {(\frac{γ}{N})}^{2},$ $\appendix \setcounter{section}{3} \begin{equation} |\Delta J_{2}^{\exp}| \sim \dfrac{\pi \langle P_{2}\rangle }{252}*\dfrac{\gamma^2}{ N^2} \approx 0.000449 \left(\dfrac{\gamma}{N}\right)^2, \end{equation}$ (C.36)with γ ∈ [5−10].

Appendix D: Neglecting the high atmosphere

The neglected potential on a point (μ, r_j) of the jth spheroid reads: $\begin{matrix} φ_{neg} = 4 πG \sum_{k = 0}^{\infty} (r_{j})^{2 k} P_{2 k} (μ) \int_{a_{1 bar}}^{a_{ext}} \int_{0}^{1} \frac{ρ (r^{'})}{r^{' 2 k - 1}} P_{2 k} (μ^{'}) d μ^{'} d r^{'} . \end{matrix}$ $\appendix \setcounter{section}{4} \begin{eqnarray} \phi_{\rm neg} = 4 \pi G \sum_{k=0}^{\infty} (r_j)^{2k} P_{2k}(\mu) \int_{a_{\rm 1\,bar}}^{a_{\rm ext}}\int_{0}^{1} \dfrac{\rho(\vec{r'})}{r'^{2k-1}} P_{2k}(\mu'){\rm d}\mu'\, {\rm d}r'. \end{eqnarray}$ (D.1)With ρ = ρ_ext = const. and Appendix A: $\begin{matrix} I & = & \int \begin{matrix} a_{ext} \\ a_{1 bar} \end{matrix} \int_{0}^{1} ρ (r^{'}) \frac{1}{r^{' 2 k - 1}} P_{2 k} (μ^{'}) d μ^{'} d r^{'} \\ ≃ & ρ_{ext} \int_{a_{1 bar}}^{a_{ext}} \frac{1}{r^{' 2 k - 1}} d r^{'} \int_{0}^{1} P_{2 k} (μ^{'}) (1 + e^{2} μ^{' 2})^{k - 1} d μ^{'} . \end{matrix}$ $\appendix \setcounter{section}{4} \begin{eqnarray} I &=& \int_{a_{\rm 1\,bar}}^{a_{\rm ext}}\int_{0}^{1} \rho(\vec{r'})\dfrac{1}{r'^{2k-1}} P_{2k}(\mu'){\rm d}\mu'\, {\rm d}r' \nonumber \\ &\simeq& \rho_{\rm ext}\int_{a_{\rm 1\,bar}}^{a_{\rm ext}}\dfrac{1}{r'^{2k-1}}{\rm d}r' \int_{0}^{1} P_{2k}(\mu') (1+{\rm e}^2 \mu'^2)^{k-1}{\rm d}\mu'. \end{eqnarray}$ (D.2)For any k> 0, we integrate a Legendre polynomial of degree 2k with a polynomial of order 2(k−1). By orthogonality of the Legendre base, only the k = 0 term is non zero. Thus $\begin{matrix} φ_{neg} ≃ 4 πG ρ_{ext} [\frac{1}{2} (a_{ext}^{2} - a_{1 bar}^{2}) \frac{\arctan (e)}{e}] \cdot \end{matrix}$ $\appendix \setcounter{section}{4} \begin{eqnarray} \phi_{\rm neg} \simeq 4 \pi G \rho_{\rm ext} \left[\dfrac{1}{2} (a_{\rm ext}^2 - a_{\rm 1\,bar}^2) \dfrac{\arctan(e)}{e} \right]\cdot \end{eqnarray}$ (D.3)With a_ext = a_{1 bar} + Δa: $\begin{matrix} (a_{ext}^{2} - a_{1 bar}^{2}) ≃ 2 a_{1 bar} Δ a,andso \\ φ_{neg} ≃ 4 πG ρ_{ext} a_{1 bar} Δ a \frac{\arctan (e)}{e} \cdot \end{matrix}$ $\appendix \setcounter{section}{4} \begin{eqnarray} &&(a_{\rm ext}^2 - a_{\rm 1\,bar}^2) \simeq 2 a_{\rm 1\,bar}\Delta a \text{, and so}\\ &&\phi_{\rm neg} \simeq 4 \pi G \rho_{\rm ext} a_{\rm 1\,bar} \Delta a \dfrac{\arctan(e)}{e}\cdot \end{eqnarray}$ Using the exterior e ≃ 0.38 and defining $\bar{ρ_{S}}$ $\hbox{$\bar{\rho_{\rm S}}$}$ as the density of Jupiter if it was a sphere of radius a_{1 bar}, the neglected potential simply reads: $\begin{matrix} φ_{neg} ≃ 2.87 \frac{GM}{a_{1 bar}} (\frac{ρ_{ext}}{\bar{ρ_{S}}}) (\frac{Δ a}{a_{1 bar}}) \cdot \end{matrix}$ $\appendix \setcounter{section}{4} \begin{eqnarray} \phi_{\rm neg} \simeq 2.87 \dfrac{G M}{a_{\rm 1\,bar}} \left(\dfrac{\rho_{\rm ext}}{\bar{\rho_{\rm S}}} \right) \left(\dfrac{\Delta a}{a_{\rm 1\,bar}} \right) \cdot \end{eqnarray}$ (D.6)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_{1 bar}/GM and U_j the dimensionless potential of the jth spheroid, we derived the first order perturbation of Eq. (51) of H13: $\begin{matrix} - \frac{1}{ζ_{j} (μ)} & (\sum_{i = j}^{N - 1} \sum_{k = 0}^{\infty} \begin{matrix} ˜ \\ J_{i, 2 k} \end{matrix} {(\frac{λ_{i}}{λ_{j}})}^{2 k} ζ_{j} (μ)^{- 2 k} P_{2 k} (μ) \\ + \sum_{i = 0}^{j - 1} \sum_{k = 0}^{\infty} {\begin{matrix} ˜ \\ J^{'} \end{matrix}}_{i, 2 k} {(\frac{λ_{j}}{λ_{i}})}^{2 k + 1} ζ_{j} (μ)^{2 k + 1} P_{2 k} (μ) \end{matrix}$ $\appendix \setcounter{section}{4} \begin{eqnarray} -\dfrac{1}{\zeta_j (\mu)} &&\left( \sum_{i=j}^{N-1} \sum_{k=0}^{\infty} \tilde{J}_{i,2k} \left(\dfrac{\lambda_i}{\lambda_j}\right)^{2k} \zeta_j(\mu)^{-2k} P_{2k} (\mu) \right. \nonumber \\ &&+\sum_{i=0}^{j-1} \sum_{k=0}^{\infty} \tilde{J'}_{i,2k} \left(\dfrac{\lambda_j}{\lambda_i}\right)^{2k+1} \zeta_j(\mu)^{2k+1} P_{2k} (\mu) \nonumber \\ &&\left. + \sum_{i=0}^{j-1} \tilde{J''}_{i,0} \lambda_j^3 \zeta_j(\mu)^3 \right) + \dfrac{1}{2} \lambda_j^3 \zeta_j(\mu)^2 (1 - \mu^2) = U_{j} + \Delta U. \label{eq:51_HM13} \end{eqnarray}$ (D.7)For simplicity, we set $\hbox{$x = \bar{\zeta}_j (\mu)$}$ 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. (D.7)at the first order: $\begin{matrix} - \frac{1}{x} (- \sum_{i = j}^{N - 1} \sum_{k = 0}^{\infty} \begin{matrix} ˜ \\ J_{i, 2 k} \end{matrix} {(\frac{λ_{i}}{λ_{j}})}^{2 k} x^{- 2 k} (2 k \frac{Δ x}{x}) P_{2 k} (μ) \\ + \sum_{i = 0}^{j - 1} \sum_{k = 0}^{\infty} {\begin{matrix} ˜ \\ J^{'} \end{matrix}}_{i, 2 k} {(\frac{λ_{j}}{λ_{i}})}^{2 k + 1} x^{2 k + 1} ((2 k + 1) \frac{Δ x}{x}) P_{2 k} (μ) \\ + \sum_{i = 0}^{j - 1} {\begin{matrix} ˜ \\ J^{′′} \end{matrix}}_{i, 0} λ_{j}^{3} x^{3} (3 \frac{Δ x}{x})) \\ + \frac{1}{x} \frac{Δ x}{x} (\sum_{i = j}^{N - 1} \sum_{k = 0}^{\infty} \begin{matrix} ˜ \\ J_{i, 2 k} \end{matrix} {(\frac{λ_{i}}{λ_{j}})}^{2 k} x^{- 2 k} P_{2 k} (μ) \\ + \sum_{i = 0}^{j - 1} \sum_{k = 0}^{\infty} {\begin{matrix} ˜ \\ J^{'} \end{matrix}}_{i, 2 k} {(\frac{λ_{j}}{λ_{i}})}^{2 k + 1} x^{2 k + 1} P_{2 k} (μ) \\ + \sum_{i = 0}^{j - 1} {\begin{matrix} ˜ \\ J^{′′} \end{matrix}}_{i, 0} λ_{j}^{3} x^{3}) + λ_{j}^{3} x^{2} (\frac{Δ x}{x}) (1 - μ^{2}) = Δ U . \end{matrix}$ $\appendix \setcounter{section}{4} \begin{eqnarray} &&-\dfrac{1}{x} \Bigg (- \sum_{i=j}^{N-1} \sum_{k=0}^{\infty} \tilde{J}_{i,2k} \left(\dfrac{\lambda_i}{\lambda_j}\right)^{2k} x^{-2k} \left(2k \dfrac{\Delta x}{x} \right) P_{2k} (\mu) \nonumber \\ &&\qquad +\sum_{i=0}^{j-1} \sum_{k=0}^{\infty} \tilde{J'}_{i,2k} \left(\dfrac{\lambda_j}{\lambda_i}\right)^{2k+1} x^{2k+1} \left((2k+1) \dfrac{\Delta x}{x} \right) P_{2k} (\mu) \nonumber \\ &&\qquad+ \sum_{i=0}^{j-1} \tilde{J''}_{i,0} \lambda_j^3 x^3 \left(3\dfrac{\Delta x}{x} \right) \Bigg) \nonumber \\ &&\qquad +\dfrac{1}{x} \dfrac{\Delta x}{x} \Bigg ( \sum_{i=j}^{N-1} \sum_{k=0}^{\infty} \tilde{J}_{i,2k} \left(\dfrac{\lambda_i}{\lambda_j}\right)^{2k} x^{-2k} P_{2k} (\mu) \nonumber \\ &&\qquad+ \sum_{i=0}^{j-1} \sum_{k=0}^{\infty} \tilde{J'}_{i,2k} \left(\dfrac{\lambda_j}{\lambda_i}\right)^{2k+1} x^{2k+1} P_{2k} (\mu) \nonumber \\ &&\qquad + \sum_{i=0}^{j-1} \tilde{J''}_{i,0} \lambda_j^3 x^3 \Bigg) + \lambda_j^3 x^2 \left(\dfrac{\Delta x}{x} \right) (1 - \mu^2) = \Delta U. \end{eqnarray}$ (D.8)The second bracket could be simplified by the fact that x satisifies Eq. (51) of H13: $\begin{matrix} - \frac{Δ x}{x} (- \sum_{i = j}^{N - 1} \sum_{k = 0}^{\infty} 2 k \begin{matrix} ˜ \\ J_{i, 2 k} \end{matrix} {(\frac{λ_{i}}{λ_{j}})}^{2 k} x^{- 2 k - 1} P_{2 k} (μ) \\ + \sum_{i = 0}^{j - 1} \sum_{k = 0}^{\infty} (2 k + 1) {\begin{matrix} ˜ \\ J^{'} \end{matrix}}_{i, 2 k} {(\frac{λ_{j}}{λ_{i}})}^{2 k + 1} x^{2 k} P_{2 k} (μ) \\ + \sum_{i = 0}^{j - 1} 3 {\begin{matrix} ˜ \\ J^{′′} \end{matrix}}_{i, 0} λ_{j}^{3} x^{2}) - \frac{Δ x}{x} (U_{j} - \frac{1}{2} λ_{j}^{3} x^{2} (1 - μ^{2})) \\ + (\frac{Δ x}{x}) λ_{j}^{3} x^{2} (1 - μ^{2}) = Δ U . \end{matrix}$ $\appendix \setcounter{section}{4} \begin{eqnarray} &&-\dfrac{\Delta x}{x} \Bigg (- \sum_{i=j}^{N-1} \sum_{k=0}^{\infty} 2k \tilde{J}_{i,2k} \left(\dfrac{\lambda_i}{\lambda_j}\right)^{2k} x^{-2k-1} P_{2k} (\mu) \nonumber \\ &&\qquad +\sum_{i=0}^{j-1} \sum_{k=0}^{\infty}(2k+1) \tilde{J'}_{i,2k} \left(\dfrac{\lambda_j}{\lambda_i}\right)^{2k+1} x^{2k} P_{2k} (\mu) \nonumber \\ &&\qquad + \sum_{i=0}^{j-1} 3\tilde{J''}_{i,0} \lambda_j^3 x^2 \Bigg) -\dfrac{\Delta x}{x} \Bigg (U_j - \dfrac{1}{2} \lambda_j^3 x^2 (1 - \mu^2)\Bigg) \nonumber \\ &&\qquad+ \left(\dfrac{\Delta x}{x} \right)\lambda_j^3 x^2 (1 - \mu^2) = \Delta U. \end{eqnarray}$ (D.9)Now, we restricted ourselves to the outermost spheroid j = 0, λ_j = 1. Only the first sum remained, and we recognized Eq. (1): $\begin{matrix} \frac{Δ x}{x} (\sum_{k = 0}^{\infty} 2 k J_{2 k}^{ext} x^{- 2 k - 1} P_{2 k} (μ) - U_{j} + \frac{3}{2} x^{2} (1 - μ^{2})) = Δ U . \end{matrix}$ $\appendix \setcounter{section}{4} \begin{eqnarray} \dfrac{\Delta x}{x} \left(\sum_{k=0}^{\infty} 2k J_{2k}^{\rm ext} x^{-2k-1} P_{2k} (\mu) - U_j + \dfrac{3}{2}x^2\left(1 - \mu^2\right) \right)= \Delta U. \end{eqnarray}$ (D.10)

With Juno data, we know the first few $J_{2 k}^{ext}$ $\hbox{$J_{2k}^{\rm ext}$}$ 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 Eq. (40) of H13: $\begin{matrix} \begin{matrix} ˜ \\ J_{i, 2 k} \end{matrix} = - \frac{3}{2 k + 3} \frac{δ ρ_{i} λ_{i}^{3} \int_{0}^{1} P_{2 k} (μ) ζ_{i} (μ)^{2 k + 3} d μ}{\sum_{j = 0}^{N - 1} δ ρ_{j} λ_{j}^{3} \int_{0}^{1} ζ_{j} (μ)^{3} d μ} \cdot \end{matrix}$ $\appendix \setcounter{section}{4} \begin{eqnarray} \tilde{J}_{i,2k} = -\dfrac{3}{2k+3} \dfrac{\delta \rho_i \lambda_i^3 \int_0^1 P_{2k}(\mu) \zeta_i(\mu)^{2k+3}{\rm d}\mu}{\sum_{j=0}^{N-1} \delta \rho_j \lambda_j^3 \int_0^1 \zeta_j(\mu)^3 {\rm d}\mu}\cdot \label{eq:delta_atmo} \end{eqnarray}$ (D.11)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₂. By comparing for different number of points in the Riemann sum, we obtained: $\frac{Δ J_{2, 0}}{J_{2, 0}} ≃ 100 \times Δ U ≃ 287 \times (\frac{ρ_{ext}}{\bar{ρ_{S}}}) (\frac{Δ a}{a_{1 bar}}) \cdot$ $\appendix \setcounter{section}{4} \begin{equation} \dfrac{\Delta J_{2,0}}{J_{2,0}} \simeq 100 \times \Delta U \simeq 287 \times \left(\dfrac{\rho_{\rm ext}}{\bar{\rho_{\rm S}}} \right) \left(\dfrac{\Delta a}{a_{\rm 1\,bar}} \right)\cdot \end{equation}$ (D.12)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: $\begin{matrix} \frac{Δ J_{4, 0}}{J_{4, 0}} ≃ 1 \times 10^{3} \times Δ U \\ \frac{Δ J_{6, 0}}{J_{6, 0}} ≃ 6 \times 10^{3} \times Δ U \\ \frac{Δ J_{8, 0}}{J_{8, 0}} ≃ 1 \times 10^{5} \times Δ U . \end{matrix}$ $\appendix \setcounter{section}{4} \begin{eqnarray} &&\dfrac{\Delta J_{4,0}}{J_{4,0}} \simeq 1 \times 10^3 \times \Delta U \\ &&\dfrac{\Delta J_{6,0}}{J_{6,0}} \simeq 6 \times 10^3\times \Delta U \\ &&\dfrac{\Delta J_{8,0}}{J_{8,0}} \simeq 1 \times 10^5 \times \Delta U . \end{eqnarray}$

All Tables

Table 1

Values of the planetary parameters of Jupiter.

In the text

Table 2

Difference between the values calculated with our code and the results from Hubbard (2013).

In the text

Table 3

Value of J₂ × 10⁶ for an exponential γ = 7 spacing with 1500 spheroids and the SCvH eos with no core.

In the text

Table 4

Difference between square, exponential γ = 7 and HM16 spacing for 3 different tests with 1500 spheroids.

In the text

Table 5

Effective helium mass fraction above (Y′) and underneath ( $Y_{2}^{'}$ $\hbox{$Y^\prime_2$}$ ) the entropy jump, and mass of the core (M_c in units of 10²⁴ kg) for different repartitions and numbers of spheroids.

In the text

Table 6

Same as Table 5 when neglecting the high atmosphere.

In the text

Table 7

Same as Table 6 except that J₂ × 10⁶ and J₄ × 10⁶ have been changed by ± 1 × 10^-1.

In the text

All Figures

	Fig. 1 Contribution of each part of the planet to the total value of J₂, normalized to its maximum value, for a polytrope of index n = 1.
In the text

	Fig. 2 Equatorial radius as a function of the spheroid number for various functional forms for the repartition of the spheroids. a) from 0 to a_J; b) zoom of the external layers.
In the text

Fig. 3

Absolute value of the difference between the expected and numerical results on J₂ × 10⁶ 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.

In the text

	Fig. 4 Numerical error as a function of the number of spheroids for the cubic and square repartitions. a) Error between the analytical and measured J₂ × 10⁶; b) difference on J₂ × 10⁶ when the first layer has 0 or non-0 density.
In the text

	Fig. 5 Numerical error as a function of the number of spheroids for exponential repartitions. a) Error between the analytical and measured J₂ × 10⁶; b) difference on the J₂ × 10⁶ when the first layer has 0 or non-0 density.
In the text

	Fig. 6 Normalized value of $δ ρ_{i} λ_{i}^{5}$ $\hbox{$\delta \rho_i \lambda_i^5$}$ as a function of altitude for linear, cubic and exponential (γ = 8) repartitions of spheroids. a) From center to exterior; b) zoom on the external layers.
In the text

	Fig. 7 Contribution of each spheroid to the total value of J₂ × 10⁶ and to the truncated J₂ × 10⁶, 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.
In the text

	Fig. 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₄ × 10⁶ vs. J₂ × 10⁶ from test 1; b) J₆ × 10⁶ vs. −J₄ × 10⁶ from test 3, when J₂ is fit.
In the text

	Fig. 9 Comparison of the errors of the CMS arising from neglecting the high atmosphere with Juno’s error bar, shifted to be centered on the observed J_n. a) −J₄ × 10⁶ vs. J₂ × 10⁶; b) J₆ × 10⁶ vs. −J₄ × 10⁶.
In the text

	Fig. 10 Numerical error as a function of the number of spheroids when neglecting the high atmosphere for a polytropic eos.
In the text

	Fig. B.1 Normalized value of C2 + 4C1 /R_i and of its approximation $5 R_{i} / a_{J}^{3}$ $\hbox{$5R_i/a_J^3$}$ (Eq. (B.20)) as a function of altitude R_i. a) From center to exterior; b) zoom on the external layers
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.

[1] Archinal, B. A., A’Hearn, M. F., Bowell, E., et al. 2011, Celes. Mech. Dyn. Astron., 109, 101 [Google Scholar]

[2] Bolton, S. J., Adriani, A., Adumitroaie, V., et al. 2017, Science, 356, 821 [NASA ADS] [CrossRef] [Google Scholar]

[3] Chabrier, G., Saumon, D., Hubbard, W. B., & Lunine, J. I. 1992, ApJ, 391, 817 [NASA ADS] [CrossRef] [Google Scholar]

[4] Cohen, E. R., & Taylor, B. N. 1987, Rev. Mod. Phys., 59, 1121 [NASA ADS] [CrossRef] [Google Scholar]

[5] Folkner, W. M., Iess, L., Anderson, J. D., et al. 2017, Geophys. Res. Lett., 44, 4694 [NASA ADS] [CrossRef] [Google Scholar]

[6] Hubbard, W. B. 1975, Sov. Astron., 18, 621 [NASA ADS] [Google Scholar]

[7] Hubbard, W. B. 2012, ApJ, 756, L15 [NASA ADS] [CrossRef] [Google Scholar]

[8] Hubbard, W. B. 2013, ApJ, 768, 43 [NASA ADS] [CrossRef] [Google Scholar]

[9] Hubbard, W. B., & Militzer, B. 2016, ApJ, 820, 80 [NASA ADS] [CrossRef] [Google Scholar]

[10] Hubbard, W. B., Schubert, G., Kong, D., & Zhang, K. 2014, Icarus, 242, 138 [NASA ADS] [CrossRef] [Google Scholar]

[11] Nettelmann, N. 2017, A&A, 606, A139 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[12] Riddle, A. C., & Warwick, J. W. 1976, Icarus, 27, 457 [NASA ADS] [CrossRef] [Google Scholar]

[13] Saumon, D., Chabrier, G., & van Horn, H. M. 1995, ApJS, 99, 713 [NASA ADS] [CrossRef] [Google Scholar]

[14] Wahl, S. M., Hubbard, W. B., Militzer, B., et al. 2017, Geophys. Res. Lett., 44, 4649 [NASA ADS] [CrossRef] [Google Scholar]

[15] Wisdom, J. 1996, http://web.mit.edu/wisdom/www/interior.pdf [Google Scholar]

[16] Wisdom, J., & Hubbard, W. B. 2016, Icarus, 267, 315 [NASA ADS] [CrossRef] [Google Scholar]

[17] Zharkov, V. N., & Trubitsyn, V. P. 1978, The physics of planetary interiors (Tuscon: Pachart Publishing House) [Google Scholar]

A complete study of the precision of the concentric MacLaurin spheroid method to calculate Jupiter’s gravitational moments

1. Introduction

2. Analytical evaluation of the errors in the CMS method

2.1. Evaluation of the errors with 512 spheroids

2.2. Spacing as a power of k

2.3. Exponential spacing of the layers

3. Numerical calculations

3.1. How to match Juno’s error bars

3.2. Importance of the first layers

3.3. One bar radius and external radius

4. Calculations with a realistic equation of state

4.1. Impact of the high atmospheric layers on the CMS method

4.2. Errors arising from a (quasi) linear repartition

4.3. Intrinsic uncertainties on the Jupiter models

5. Taking the 1 bar level as the outer boundary condition

5.1. Irreducible errors due to the high atmosphere region (<1 bar)

5.2. Error from the finite number of spheroids

6. Conclusion

Acknowledgments

References

Appendix A: Simplifications for the angular part

Appendix B: Supplementary calculations for the general case

Appendix C: Supplementary calculations for J2

Appendix C.1: Linear spacing

Appendix C.2: Hubbard and Militzer spacing

Appendix C.3: Cubic spacing

Appendix C.4: Exponential spacing

Appendix D: Neglecting the high atmosphere

All Tables

All Figures

Appendix C: Supplementary calculations for J₂