Linking Zonal Winds and Gravity: The Relative Importance of Dynamic Self Gravity

Recent precise measurements at Jupiter’s and Saturn’s gravity fields constrain the properties of the zonal flows in the outer envelopes of these planets. A simplified dynamic equation, sometimes called the thermal wind or thermo-gravitational wind equation, establishes a link between zonal flows and the related buoyancy perturbation, which in turn can be exploited to yield the dynamic gravity perturbation. Whether or not the action of the dynamic gravity perturbation needs to be explicitly included in this equation, an effect we call the Dynamic Self Gravity (DSG), has been a matter of intense debate. We show that, under reasonable assumptions, the equation can be solved (semi) analytically. This allows us to quantify the impact of the DSG on each gravity harmonic, practically independent of the zonal flow or the details of the planetary interior model. The impact decreases with growing spherical harmonic degree `. For degrees ` = 2 to about ` = 4, the DSG is a first order effect and should be taken into account in any attempt of inverting gravity measurements for zonal flow properties. For degrees of about ` = 5 to roughly ` = 10, the relative impact of DSG is about 10% and thus seems worthwhile to include, in particular since this comes at little extra costs with the method presented here. For yet higher degrees, is seems questionable whether gravity measurements or interior models will ever reach the required precision equivalent of the DSG impact of only a few percent of less.


Introduction
For the first time, the high precision of gravity measurements by the Juno mission at Jupiter and the Cassini Extended Mission at Saturn allow the detection of the tiny perturbations related to the fierce zonal winds in the outer envelopes. However, there is an ongoing dispute about the appropriate equation for linking gravity perturbations and zonal flows (Cao and Stevenson, 2017;Kong et al., 2018;Kaspi et al., 2018). A particular matter of debate is whether the backreaction of the gravity perturbations on the flow dynamics has to be taken into account. This article addresses the question with a new semi-analytical approach.
The impact of gravity on the flow dynamics is generally given by the Navier-Stokes equation. The hydrostatic solution decribes the zero order balance between pressure gradient and effective gravity that defines the fundamental background state. The effective gravity is the sum of gravity and the centrifugal force due to the planetary rotation. Respective equipotential surfaces coincide with surfaces of constant pressure and density and different methods have to devised for finding the respective solution (Zharkov and Trubitsyn, 1978;Wisdom, 1996;Hubbard, 2013;Nettelmann, 2017).
The centrifugal forces lead to a spheroidal deformation of equipotential surfaces and density distribution ρ. The gravity potential thus acquires equatorially symmetric contributions of even degree = 2n with n = 1, 2, 3, .... Here G is the gravity constant, M the planetary mass, R the planetary radius, θ the colatitude, and P a Schmitt-normalized Legendre Polynomial of degree . The gravity harmonics J are given by the volume integral and describe deviations from the spherically symmetric zero order contribution. The degree of rotational deformation depends on the relative importance of centrifugal forces to gravity, which can be quantified by q = Ω 2 /(Gρ), where Ω is the planetary rotation rate. For Jupiter, q remains below 0.1 and deviations from the spherically symmetric gravity thus amount to only about 5%. For Saturn, q is about two times larger than for Jupiter, which is consistent with the stronger deformation of the planet. Since gravity mostly originates from the higher densities in the deep interior, where the deformation is smaller, the deviation of spherical gravity is only slightly larger than for Jupiter. Some of the classical methods for solving the rotationally deformed hydrostatic solution can be extended to include geostrophic zonal flows, which depend only on the distance to the rotation axis (Hubbard, 1982;Kaspi et al., 2016;Wisdom and Hubbard, 2016;Galanti et al., 2017;Cao and Stevenson, 2017). Cao and Stevenson (2017) explore geostrophic zonal flows that are reminiscent of Jupiter's equatorial jet. They report that the zonal wind induced gravity amounts to only three permil of the gravity induced by the planetary rotation for J 2 . For J 8 , both effects have a comparable magnitude, while zonal wind effects dominate for larger degrees. For J 20 , the related contribution is ten orders of magnitude larger than its rotational counterpart. Cao and Stevenson (2017) point out the the small contributions at low degrees can easily be offset by uncertainties in the background model, for example the composition, the equation of state, or the presence of stably stratified layers (Debras and Chabrier, 2019). In practice, the even harmonics up to J 4 , possibly even J 6 , serve to constrain the zero order background state. Only contributions beyond J 6 could thus reliably be exploited to gain information on the equatorially symmetric zonal flows.
The situation changes for the equatorially antisymmetric gravity harmonics, which can be interpreted directly in terms of a first order dynamic perturbation. (The hydrostatic background state being equatorially symmetric and of zero order.) The effect of non-geostrophic flows is estimated based on a simplified dynamic balance. Viscous forces are negligible in the Gas giant atmospheres. Since the zonal winds are rather stable and significantly slower than the planetary rotation, inertial forces are also significantly small than Coriolis forces, buoyancy, or pressure gradients. When taking the curl of the force balance, the pressure gradient also drops out and the first order balance reads where z is the distance to the equatorial plane, Ψ e the effective background potential, ρ the background density, ρ the density perturbation and Ψ the gravity perturbation. Note that we have also neglected the Lorentz-force related term here. While Lorentz forces may play a significant role at depth where electrical conductivities are higher, the are much less important in the outer envelope where zonal flows are fast but electrical conductivities drop to zero. An important point of debate is whether the term involving the gravity perturbation Ψ yields a significant contribution or can be neglected. We refer to this term as the Dynamic Self Gravity (DSG) here. When the DSG can be neglected, the balance (3) reduces to the classical Thermal Wind Equation (TWE). The full balance including DSG has thus been called Thermo-Gravitational Wind Equation (TGWE) by Zhang et al. (2015).
One group of authors insists that the DSG term can be as large as the term involving ρ (Zhang et al., 2015;Kong et al., 2016Kong et al., , 2017Kong et al., , 2018. They also point out that neglecting the DSG would fundamentally change the mathematical nature of the solution. To explore the DSG impact, Kong et al. (2017) assume a zonal wind system that reproduces the observed equatorially antisymmetric winds at Jupiter's cloud level and retains a geostrophic wind morphology at depth, i.e. the morphology is continued downwards along the direction of the rotation axis. Their amplitude, however, is supposed to decay linearly with the distance to the equatorial plane z. They report that neglecting the DSG has a surprisingly large impact on J 1 , and reduces J 3 , J 5 , and J 7 by 25%, 15%, and 7%, respectively.
A second group of authors argues that the DSG can be neglected (Kaspi et al., 2016;Galanti et al., 2017;Kaspi et al., 2018;Iess et al., 2019). Galanti et al. (2017) explore a simplified equatorially symmetric zonal flow system that matches the main features of the respective flows at cloud level. The wind structure is again continued downward along the rotation axis, but assuming an additional exponential decay with depth. They conclude that the DSG has only a minor impact. However, their figure 6 suggests that the zonal-flow-related J 2 decreases by up to 100% when neglecting the DSG. Guillot et al. (2018) use Jupiter's even gravity harmonics up to J 10 measured by the Juno mission to constrain the planets equatorially symmetric zonal winds. Analyzing a suit of possible background models, they report that J 6 , J 8 and J 1 0 can only be explained when the perturbation related to the zonal winds is taken into account. Using the TWE and assuming the exponentially decaying wind structure by Galanti et al. (2017), Guillot et al. (2018) report that the e-folding depth lies somewhere between 2000 and 3500 km.
The odd gravity harmonics J 3 to J 9 based on Juno measurements were also recently used to constrain the depth of the zonal winds. Kong et al. (2018) use the full TGWE equation while Kaspi et al. (2018) neglected the DSG. Both articles where roughly able to explain the gravity harmonics with equatorially antisymmetric zonal winds that reproduce the observed surface winds. Both also conclude that the winds must be significantly slower than observed at the surface below a depth of about 3000 km. However, the suggested radial profiles differ significantly. Since the results rely on different interior models, methods, and assumed zonal flow profiles, it is difficult to judge to which to degree the results are influenced by the DSG. Iess et al. (2019) explore Saturn's even gravity harmonics J 2 to J 10 measured by the Cassini mission. Like for Jupiter, J 6 , J 8 and J 10 can only be explained when considering the zonal wind impact. However, unlike for Jupiter, a slight modification of the surface wind structure is required. Iess et al. (2019) report that these modified winds reach down to a depth of about 9000 km. While generally using they TWE approximation, Galanti et al. (2019) report that J 8 and J 10 increase by about 10% when including DSG in the TGWE approach. Galanti et al. (2019) in addition also analyze the odd harmonics J 3 to J 9 and confirm the inferred depth of Saturn's zonal winds.
Here we explore the relative importance of the DSG with a new (semi) analytical method. Sect. 2 introduces the differential equations that define the gravity potential. Sect. 3 then develops the solution method. Sect. 4 discusses solvability aspects with some illustrative solutions and Sect. 5 quantifies the relative impact of DSG. The paper closes with a discussion in Sect. 6.

From Navier-Stokes Equation to
Inhomogeneous Helmholtz Equation The link between the dynamics and gravity is provided by the Navier-Stokes equation where u is velocity,ẑ the unit vector in the direction of the rotation axis, p the pressure, j the electric current, B the magnetic field, ν the kinematic viscosity, and S the traceless rate-of-strain tensor for constant kinematic viscosity: The effective gravity g e can be expressed by an effective gravity potential, which is the sum of the gravity potential obeying the Poisson equation and the centrifugal potential with s = r sin θ being the distance to the rotation axis. The zero order force balance is given by the hydrostatic equilibrium with vanishing flow and magnetic field: Overbars mark the hydrostatic and non-magnetic background state, while primes denote the perturbation, except for flow and magnetic field. Linearizing with respect to the perturbations yields The linearized buoyancy term has two contributions, one due to the density perturbation and a second one due to the perturbation in gravity. The latter can be separated into a conservative part, written as a gradient, and the remaining contribution: In order to address the zonal-wind related effects, one considers the curl of the Navier-Stokes equation (11) where the pressure gradient and the conservative part of (13) drop out. The approximation motivated in the introduction suggest to neglect inertia, viscous effects, and the Lorentz force contribution: The next step is to assume that ψ Ω can be neglected in comparison to the background gravity contribution Ψ , as discussed in the introduction. The background state then becomes spherically symmetric and equation (14) simplifies to This is the thermo-gravitational wind equation (TGWE) solved for a given U φ for example by Zhang et al. (2015) or Kong et al. (2018). The equation assumes the form of a classical thermal wind equation (TWE) when neglecting the DSG, ρ∇Ψ , or more precisely its non-conservative contribution.
Integrating equation (15) in latitude, dividing by background gravity g = −∂Ψ /∂r, and using equation (12) finally yields an equation that connects the perturbation in the gravity potential to the z-gradient of the zonal winds: with and the dynamic density perturbation as a source term. Note that ρ U is an auxiliary variable different from ρ . We will refer to µ(r) as the DSG coefficient. This second order differential equation must be supplemented by boundary conditions. Solving for solutions in a full sphere, we demand that Ψ vanishes at r = 0. Outside of the source, the solutions must obey A respective matching condition at the outer radius R yields the second boundary condition that we provide further below. Because ρ U is axisymmetric, we will only consider axisymmetric solutions. The integration in latitude means that equation (16) is only determined up to an arbitrary function of radius. This function could only contribute to the spherical symmetric gravity contribution which, outside of the planet, is determined by its total mass and thus carries no information on the dynamics.
The case of the TWE is easy to deal with. Neglecting the DSG implies ρ U = ρ and one simply has to solve the classical Poisson equation (12). The case of the TGWE is more complicated. Using equation (1) and equation (2) transforms the TGWE into the complicated integro-differential equation for ρ derived by Zhang et al. (2015) and the Possion equation for Ψ is then solved in a second step. Their solution is cumbersome and numerically time-consuming.
We avoid this complication by directly solving the inhomogeneous Helmholtztype equation (16) to obtain ψ . The true density perturbation can be recovered by which is obtained from equation (12) and equation (16). We note that ρ U is identical to the 'effective density' that had been introduced by Braginsky and Roberts (1995) in the context of geodynamo equations. They showed that using this variable is an elegant way of dealing with self-gravity, which greatly simplifies that system of equations to be solved. What would be a realistic DSG coefficient µ? Typical textbook density and pressure profiles consider polytropes with index unity. They not only seem to provide reasonable approximations for Jupiter's interior, as is illustrated in Fig. 1, but also yield an analytical expression of the background density and gravity. The former is given by where ρ c is the density at r = 0, and χ a rescaled radius: The gravity profile is then and the DSG coefficient becomes constant: Panel a) of Fig. 1 compares the pressure profile in the Jupiter model by Nettelmann et al. (2012) and French et al. (2012) with a polytrope with index unity, illustrating that this indeed provides a good approximation. More generally, for an adiabatic background state, the density gradient can be written in terms of a pressure gradient: with being the compressibility at constant entropy. Combining equation (25) and equation (9) shows that the gradient in the background density is given by The DSG coefficient is thus given by Panel b) of Fig. 1 compares the constant expression (24) for the index-unity polytrope (dashed line) with the profile (28) based on ab-initio equation-of-state simulations and pre-Juno gravity data (French et al., 2012). Considering the strong variation of other thermodynamic quantities, the µ(r) variations remain remarkable small. In the lower layer r < 0.25R, µ(r) is nearly constant and close to π 2 /R 2 . In the outer envelope r > 0.85 R, µ becomes more variable, reaching amplitudes 40% larger than π 2 /R 2 . A constant µ value thus seem to provide a decent approximation and will considerably ease the task of solving the inhomogeneous Helmholtz equation, as we will discuss in Sect. 3.

Solving Poisson and Inhomogeneous Helmholtz Equations
We start with briefly recapitulating the Green's function method for solving the Poisson equation in Sect. 3.1. Sect. 3.2 then discusses the adapted approach for solving the inhomogeneous Helmholtz equation with constant DSG coefficient µ. The involved methods represent textbook knowledge, but their application to the specific gravity problem is new however, we nevertheless discuss them in some detail.

The Classic Green's-Function Solution
A common way of solving the Poisson equation (7) is the Green's function method. The respective Green's function Γ is defined by where vectors r andr denote the location of potential and density, respectively. The Green's function also has to fulfill the same boundary conditions as the gravity potential. The solution is then given by the integral where denotes the integration over the spherical volume. The classical Green's function for the Poisson problem is given by  2012) is a three layer model with a rocky core that occupies the inner 10% in radius and two gaseous envelopes, above and below 0.625 R, which differ in the metallicity (fraction of elements heavier then helium). Panel b) compares the normalized DSG profile µ(r) R 2 (solid line) suggested by the ab-initio data points by French et al. (2012) (circles) with the constant value π 2 expected for the polytrope (dashed line).
but of more practical use is the representation where Γ is expanded in eigenfunctions of the Laplace operator. Since the Legendre polynomials are eigenfunctions of the horizontal part of the Laplace operator, they are a natural choice to describe the latitudinal dependence: The Schmitt normalization assumed here means that The two possibilities for the radial function are f (r) = r and f (r) = r −( +1) . The expanded Green's function then reads where r > (r < ) denotes that larger (smaller) of the two radii r andr. The matching condition to the field for r > R reduces to the mixed boundary condition which is obviously fulfilled by the radial ansatz functions and thus by the Green's function.
Plugging the Green's function into equation (30) then shows that the potential field for r > R is given by with the expansion coefficients This is equivalent to the differently normalized classical expansion equation (1) and equation (2). The same solution applies to Ψ when replacing ρ by ρ . Should the impact of DSG µ be negligible, we could simply use ρ ≈ ρ U , an approach generally followed by one group of authors mentioned in the introduction (Kaspi et al., 2016;Galanti et al., 2017;Kaspi et al., 2018;Iess et al., 2019;Galanti et al., 2019).

Solving the Inhomogeneous Helmholtz equation
For constant µ(r) = K 2 , the modified potential field equation becomes an inhomogeneous Helmholtz equation The respective Green's function is now defined by and has to fulfill the boundary conditions. Like for the classical Green's function solution discussed in Sect. 3.1, we are looking for a solution in terms of orthonormal functions. While Legendre polynomial can once more be used for the horizontal dependencies, the radial functions have to be different. We will rely on eigenfunctions f (r)P (θ) of the Laplace operator where the f (r) fulfill the boundary conditions.
An orthonormal set of such radial functions can be constructed from spherical Bessel functions (Abramowitz and Stegun, 1984), which solve the differential equation We only use the spherical Bessel functions of the first kind, j , with > 0 that all vanish at r = 0. Spherical Bessel functions of the second kind diverge at the origin, while j 0 (r = 0) = 1. Simple rescaling of the argument yields eigenfunctions of the Laplace operator: with eigenvalues λ = −k 2 n .
The different k n are chosen so that j (k n R) fulfills the boundary condition (36). Because of recurrence relation (88) (see App. C), this condition reduces to which means that the k n are the roots of j −1 (x) divided by the outer boundary radius R. We start the numbering at the smallest root larger than zero so that 0 < k 1 < k 2 < k 3 < .... Panel (a) of Fig. 2 illustrates the spherical Bessel functions j for different degrees . Table 1 list the first five roots for ≤ 5.
Since the Laplace operator is hermitian (adjoint) and our radial ansatz functions fulfill the boundary conditions, the eigenvalues are real and the eigenfunctions for different eigenvalues are orthogonal. For completeness, we include this textbook knowledge is App. A. The orthonormality condition thus reads N n N n ∫ dr r 2 j (k n r) j (k n r) = δ n,n , where the N n are normalization constants derived analytically in Sect. B: Panel (b) of Fig. 2 shows the first five normalized functions, for = 2.
We can now expand the potential field perturbation in Legendre polynomials and the new orthonormal radial functions: Using this expansion in equation (39), multiplying with the ansatz functions j n (r)P (θ) and integrating over the volume yields a spectral equation for the expansion coefficients: The coefficients are thus simply given by A comparison with equation (35) shows that the Green's function for the inhomogeneous Helmholz equation is then The potential field for r > R has to decay like (R/r) +1 . The respective solution is thus given by with As expected, this solution is identical to the classical result (37) for K 2 = 0. We show this analytically in Sect. D.

Illustrative Examples
We can easily convince ourselves that equation (48) with coefficients (50) provides a correct solution when assuming that the source is given by only one ansatz function: Only the respective potential field coefficient thus has to be considered and the solution for r < R is Solving for a more general source thus boils down to the question: How well can ρ U be expanded in the ansatz functions? A special situation arises when K 2 = k 2 n . For the polytropic density distribution with polytropic index unity, this happens for = 1 and n = 1 where K = k 1,1 = π/R. The two non-conservative buoyancy terms then cancel exactly, because of matching radial functions in the background profiles and the primed perturbations. Nothing is left to balance the respective left hand side of the simplified dynamic equation (15) or the related ρ U contributions in (16). The respective potential field perturbation thus decouples from the simplified dynamical equation. Even when K is not identical but close to k 1,1 , the dynamic equation requires an unrealistically large potential field perturbation and the precise value of K would have an enormous effect. It thus seems a good idea to generally avoid these resonance conditions and we will simply not interpret respective Ψ 1,1 contributions. Since the = 1 gravity contribution generally vanishes due to the choice of origin r = 0, these considerations are of little practical use.
Partial integration of the dynamic density perturbation yields   Zhang et al. (2015). Column 2 and 3 list TWE and TGWE results for the slower decaying flow with h = 0.143. Column 4 and 5 list respective values for the faster decaying case h = 1.143 also illustrated in Fig. 3. The last line lists the values published by Zhang et al. (2015).
While latitude-dependence is this purely determined by the zonal flow, ρ, U φ and their radial derivatives influence the radial profile of ρ U . Since the expansion of the latitude-dependence in Legendre polynomials is not specific to solutions with or without DSG, we concentrate on discussing the expansion in radius. The steep radial gradients in density and zonal flows characteristic for gas planets may prove challenging here.
Choosing a truncation N for the radial expansion defines the numerical representation of ρ U : with and The quality of the representation is quantified by the misfit We start with exploring a test case suggested by Zhang et al. (2015). They assume the polytrope index unity density profile (21) and a zonal flow defined by with amplitude U 0 = RΩ/100 and radial dependence Jupiter values used to define flow and gravity are R = 6.9894×10 7 m, Ω = 1.759×10 −4 s −1 , and M = 1.898×10 27 kg. Two relative decay scale heights h = 0.143 and h = 1.143 are explored. The flow yields = 0 and = 2 gravity perturbations, but since the former would be nonphysical in a real gravity problem we only consider the latter. Table 2 compares the respective J 2 coefficients published by Zhang et al. (2015) with values for different truncations N. While the results for h = 0.143 exactly match those of Zhang et al. (2015), those for h = 1.143 already differ in the second figure. We attribute this to convergence problems reported by Zhang et al. (2015).
The well behaved convergence for the expansion of f 1 (r) is documented in Table 2 and illustrated in Fig. 3. Panel a) and b) demonstrate that the function is already almost perfectly represented with a truncation of N = 40. Small differences tend to remain close to the outer boundary and at small radii due to the specific properties of the j n . Spectrum and misfit M, depicted in panels c) and d) respectively, decay continuously with truncation but with a slower rate at higher degrees because of the difficulties in exactly capturing the vanishing values for r → 0.
As a second example we explore the function used in the classical potential field solution for K 2 = 0. This is an ideal test case, since the expansion coefficients are known analytically (see App. D). Fig. 4 illustrates the quality of the expansion for = 3. Panels a) and b) once more illustrate the difficulties of representing the function at the boundaries. The last example is the radial function that determines the radial dependence of one term in ρ U according to equation (57). Following the example of Kong et al. (2018), we assume a polytrope of index one and the Gaussian-like flow profile: where d = R − r is the depth, D = 0.15 R is the maximum depth of ρ U , and h = 0.22 determines the decay rate. Fig. 5 demonstrates that the resulting highly localized function is also already well represented for a truncation of N = 40. Overall, spectrum and misfit once more decay with growing N, which confirms that there are no principal numerical problems with expand this demanding function into the j n . The pronounced length scale defined by the width of the function peak leads to the local minima in the spectrum where they match the distance between the zero intercepts in the j n .

Relative Importance of Dynamic Self Gravity
The analytical solution shows that the impact of the DSG simply depends on the ratio k 2 n /K 2 . The relative importance of K 2 in the inhomogenous Helmholtz equation for a given spherical harmonic degree and radial index n can be quantified by Table 3 lists S n for spherical harmonic degrees up to = 30 and n up to 5, assuming K = π. The values indicate that the DSG should be considered a first order effect for ≤ 4, reaches the 10% level at = 5 or = 6 and amounts to only about 1% for ≥ 20.

Discussion and Conclusion
The dominant balance between the Coriolis force and buoyancy terms in the azimuthal component of the vorticity equation establishes a connection between zonal flows and gravity. Simple manipulations lead to what has been called the thermo-gravitational wind equation (TGWE) by Zhang et al. (2015). This contains two buoyancy contributions: one related to the density perturbation and a second that we named dynamics self gravity (DSG) since it directly links the disturbed gravity potential and zonal flows.
The dynamic perturbation of the gravity potential Ψ is defined by the inhomogeneous differential equation where µ is the DSG factor and ρ U is the source term describing the impact of the zonal flows. The only difference to the classical Poisson equation for a gravity potential is the DSG term. The dynamic density perturbation ρ U , which is identical to the effective density introduced by Braginsky and Roberts (1995), is obtained from zonal flow and background density by a simple integral. A polytrope of index unity offers a reasonable approximation for the interior of Jupiter and other gas planets. This implies that µ = π 2 /R 2 is constant, which considerably eases the task of solving equation (69). The problem then assumes the form of an inhomogeneous Helmholtz equation and the solution becomes particularly simple when expanding the radial dependence in modified spherical Bessel functions that fulfill the boundary conditions. Like in the classical gravity problem, Legendre polynomials remain the representation of choice for the latitudinal dependence. These basis functions allow a very efficient (semi) analytical solution to the problem. Each of the calculations presented here required only a few seconds of run time on a standard 4-core notebook.
There has been a discussion whether the DSG term could be neglected when inverting high precision gravity observations at Jupiter and Saturn for zonal flow properties. Our new formulation allows us to quantify the relative impact of the DSG for each gravity harmonic, practically independent of the considered zonal flow or background state.
A special case arises for degree = 1. For the background density with polytropic index unity, the = 1 solution comprises the case where the two buoyancy contributions in the TGWE cancel. This corresponds to the homogeneous solution of the Helmholtz equation. Zonal flow and gravity perturbation then decouple, and it becomes impossible to draw on the zonal flows from the respective gravity contribution. Kong et al. (2017) seem to have noticed the related problems without realizing their origin. However, this is of little practical interest since the origin is generally chosen to coincide with the center of gravity so that = 1 contributions vanish.  . The fourth column shows S , the relative impact of the DSG for radial profile f 3 also listed in Table 4.
Jupiter's J 3 , J 5 and J 7 coefficients, the relative impact of DSG is comparable to the error and should thus be taken into account when inverting gravity harmonics for zonal flow properties. This agrees with the results and conclusion by Kong et al. (2017). The error of the higher order harmonics may decrease as the Juno mission progresses. For Saturn, J 3 , J 5 and J 10 seem precise enough to warrant including DSG effects. The estimates of Kong et al. (2017) and Galanti et al. (2019) about the relative impact of the DSG is compatible with our results. Including the DSG term generally increases the amplitude of the gravity coefficients. As pointed out by Galanti et al. (2017) and Cao and Stevenson (2017), including the rotational deformation of the background density in the TWE or TWGE approaches may have a similar relative impact on the odd gravity harmonics as the DSG. Both effects may thus have to be taken into account when trying to explain these harmonics by the zonal wind dynamics.

A Orthogonality
In this section we show that the spherical Bessel functions for different k n are orthogonal and that k 2 n is real. We start by recalling the properties of a selfadjoint or Hemitian linear operator L. Let f and g be eigenvectors (functions) of L with eigenvalues λ and µ: For a self-adjoint operator we have It follows that λ g, f = µ g, f and thus λ = µ . The eigenvalue is thus real and for λ µ we must have Here the angular brackets denote the integration over the interval of interest, in our case To show under which conditions an operator is Hermitian, we chose a somewhat more general textbook example: L = a(r) ∂ 2 ∂r 2 + b(r) ∂ ∂r + c(r) .
Partial integration yields f , L g = r 2 a f ∂g ∂r + r 2 b f rg − g ∂(r 2 a f ) ∂r R 0 + ∫ R 0 dr g ∂ 2 (r 2 a f ) ∂r 2 − ∂(r 2 b f ) ∂r + r 2 f c (76) Rewriting part of the last integral in terms of the operator L leads to f , L g = L f , g + r 2 a f ∂g ∂r The remaining integral vanishes when ∂(r 2 a) ∂r = r 2 b , Using recurrence relation (88) and L'Hopital's rule yields Finally, using recurrence relations (89) leads to ∫ R 0 dr r 2 j 2 (kr) = R 3 2 j 2 (k R) and thus the normalization constant N n = 2 1/2 r 3/2 o j (k n r o ) . (87)

D Equivalence of new and classical solution
For K 2 = 0, both the classical solution equation (37) where j n = N n j (k n r).
In order to show that this is indeed true, we expand the radial dependence under the integral in the classical solution into our set of orthonormal spherical Bessel functions:r = ∞ n=1 j n (r) ∫ R 0 d r r +2 j n (r) .
Partial integration and using the boundary conditions (36) yields Using recurrence relation (88) and performing another partial integration finally gives ∫ R 0 d r r +2 j n (r) = (2 + 1) Plugging this into equation (92) and then the result into the left hand side of equation (91) finally proves equation (91).