Linking Zonal Winds and Gravity II: explaining the equatorially antisymmetric gravity moments of Jupiter

The recent gravity field measurements of Jupiter (Juno) and Saturn (Cassini) confirm the existence of deep zonal flows reaching to a depth of 5\% and 15\% of the respective radius. Relating the zonal wind induced density perturbations to the gravity moments has become a major tool to characterise the interior dynamics of gas giants. Previous studies differ with respect to the assumptions made on how the wind velocity relates to density anomalies, on the functional form of its decay with depth, and on the continuity of antisymmetric winds across the equatorial plane. Most of the suggested vertical structures exhibit a rather smooth radial decay of the zonal wind, which seems at odds with the observed secular variation of the magnetic field and the prevailing geostrophy of the zonal winds. Moreover, the results relied on an artificial equatorial regularisation or ignored the equatorial discontinuity altogether. We favour an alternative structure, where the equatorially antisymmetric zonal wind in an equatorial latitude belt between $\pm 21^\circ$ remains so shallow that it does not contribute to the gravity signal. The winds at higher latitudes suffice to convincingly explain the measured gravity moments. Our results indicate that the winds are geostrophic, i.e. constant along cylinders, in the outer $3000\,$ km and decay rapidly below. The preferred wind structure is 50\% deeper than previously thought, agrees with the measured gravity moment, is compliant with the magnetic constraints and the requirement of an adiabatic atmosphere and unbiased by the treatment of the equatorial discontinuity.


INTRODUCTION
The spacecrafts Juno and Cassini delivered high-accuracy measurements of the gravity potentials for Jupiter and Saturn (Kaspi et al. 2018;Iess et al. 2018;Guillot et al. 2018;Galanti et al. 2019;Iess et al. 2019), which provide valuable constraints on the interior structure and dynamics of their atmospheres. For the first time, it has been possible to resolve the tiny undulations in the gravity potential induced by zonal winds. This has ended the long-standing debate on whether the zonal winds on Jupiter and Saturn are shallow weather phenomena or reach deeper into the planets' convective envelopes. The gravity data suggest that the winds extend down to 5% and 15% of Jupiter's and Saturn's radii, respectively (Kaspi et al. 2018;Galanti et al. 2019;Kaspi et al. 2020). However, modelling the deep-reaching zonal mass fluxes on a gas planet with all their complexity and relating them to anomalies in the gravity field is a difficult, non-unique problem.
The recent attempts of constraining the deep-reaching winds with gravity data are based on simple parametrizations of the vertical wind structure. The observed surface winds are thus con-★ Contact e-mail: dietrichw@mps.mpg.de tinued downward along cylinders that are aligned with the rotation axis. The resulting geostrophic flow is then multiplied with a parametrized radial profile to model the decay with depth. The parameters of this profile are constrained by the gravity measurements (Kaspi et al. 2010(Kaspi et al. , 2013(Kaspi et al. , 2018Kong et al. 2018;Galanti et al. 2019). More recently, constraints from Jupiter's magnetic field (Duer et al. 2019) or it's temporal evolution  have also been used in addition to this. However, some studies found it necessary to partially modify the zonal flow compared to what is observed at cloud level (Kong et al. 2018; in order to match the modelled gravity anomalies to the observed ones. Gravity moments of odd degree exclusively contain the impact of equatorially antisymmetric zonal flows, since they are not obscured by the rotational deformation of the planet as opposed to their even counterparts. The geostrophic, downward continuation of the antisymmetric zonal flows from each hemisphere yields opposite signs when they each reach the equatorial plane and hence introduces a problematic discontinuity or equatorial 'step'. This has been reported to potentially bias the results (Kong et al. 2016). To circumvent this issue, several published studies, such as Kaspi et al. (2018);  calculate the antisymmetric gravity signal in each hemisphere, but ignore the resulting equatorial shear region. Alternatively Kong et al. (2017); Kong et al. (2018) suggested to smooth the equatorial region artificially.
Zonal winds induce perturbations in the gravity potential via pressure perturbations that change the density distribution. The governing leading order Navier-Stokes (or Euler) equation, a balance between pressure gradient; Coriolis force; and gravity forces; establishes the link between zonal flows and the induced density perturbation. The gravity term has two contributions, one due to the zonal-flow induced density perturbation and a second due to the dynamic gravity perturbation. If the latter, termed dynamic self gravity (DSG) by , is ignored, the balance reduces to the classic thermal wind equation (TWE) which can directly be solved for the density perturbation and hence the gravity perturbation. This approach has been commonly used for modelling the gravity anomalies of Jupiter and Saturn and characterising their deep interior structure (Guillot et al. 2018;Kaspi et al. 2018;Iess et al. 2019;Galanti et al. 2019). Several authors claim that the DSGterm is indeed negligible (Galanti et al. 2017;Kaspi et al. 2018), but this has been disputed by Zhang et al. (2015), Kong et al. (2017), and .
Retaining the DSG-related term yields the so called thermogravitational wind equation (TGWE), which is harder to handle mathematically and numerically. Zhang et al. (2015) reformulate the TGWE into an integro-differential equation for the density perturbation. Their approach is restricted to polytropic interiors with index unity. More recently,  show that the TGWE can be formulated as an inhomogeneous Helmholtz equation for the gravity potential. The results show that the relative impact of the DSG decreases with the spherical harmonic degree ℓ of the gravity perturbation. Being of order one for ℓ = 2, it decreases to 10% at about ℓ = 5 and becomes smaller than 1% for ℓ > 25. However, the polytrope provides only a crude approximation of Jupiter's interior.
Here we introduce a solution method for the TGWE that is based on formalism of  but can also handle more realistic interior models that are reflected by a radius-dependent DSG coefficient.
The discrepancy between existing studies of Jupiter's odd gravity moments is related to differences in the modelling assumptions. Kaspi et al. (2018) explained the first four odd gravity moments of Jupiter by extending the observed zonal wind profile obtained during the perĳove (Tollefson et al. 2017) along cylinders, inverting the TWE equation, assuming a realistic Jupiter interior model (Guillot et al. 2003) and a relatively smooth radial decay while ignoring the discontinuity at the equatorial plane (see fig. 1, green). This has been challenged by Kong et al. (2017), arguing that the thermo-gravitational wind equation (TGWE) should be used. The subsequent application to Jupiter by Kong et al. (2018) is based on a polytropic interior and the surface flow measurements from the Cassini mission (Porco et al. 2003). In order to carry this out the surface zonal flow was also modified by putting a cap on the amplitude of the spike in the antisymmetric flow at = 21 • latitude, and applied a second decay function to smooth out the discontinuity at the equatorial plane (orange profile).
Both studies favour a rather smooth decay with depth in order to match the observed gravity signal, which implies quite significant flow amplitudes below 0.95 . The radial decay profile of Kaspi et al. (2018) implies a remaining 10% flow amplitude at 0.95 , whereas Kong et al. (2018) found more than half of the surface amplitude remaining at this depth. The weak secular variation of Jupiter's observed magnetic field over the past decades, however, can only be explained if the amplitude of the zonal wind at the depth of the horizontal advection of the radial field, i.e. 4% to 6% of Jupiter's radius, is as weak as cm/s compared to tens or hundreds of meters per second at the surface (Moore et al. 2019) (black dots in the figure). More recently,  showed that a combined gravito-magnetic analysis favoured a sharper radial profile with a 2000 km deep barotropic part and 600 km decay region (red profile in fig. 1). However, also this analysis neglected the flow discontinuity at the equatorial plane.
In this study we revisit the key model assumptions and properties on which the analyses of odd-degree gravity moments are based and quantify their impact on the solutions. This includes justifying the fundamental equation, the problematic treatment of the equator, various surface zonal flows profiles and different models of Jupiter's interior.

Governing equation
The gravity potential is directly related to the density distribution via the Poisson equation where is the gravitational constant. The general solution is sin˜˜˜ (2) and the associated gravity force is = −∇ . It is useful to separate the external gravity potential into a spherically symmetric part of a non-rotating, solid body and a series of higher order terms originating from density perturbations and the rotational deformation of the planet: where is the equatorial planetary radius, the total mass and ℓ are the Legendre polynomials of degree ℓ. The ℓ = 1 contribution vanishes because the origin of the coordinate system has been chosen to coincide with the center of gravity. The gravity moments ℓ are given by: While the odd moments contain the signal of the equatorially antisymmetric component of zonal flows, the even moments are dominated by the effects of the rotational deformation of the planet. We therefore focus on the equatorially antisymmetric winds and the respective observed odd gravity moments 3 , 5 , 7 and 9 . The general relation between the zonal mass flux and a density anomaly is expressed by the reduced Navier-Stokes equation, that describes the conservation of momentum for a steady, inviscid, non-magnetic, inertia-less flow rotating around theˆ-axis with a rotation rate Ω.
In the co-rotating frame of reference this reads: where the terms (from left to right) are the Coriolis force, the pressure gradient, the gravity and the centrifugal force. This leading order force balance applies to the quasi stationary zonal flows in the outer envelope where the electrical conductivity is so low that the Lorentz force can be neglected. The centrifugal force in the leading order force balance represents the rotational deformation and renders the steady background state two-dimensional and non-spherical (e.g. Cao & Stevenson 2017). However, the rotational deformation itself is rather insignificant for the antisymmetric problem (Kong et al. 2016) and thus we ignore the centrifugal forces for now. Pressure, density and gravity are separated into a hydrostatic background that depends only on radius and a small perturbation, e.g. = ( ) + ( , ). The first order perturbation equation is then given by: where the last two terms on the right hand side are gravity force contributions due to a dynamic (i.e. flow-induced) density anomaly ( ∇ ) and the dynamic self-gravity (DSG, ∇ ) term. In the classic thermal wind approach (e.g. Kaspi et al. 2018), the DSG is neglected. The density anomaly can then simply be found from the thermal wind equation (TWE), which is the azimuthal component of the curl of eq. 6: An integration along latitude and division by the background gravity yields the anomalous density field, which is subsequently used to calculate ℓ by eq. 4. Zhang et al. (2015) and  show that the DSG term represents a first order effect and, for example, changes the Replacing by eq. 2 and integrating over latitude leads to the integro-differential equation solved by Zhang et al. (2015) 2Ω While the integration function ( ) renders mathematically nonunique, the gravity moments nevertheless remain unique (Kaspi et al. 2010;Zhang et al. 2015). The integration function ( ) would only contribute to the spherical symmetric gravity harmonic which is determined by the total mass and therefore we can set ( ) to zero. Treating eq. 7 is mathematically demanding and has so far only been solved for a simple interior model, i.e. a polytrope of index unity (Zhang et al. 2015;Kong et al. 2018). A potential work-around was devised by Braginsky & Roberts (1995) in the framework of classic geodynamo theory. It was shown that the DSG term can be absorbed into the so-called effective variables (density and pressure): The radial function is thereby characterised by the compressibility of the considered medium: where is the sound speed. For a polytropic perfect gas this can be further simplified, where is the polytropic index, and are density and pressure at the centre of the planet. They depend on the specific solution of the Lane-Emden equation. Moreover, for a polytropic index unity, is a constant and amounts to 2 / 2 . The effective variables, and , are equal to and reduced by the contribution of the local elevation or depression of the associated equipotential surface thus capturing the effect of the DSG. Using the effective variables and the definition of (eq. 13), the Navier-Stokes equation (eq. 6) simplifies to: Taking the azimuthal component of the curl, integrating along colatitude and replacing ∇ = , then yields the effective density perturbation: Since this effective density disturbance formulates the zonal flow impact, it has been called by  and we adopt this name here. Note that this is not the true density disturbance that could serve to calculate ℓ via eq. 4. In particular, and are related as defined by eq. 12. Now we can find an equation for the gravity potential by replacing with using eq. 1 in eq. 9, integrating along latitude and making use of eq. 16: This is a two-dimensional, inhomogeneous PDE of second order, which describes the wind-induced anomalies in the gravity potential of a gas planet . The effective density perturbation derived from eq. 16 acts as the source term for this Helmholtz-like equation. Only when = does this equation become an inhomogeneous Helmholtz equation and can be solved in a semi-analytical way ). In the more general case where depends on the radius we refer to the numerical methods discussed in sec. 4.
It is important to note that the 2nd order differential equation (eq. 17) and the integro-differential form of the TGWE (eq. 10) describe the same physical problem. The main difference is that eq. 17 solves for , while eq. 10 solves for the density anomaly . Eq. 17 not only directly provides the gravity potential we are interested in, but is also much easier to solve.

Geostrophy of the zonal winds
The effective variables can be further exploited to show that the flow ( ) and not the mass flux ( ) is geostrophic and should be initially extended along cylinders. To emphasise under which conditions geostrophic winds can be modulated along the rotation axis, we express the density as a function of pressure and entropy: Dividing the Navier-Stokes equation (eq. 15) by the background density and defining a reduced pressure ★ = / , eq. 15 yields: Note, the term in the brackets scales with the non-adiabaticity of the background state and hence can safely be neglected for a vigorously convecting atmosphere like Jupiter's. This equation highlights that the winds and not the mass fluxes are cylindrically invariant. Only if there are sizeable latitudinal gradients in the zonally averaged entropy, are deviations from geostrophy possible, if we restrict the consideration to regions of negligible electrical conductivity. The associated temperature fluctuations required to drive the zonal wind out of its cylindrical invariance can be estimated by azimuthal component of the curl of eq. 19: where / = has been used. Assuming that the vertical and latitudinal derivative can be approximated with the same length scale (e.g. close to the equator), thus ≈ Δ / and 1/ ≈ / . Furthermore, the entropy fluctuations can be (to first order) approximated by temperature fluctuations, thus ≈ / . This then yields For Jupiter Ω ≈ 1.76 · 10 −4 s −1 and ≈ 25 m/s 2 . Then, in order to induce a vertical variation of the wind Δ = 10 m/s at a temperature of ( = 0.95 ) = 4 · 10 3 K (Nettelmann et al. 2012), we find ≈ 0.5 K -an unrealistically high value considering that the temperature fluctuations associated with convective motions are on the order of 10 −4 K (e.g. Jones 2007). These estimates are applicable to the convective part of the atmosphere, i.e. below = 1 bar. In conclusion, whenever convection restricts the degree of non-adiabaticity, the zonal flows (and not the mass flux) are cylindrically invariant.

PARAMETRIZING THE DYNAMIC DENSITY PERTURBATION
The dynamic density perturbation (eq. 16) is governed by the assumed interior structure of the planet via the background density and gravity, and , respectively, and the -gradient of the zonal mass flux. Theoretical considerations and numerical simulations suggest a geostrophic zonal flow structure for the fast rotation and low viscosity gas planets (e.g. Taylor 1917;Dietrich & Jones 2018;Gastine & Wicht 2021). This means that the flow depends only on the distance = sin to the rotation axis, where is the colatitude. We could then simply downward continue the surface zonal flow 0 along cylinders. However, the gravity measurements and the secular variation of the Jupiter's magnetic field show that the wind speed must significantly decrease with depth and should be confined to outer 5% of the planetary radius (Kaspi et al. 2018;Moore et al. 2019;). This decrease is commonly parametrized with an additional radial decay function , with ( ) = 1. The cause for the deviation from geostrophy remains unspecified. Electromagnetic effects have been alluded to. Buoyancy forces arising in a stably stratified region with the assistance of Lorentz forces are a promising mechanism Gastine & Wicht 2021). We start with discussing the interior state, then the different models of the surface zonal profile and finally the calculation of the dynamic density source .

Interior state
The source term (eq. 26) and the DSG coefficient, , depend on the background density and gravity profile, and , respectively. We explore their influence via the chosen interior state model by comparing three commonly used models. A simple model for Jupiter's interior is a polytropic ideal gas of index unity that provides a decent description of the pressure dependence on the density (Hubbard 1999), but not of the ( − )curve. The background density is then: where is the central density. The associated background gravity is given by and is hence directly proportional to the radial density gradient as discussed in Zhang et al. (2015). This proportionality ( ∝ / ) is an exclusive property of a polytrope with index unity and implies a constant = 2 / 2 . This has been exploited by Zhang et al. (2015) and  to solve the TGWE. The grey profiles in fig. 2 illustrate the density, its radial gradient and for the polytrope of index = 1. The first analysis of the Juno gravity data by Kaspi et al. (2018) was based on a Jupiter interior model by Guillot et al. (1994); Guillot et al. (2003). This model is henceforth named 'G03' and assumes a three layer structure with a helium-depleted, molecular outer hydrogen layer, a helium-enriched metallic inner hydrogen layer and a central dense core. Fig. 2 shows the density, density gradient and for a setup with interpolated hydrogen EOS, a 1-bar temperature of 165 K, a core of 4.2 earth masses and heavy element abundance of 33 earth masses (Guillot et al. 2003). In comparison to the rather simple polytropic model (grey), the densities are substantially higher for < 0.975 and slightly lower above this radius. Thus the density gradient shows a pronounced maximum around = 0.97 ( fig. 2, middle panel). The DSG coefficient ( ) increases with radius and reaches a 2.5 times larger value than for the polytrope.
Alternatively, we use the more recent calculations from Nettelmann et al. (2012) based on the updated H-REOS2 model ('N12'). The DSG coefficient is related to the sound speed (eq. 13), which was calculated for the same Jupiter model by French et al. (2012). For this model (termed J11-8a) the depth of the molecular-metallic phase transition is at = 8 Mbar, whereas in G03 model this happens at shallower 2 Mbar. As shown in the left panel of fig. 2 the N12 model yields the highest density for < 0.92 and falls between the G03 model and the polytrope at larger radii. The density gradient and hence show two distinct maxima with smaller amplitude than the G03 model.

Jupiter's surface zonal flow
The surface flow of Jupiter is deduced from tracking cloud features, either with the Hubble Space Telescope (HST) or from space crafts. Fig. 3 displays several zonal flow measurements that have been obtained over the last decades. The oldest illustrated flows (blue and yellow) are based on observations by Voyager I and II (Limaye et al. 1982;Limaye 1986) and represent the flow in 1979. HST monitored the wind structure several times between 1995 and 1998 (green profile, ?). Cassini measurements of the wind profile stem from late 2000 (red, Porco et al. 2003). Later, various different HST-based profiles were obtained between 2009 and 206, during the perĳove nine of the Juno space craft (Tollefson et al. 2017). Kaspi et al. (2018) and  use this most recent flow profile for their gravity data analysis, while Kong et al. (2018) prefer the model based on Cassini images from late 2000 to early 2001 presented in Porco et al. (2003).
All flow models show the same principle structure but also differ in some details. The equatorially antisymmetric flow contribution shown in fig. 3 is clearly dominated by the prograde jet which starts at a latitude around = 17 • and extends to about = 23 • . The amplitude of this jet varies between 65 m/s and 75 m/s for the different models. Tollefson et al. (2017) conclude that the variations exceed the model errors, at least for some epochs and at some latitudes. Their Lomb-Scargle periodogram analysis reveals dominant variation periods of about 7 yr and 14 yr. Since the latter is close enough to Jupiter's orbital period of 11.9 yr and the former to the respective first overtone, the variations may represent seasonal cycles. These variation are unlikely to penetrate deeper into Jupiter's atmosphere and therefore play no role for the gravity signal. The arithmetic mean profile, shown as a black line in Fig. 3, should represent the deeper flows more faithfully than a single model. However, the average only covers a period of about 36 years, or about three times the seasonal cycle, with a small number of 'snapshots'. Smaller deviations from the time average are thus certainly conceivable. We mostly rely on the mean wind profile but also explore the impact of using specific snapshots on the gravity signal in sec. 5.5.

Treatment of the equatorial discontintuity
An obvious problem arises with the equatorially antisymmetric contributions to the dynamic density source. Geostrophically downward continuing the surface flow in each hemisphere separately yields a discontinuity, or step, at the equatorial plane. It can be large in case the geostrophic continuation of the surface flow hits the equatorial plane at a radius where is significantly larger than zero. This equatorial step seems unphysical but can easily be dealt with mathematically.
Using the product ansatz (eq. 22) in eq. 16, we can separate the dynamic density perturbation into two contributions: The second integral contributes only at the equator where thederivative yields a delta-peak. This can be integrated analytically: where is the Heaviside step function and Δ = 2 ( → /2). As we will show below, the equatorial step can yield an important contribution to the gravity signal which is problematic. Kaspi et al. (2018) and  therefore ignore the respective term in an approach they call 'hemispheric'. They calculate in each hemisphere but dismiss the contribution at (or very close to) the equator, which is equivalent to evaluating  This approach seems to offer a simple fix but is also inconsistent since the step is an integral part of the chosen flow model. Alternatively, Kong et al. (2018) use an additional lineardependence of , such that the vertical gradient at the equator is regular and smooth. Their zonal flow model is thus given by where | | = ( 2 − 2 ) 1/2 is the local distance from the equatorial plane to the surface. The equatorial step and the linear -dependence are the steepest and smoothest end member, respectively, of the possible functions that reconcile the northern and southern zonal flows.
Here we introduce a novel approach guided by the physical principles of atmospheric dynamics. The generation of deepreaching zonal flows powered by the statistical correlations of the convective flows (Reynolds stresses) are best captured in numerical simulations (Heimpel et al. 2005;Christensen 2002;Dietrich & Jones 2018;Gastine & Wicht 2021). Fig. 4 illustrates the zonal flow in a typical simulation of convection in a fast rotating spherical shell with an aspect ratio of / = 0.8 using the MagIC code (Wicht 2002;Gastine & Wicht 2012;Schaeffer 2013). For the chosen parameters ( = 3 · 10 −5 , = 0.25, = 7 · 10 7 , = 2.3, definition as in Gastine & Wicht (2012)), a broad equatorial prograde jet develops, flanked by a deep retrograde jet adjacent to the tangent cylinder (TC). The TC is a virtual cylinder touching the inner boundary at the equator and forms an important dynamical boundary in rotating convection. At higher latitudes, i.e. within the TC, alternating jets develop in both hemispheres. These winds show minimal variations in the direction of the rotation axis. Simulations at even more realistic parameters yield solutions with higher degrees of geostrophy but are also more numerically demanding (Heimpel et al. 2016;Gastine & Wicht 2021).
Panels b) and c) of Fig. 4 show that the equatorially symmetric flow is dominant outside the TC and that the equatorially antisymmetric zonal flow is much weaker outside than inside the TC. Outside of the TC, equatorially antisymmetric flows are at odds with geostrophy. Inside the TC, however, antisymmetric contributions can be, and indeed are, highly geostrophic. The same dynamical reasons that enforce geostrophy (dominant Coriolis force and small viscosity) also allow for only weak equatorially antisymmetric flows outside the TC.
To account for this physical property of zonal flow in rotating spherical shells, we multiply the antisymmetric surface flow with an additional attenuation function that reduces the amplitude outside of the latitude TC , i.e. where the TC touches the outer boundary: A second parameter which is introduced here is the width of the latitudinal cut-off, . Both functions, ( TC , ) and (ℎ, ℎ ), define an individual tangent cylinder and hence should be consistent. More details on selecting the best parameter combination for TC , , ℎ and ℎ are discussed in sec. 5. When TC is sufficiently large and the radial decay function drops rapidly below a certain depth, the resulting flow splits into an independent northern and southern part. Then the equatorial step contribution is identical to zero. The underlying assumption is, that part of the observed zonal flow at cloud level, in particular its antisymmetric parts at low latitude, are shallow. We term cases that apply such an attenuation of antisymmetric surface flows at low latitudes 'TC-models'.

Radial decay function
The different forms of the radial decay function ( ) that have been suggested to explain the gravity observations are illustrated in Fig. 1. Kaspi et al. (2018) use a combination of an exponential decay and a hyperbolic tangent: Here is the weight of the hyperbolic tangent contribution, ℎ the decay depth and ℎ the decay width. For explaining the gravity observations, Kaspi et al. (2018) propose a large = 0.92, a relatively large ℎ = 1570 km and a depth of ℎ = 1803 km. This results in a smoothly decaying radial profile (green) illustrated in Fig. 1. Kong et al. (2018), on the other hand, used an inverse Gauss profile for > , and = 0.0 for . Kong et al. (2018) report that the observations are best matched with a combination of = 10 484 km and ℎ = 15 377 km. Fig. 1 shows that the respective profile (orange) decays also rather smoothly.
Both the solution suggested by Kaspi et al. (2018) and by Kong et al. (2018) are not compatible with the magnetic observations (Moore et al. 2019). Furthermore, vigorous convection provides an almost adiabatic environment and hence the decay functions should be rather flat with a sharp drop at greater depth (see also sec. 2.2). Realizing this,  proposed a steeper decaying alternative with a hyperbolic tangent outer, and an exponential inner branch: They suggest = 0.45, ℎ = 2002 km, / = 0.972 and ℎ = ℎ = 204 km in conjunction with an optimised surface flow profile to explain the gravity observations. The radial profile drops almost faster than the magnetic constraints require (see fig.1).
We adopt the pure hyperbolic tangent profile: For ℎ = 200 km and ℎ = 2000 km the profile would be virtually identical to the one suggested by . However, our analysis favours a substantially deeper flow with ℎ = 3000 km and ℎ = 500 km. Fig. 1 shows the respective decay function in blue.

NUMERICAL METHOD
The semi-analytical method for solving the TGWE (eq. 17) developed by  is restricted to the case of a constant . However, as shown in fig. 2, ( ) varies strongly and reaches much higher values than 2 / 2 , particularly in the outer 10% of Jupiter's radius where the flow-induced gravity moments originate from. We have therefore developed a numerical method that can also handle radial variations in . We first formulate the dynamic density source ( , ) via eq. 26 in spatial space by choosing a radial decay function , an interior state model setting , and , and a surface zonal flow profile , that is cylindrically downward continued. We then expand the latitudinal dependence of the gravity potential perturbation and the dynamic density in Legendre polynomials, e.g. for the former this is given by Since the background state, and hence , are only a function of radius, the TGWE (eq. 17) decouples for each degree ℓ, and we are then left with a set of radial ODEs, one for each spherical harmonic degree ℓ: At the outer boundary the potential must match to the solution of a source free region, i.e. ∇ 2 = 0, with the characteristic radial dependence −(ℓ+1) . This leads to the mixed outer boundary condition: and ℓ (0) = 0 at the inner boundary. Our numerical method relies on an expansion in the modified spherical Bessel functions introduced by . These are radial eigenfunctions of the homogeneous Helmholtz equation. For each spherical harmonic degree ℓ we construct a set of N normalised orthogonal functions ★ ℓ ( ℓ ). The radial scales, ℓ , are determined by finding the first N radii (starting at the origin), at which ★ ℓ ( ℓ ) fulfils the boundary condition (eq. 38). The radial functions ℓ ( ) are expanded in ★ ℓ : where the expansion coefficients ℓ are given by Since the modified spherical Bessel functions are eigenfunctions of the Laplace operator, eq. 37 transforms into a set of linear algebraic equations for the expansion coefficients ℓ : Eq. 41 defines a coupled set of algebraic equations and solved by using a matrix formalism. The source term, ℓ , on the RHS is discretized along radial grid points, . The solution of the matrix equation, ℓ , contains the Bessel function expansion coefficients of the gravity field perturbation. We thus introduce a square matrix H ℓ defined by where , ∈ [1, ]. Then the TGWE with radially varying DSG coefficient can be written in the symbolic matrix form where ℓ represents the source vector per degree along all radial grid points and ℓ is the solution vector containing the expansion coefficients. The matrix equation is solved via LU factorization. For dealiasing, we thus ignore the highest 10% of the Bessel coefficients. We use an equidistant radial mesh with up to = 1024 grid points, Gauss-Legendre points along latitude ( = 4096) for an accurate Legendre expansion (eq. 36) and solve for the first four odd spherical harmonic degrees (ℓ = 3, 5, 7, 9).
Having obtained the solution ℓ in Bessel space, eq. 39 yields the radial representation ℓ . The gravity moments are then simply given by If required, the density perturbation can be calculated by solving ∇ 2 = 4 . To test our method, we compare the solution for constant with results using the semi-analytical method introduced by . To verify our method with radius-dependent , we solve eq. 37 with an independent solver based on the shooting method. This is a standard tool for solving initial value problems and can readily be applied here, making use of a variable transformation to account for the asymptotic behaviour at the origin when applying this method.

MODELLING THE GRAVITY PERTURBATIONS
We start by exploring the challenges of modelling the odd gravity perturbations for a reference model that combines the interior model by Guillot et al. (2003) with the average flow introduced in sec. 3.2 and employs the TGWE method (eq. 37 or 41). Furthermore the dynamic density source is calculated via eq. 26 and thus we assume a hemispherically geostrophic flow, keep the equatorial step and apply the hyperbolic tangent radial decay profile for in accordance with eq. 34. Fig. 5 shows an attempt to model the observed gravity harmonics by varying the two parameters in , the decay depth ℎ and the decay width ℎ . A latitude-dependent attenuation in the form of eq. 28 or eq. 29 is not applied at this point. To assess the quality of the modelled gravity moments, we compare them individually for each degree ℓ to the observations rather than minimising an ℓ-independent, global cost function (Kaspi et al. 2018). Semi-transparent, coloured, regions indicate the parameter combinations where the respective gravity harmonic stays within 3 of the observations (Durante et al. 2020). There is a trade-off between both parameters since increasing either boosts the impact of deeper regions. Unfortunately, there is no parameter combination where all stripes overlap, i.e. no combination of ℎ and ℎ that would explain all four odd harmonics with one single decay function. Gravity harmonic 3 seems to be particularly problematic as it always maintains some distance from the other three, since to match it requires particularly deep sources.

Reference model and equatorial treatment
The grey background contours show at which radius the radial decay function would have decayed by three orders of magnitude. The SV constraint by Moore et al. (2019) suggests that this should happen somewhere between = 0.93 and = 0.96 . Modelling the observed 3 always requires deeper sources that are incompatible with this constraint. For the smooth decay functions suggested by Kaspi et al. (2018) or Kong et al. (2018) the three orders of magnitude drop lies far below = 0.90 (see also fig. 1).
Numerical simulations and theoretical consideration suggest that the flow remains geostrophic until a stably stratified layer, Lorentz forces or a combination of both, drastically quench the amplitude over a few hundred kilometres Wicht & Gastine 2020). We therefore restrict our analysis to ℎ = 500 km in the following, a value that represents the scale height of the electrical conductivity in the outer atmosphere of Jupiter quite well (French et al. 2012;Wicht et al. 2019) and will be further validated from the results of the TC model (see sec. 5.2). Fig. 6 illustrates how the gravity harmonics change when varying ℎ while keeping ℎ = 500 km fixed. For the reference model (Fig. ??, a), the harmonics 5 , 7 , and 9 all agree with the respective observations for ℎ values between 1500 km and 2300 km. However, 3 requires much deeper flows with ℎ ≈ 6500 km. Tab.1 gives the values for ℎ and the respective uncertainties for all considered models. We use the spherical harmonic degree ℓ as an index to denote the different values of ℎ ℓ required to explain the different observations. For a successful model, all four ℎ-values must agree within the errors. This is clearly not the case for our reference model. Only for degree ℓ = 7 and ℓ = 9 the decay depths match within the uncertainties at ℎ ≈ 2000 km. For reasonable values of ℎ between 2000 km and 4000 km the modelled gravity moment 3 is of the wrong sign, whereas 5 is much larger than the observed value in that range of ℎ.
Panel b) in fig. 6 illustrates that the situation improves when we follow the approach by Kaspi et al. (2018) and ignore the contribution due to the equatorial step in eq. 26. While the values of ℎ 5 , ℎ 7 and ℎ 9 remain nearly unchanged, ℎ 3 decreases by half to about 3200 km (see also tab. 1). This shows the large impact of the equatorial discontinuity, in particular on degree ℓ = 3. If in addition to this we assume the smoothly decaying radial function suggested by Kaspi et al. (2018), ignore the DSG term and utilise the perĳove surface flow model, we can largely reproduce their results and find a reasonably good agreement of the decay depths across the degrees.
Panel c) of fig. 6) shows the results when using the additional linear z-dependence suggested by Kong et al. (2018). This leads to shallower winds closer to the equator, a smooth vertical gradient and thus no equatorial step. Again, ℎ 3 decreases significantly, but not as much as for the model that ignores the equatorial step. The other harmonics are somewhat more affected, such that all ℎ ℓ increase by roughly 10% with respect to the reference model (tab. 1). Though this additional modification of the vertical flow structure avoids the equatorial discontinuity, it is hard to justify physically. Again we would have to also adopt the other model ingredients (TGWE model with polytropic interior, Cassini flow with capped spike, inverse Gauss decay profile) to reproduce the results by Kong et al. (2018) and match all harmonics.

TC model
Finally we explore our TC model where the flow amplitude outside of an assumed tangent cylinder is reduced by applying an additional attenuation function (eq. 30). The main parameters to examine are the latitude of the TC, , and the decay depth ℎ (fig. 7). Like in fig. 5 the transparent stripes show where individual gravity harmonics are in agreement with the observed values (Durante et al. 2020). The TC decay width and the radial decay width are fixed to = 2 • and = 500 km respectively. For small the solution is not affected, implying that the antisymmetric zonal flow between a latitude of ±15 • is irrelevant for the gravity signal. This changes once the TC angle is increased beyond = 15 • since we start reducing the amplitude of the dominant prograde jet at about = 21 • latitude. The required decay depth ℎ ℓ decreases for ℓ = 3 but increases for ℓ = 5, 7, and 9. At ℎ = 2975 km and = 20.9 • the model can finally explain all observed gravity harmonics (see also tab. 1). The corresponding surface flow of the TC model is shown in the inset of fig. 7 (blue profile) indicating that the dominant jet around 20 • is thinner than in the reference model (black) and reduced to a peak amplitude of 35 m/s. Also Kong et al. (2018) found it necessary to reduce the amplitude of the prograde jet.
Reducing the flow in an equatorial latitude band reduces the impact of the equatorial step on the density anomaly. It even vanishes when the TC is positioned at sufficiently large latitudes such that the flow is split into two independent parts. This must be consistent with the assumed radial decay function : a deeper-reaching flow (larger ℎ or ℎ ) requires a larger TC to separate the northern and southern hemisphere. Thus, a geometric relation between radial decay function and the width of the equatorial cut-off TC can be formulated. The dark curves in fig. 7 indicate the position of a TC attached to various depths; i.e. where the associated drop-off equals 10 −1 (light grey), 10 −2 (dark grey) or 10 −3 (black). Interestingly, the black line exactly hits the intersection point of the coloured isocontours. This means that the TC defined by reducing the flow amplitude in an equatorial band of ±20.9 • coincides with the TC defined by the 10 −3 -drop-off of the radial decay function . Thus this marks the point in the ( TC -ℎ)-parameter space, where the flow is just split into two separate hemispheric flows at the minimum possible TC .
Moreover, the other two parameters, ℎ and , can be constrained with this result. Higher values of do not lead to an intersection of the individual gravity solutions (colours in fig. 7), whereas larger or smaller values of ℎ shift the black line out of the intersection point.
Consequently, the TC gravity model not only explains all grav- ity measurements; it is also independent of the handling of the equatorial step since the geostrophic extension of the antisymmetric flows in either hemisphere do not reach the equatorial plane. Fig. 6 d) illustrates how the gravity harmonics change for = 20.9 • when increasing ℎ. The TC scenario most strongly affects 3 which now remains negative for all ℎ. However, the other harmonics are more significantly changed than by the different approaches to treat the equatorial discontinuity illustrated in fig. 6. The subsequent ℎ 3 is much smaller than in other models, whereas the other ℎ increase drastically. All ℎ ℓ are in agreement within the uncertainties at approximately ℎ = 2975 km (see also tab. 1).

High sensitivity of 3
Evidently the 3 gravity moment is highly sensitive to different models for the equatorial treatment. To understand this, we analyse the respective source -the dynamic density perturbation 3 -in more detail. Separating the first term in eq. 26 into two and transforming in Legendre space yields three source contributions for ℓ = 3: where F 3 and H 3 exhibit the Legendre-transformations of the integrated flow and the step contribution: The first part, 3, mainly represents the geostrophic part of the flow (where is constant), the second term, 3, is dominated by the decay region and the last, 3, shows the influence of the equatorial step. In fig. 8 the left panel shows the resulting profiles for the reference model with ℎ = 6541 km, i.e. where the gravity model reproduces the observed 3 . The right panel displays the TC model with ℎ = 2975 km. The radial functions F 3 and H 3 are plotted in the insets of the respective model. For the reference model, the first contribution (blue) is proportional to the radial derivative of and is therefore negative since F 3 is also negative. It dominates at shallow depths > 0.95 . The second source contribution (orange) is proportional to the radial derivative of the depth profile and therefore positive. It is only significant between 0.90 and 0.92 where the large ℎ-value positions the steep decrease in . The flow amplitude represented by F 3 is rather small at this depth hence the small positive bump in .
The third contribution, (green) is mostly positive (due to H 3 ) and thus has the wrong sign for explaining the negative 3 . For a successful model, the negative has to overcompensate and which only succeeds when reaches deep enough. Both and individually would yield much larger 3 amplitudes than required. Only the delicate balance between the three contribution yields the correct sign and amplitude. This explains the sensitivity of 3 to the reference model setup.
The source analysis also readily explains why ignoring in the hemispherical model (panel b) of fig. 6) allows one to explain 3 with a much smaller value of ℎ 3 . The linear -dependence assumed for panel c) of fig.6 only decreases the amplitude of and is therefore less effective.
For the TC model (right panel of fig 8) the equatorial step contribution, 3, is identical to zero validating the results from the previous chapter (see also fig. 6, d). The density anomaly entirely comprises of a shallow, negative and deep-rooted positive decay part . Evidently, is larger because F 3 is larger at the depth of the decay (close to 0.96 ).

Influence of dynamic self-gravity
We have mentioned above that some authors ignore the dynamic self-gravity or DSG-term (Kaspi et al. 2018;, which allows them to solve the TWE rather than the TGWE equation. We model TWE-based gravity moments by setting = 0 in eq. 17. Further, in order to isolate the importance of a realistic radial -profile we model the constant coefficient TGWE equation used by Zhang et al. (2015); Kong et al. (2018);  by calculating gravity moments with a constant = 2 / 2 .
Tab. 1 lists the decay depth ℎ ℓ that is required to to match the observed gravity signal. Ignoring the DSG term (TWE model, or = 0) shows a reduction of ℎ 3 by almost 1000 km (or 15%) and an increase of ℎ 5 by 170 km (10%). In general, the absolute difference shrinks with the degree ℓ. The model with = 2 / 2 shows an intermediate behaviour indicating that a more realistic radial profile of is an important part of the model.
Alternatively, we calculate the deviation of the gravity moments ℓ between different models for fixed degree ℓ and decay depth ℎ. We quantify the error in ℓ -error of neglecting the DSG term as a function of decay depth ℎ in terms of TWE ℓ The measure cTGWE ℓ , defined in analogy to eq. 48, quantifies the error for the constant coefficient TGWE. Fig. 9 shows how both errors change for our reference model when varying ℎ between 1000 km and 5000 km (we exclude larger ℎ values because here 2 changes sign and the definition in eq. 48 becomes useless). Except reference model TC model  for ℓ = 9 the errors depend little on ℎ in this range. For ℓ = 3 the error TWE 3 always exceeds 60%. For ℓ = 5 the error is still significant at about 20%. For ℓ = 7 and ℓ = 9, the error reaches about 10% for larger values of ℎ. These errors are even larger than quantified by Zhang et al. (2015) or  and show that the DSG term with a radius dependent DSG coefficient should not be neglected. The cTGWE approach, i.e. where is a constant, captures roughly half of the effect of the full TGWE approach. For example, for ℓ = 3 the error cTGWE 3 amounts to 40%, whereas TWE 3 exceeds 70% at ℎ = 3000 km We can conclude that the DSG term is a first order effect and should be taken into consideration in its radius dependent form for gravity inversion problems.

Sensitivity to surface zonal flow model
We have discussed above that the different zonal wind models published over the years likely are likely to reflect seasonal variations, which may be restricted to the very shallow atmosphere and therefore do not affect the gravity signal. This motivates us to use a mean flow averaged over the available profiles from different epochs (see sec. 3.2 and fig. 4). Kaspi et al. (2018) relied on the HST observations obtained during Juno's perĳove 9 pass (Tollefson et al. 2017), while Kong et al. (2018) based their analysis on the Cassini measurements (Porco et al. 2003). We implement both, the Perĳove and the Cassini surface flow model, in our standard setup and additionally apply the amplitude cap by Kong et al. (2018) in a model we call 'cap'. Tab. 1 lists the change in ℎ ℓ required when using either flow and shows that the largest influence is again found for degree ℓ = 3, where in particular the Cassini measurements require deeper sources. The Perĳove profile, on the other hand yields larger offsets in the higher degrees. The results of Kong et al. (2018) are based in a modified version of the Cassini profile (Porco et al. 2003), i.e. capping the large spike at 20 • -latitude to a maximum amplitude of spike = 43.2 m/s. This flow modification is a necessary part of their model in order to match the gravity moments. Tab. 1 show that the major effect is once more a reduction of ℎ 3 and a smaller increase for the higher degrees. Note once again that the successful reproduction of Jupiter's odd-degree gravity moments by Kong et al. (2018) is due to a combination of applying this particular flow modification, the additional linear -dependence, the smooth radial decay function, assuming a polytropic interior and solving the TGWE.
Generally, the impact of different surface flow profiles reaches 10% and is thus smaller than, for example, the effect of the equatorial step or neglecting the DSG.

Influence of interior model
Finally we also test the impact of using different interior structure models. We replace the model by Guillot et al. (2003) in our reference setup by the simple polytrope of index unity or the model by Nettelmann et al. (2012). Details of the models are outlined in sec. 3.1. The chosen interior model factors in the dynamic density source (eq. 26) as the background density and gravity, and , and in the DSG coefficient (eq. 13). Fig. 2 shows the relevant profiles of density, density gradient and . Tab. 1 shows that the impact on the ℎ ℓ is even smaller than the impact of the surface flow model, always remaining below 10%. This indicates that the odd-degree gravity moments are merely insensitive to details of the Jupiter's interior structure, such as properties of the core or the heavy element abundance.

DISCUSSION
For the first time, the precise measurements of Jupiter's gravity field by the Juno mission allowed quantifying the tiny undulations caused by the zonal winds. First attempts to invert the odd-degree gravity moments suggested a rather smooth decay with depth, reaching about 2 m/s or even 20 m/s at a radius of = 0.94 (Kaspi et al. 2018;Kong et al. 2018). This seems at odds with the magnetic observations, suggesting that the speed at this depth should not exceed the centimetre per second level (Moore et al. 2019). The gradual decay is also difficult to reconcile with the results from numerical simulations and theoretical considerations. They suggest the winds should remain geostrophic in the outer part of the atmosphere until they are abruptly quenched at some depth Gastine & Wicht 2021). In addition to the too smooth radial decay, the handling of the discontinuity at the equatorial plane that can arise when asymmetric surface winds are downward continued in z-direction (parallel to the rotation axis). Such a discontinuity is unphysical, but if it existed it would contribute substantially to the modelled gravity signal. To avoid this problem, Kong et al. (2018) therefore somewhat arbitrarily used an additional linear variation of the zonal with to eliminate the discontinuity, which is obviously at odds with a largely geostropic flow. Kaspi et al. (2018) chose an alternative approach and simply ignore the contribution originating from the discontinuity. Our analysis suggest that this was essential to successfully model the observed odd gravity harmonics but at the expense of physical consistency. More recently,  suggest a different radial profile that seems consistent with the magnetic observations consisting of a geostrophic outer envelope and a steep rapid decay at = 0.972 . They also adopt the arguable approach of ignoring the equatorial discontinuity from Kaspi et al. (2018). The differences in the findings might be caused by various different modelling assumptions regarding how the zonal flow relates to anomalies on the gravity potential, the treatment of the equatorial step, the applied surface zonal flow model and the interior model. Here we aim to quantify the impact of each of those. more faithfully. Our results suggest that using a specific individual profile rather than the mean significantly biases the analysis.
We also quantify the impact of different interior models and find that its influence is smaller than the equatorial treatment, the DSG term or a specific flow profile. Hence constraining the properties of the deep interior structure (Guillot et al. 2018), such as core mass or composition, on the basis of the odd gravity moments seems rather challenging as they are far more sensitive to other model assumptions, such as the equatorial treatment, considering the DSG term or assumptions about the structure of the zonal winds.
Here we suggest an alternative scenario that is motivated by the observation that the equatorially antisymmetric flow contributions are much weaker outside than inside the tangent cylinder (TC) in numerical simulations. The TC is the imaginary cylinder with a radial distance from the rotation axis that is defined by the depth at the equator where deviations from a barotropic state and/or Lorentz forces enforce a sharp drop of the zonal flow. This seems consistent with the simple fact that equatorially antisymmetric flows are necessarily ageostrophic outside the TC, while there is no such conflict inside the TC. In our model we therefore considerably damp the antisymmetric flow outside the TC, smoothing the transition with a tangent hyperbolic function. This relieves us from the problem of the equatorial discontinuity without having to assume an unrealistic -dependence of the flow. We thus favour a model where the antisymmetric surface flow in an equatorial latitude band of ±21 • is a shallow phenomena, not contributing to the gravity signal. This, in turn, requires strong baroclinicity in the first few hundred kilometers below the cloud level generated from latitudinal gradients in temperature or composition. Interestingly, such pronounced equatorial variations in the nadir brightness temperature and the ammonia mixing ratio have been detected by the Juno mission (Bolton et al. 2017).
We use the four significant odd moments, 3 to 9 to individually constrain the vertical and latitudinal structure of the zonal winds. We can match the modelled gravity moments to the observations within the 3 uncertainties (Durante et al. 2020) when using a radial decay profile that models a geostrophic flow inside TC for > 0.958 but then decays rapidly with depth (see also tab. 2). Outside TC the flow must be shallow. Our preferred decay profile suggests that the winds are 50% deeper than suggested by  and hence the geostrophic part reaches to substantially larger pressures (830 kbar rather than 300 kbar). At = 0.935 the profile reaches finally 10 −3 what is in agreement with the constraints from Moore et al. (2019). This depth equals also to the lower boundary of a TC touching the surface at ±21 • implying that our preferred wind structure exhibits a splitting of the gravity perturbation into an independent northern and southern hemisphere, thus the contribution from the equatorial step naturally vanishes.
In order to explain the odd gravity moments with a physically consistent handling of the flow near the equatorial plane, we find it necessary to exclude part of the observed cloud level zonal winds from the analysis by assuming that it is restricted to shallow depth. Similarly, Kong et al. (2018) had to cap the amplitude of the strong jet centred at a latitude of = 21 • N to fit the gravity moments in their model. While it is physically reasonable that deep antisymmetric winds do not exist outside the tangent cylinder (defined with the depth at the equator where deep winds drop sharply in amplitude), there is obviously no guarantee that the antisymmetric cloud-level winds inside the tangent cylinder are fully representative of deep flow and do not also comprise a shallow component. This uncertainty is probably the severest limitation for utilising the gravity anomalies to infer the precise depth extent of the zonal flow. A better understanding of the atmospheric dynamics at different latitudes in the top few hundred kilometres below the clouds is necessary to overcome this limitation.