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

J. Wicht  Max Planck Institute for Solar System Research, Justus-von-Liebig-
Weg 3, 37077 Göttingen, Germany  W. Dietrich  Max Planck Institute for Solar
System Research, Justus-von-Liebig-Weg 3, 37077 Göttingen, Germany  P. Wulff
Max Planck Institute for Solar System Research, Justus-von-Liebig-Weg 3, 37077
Göttingen, Germany  U. R. Christensen  Max Planck Institute for Solar System
Research, Justus-von-Liebig-Weg 3, 37077 Göttingen, Germany

###### Abstract

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 $\ell$. For degrees $\ell=2$ to about $\ell=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 $\ell=5$
to roughly $\ell=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.

###### Acknowledgements

This work was supported by the German Research Foundation (DFG) in the
framework of the special priority programs ’Exploring the Diversity of
Extrasolar Planets’ (SPP 1992).

##  1 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 back-reaction 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 $\rho$. The gravity potential

$\varPsi(\mathbf{r})=-\frac{GM}{r}\;\left[1-\sum_{\ell=2}^{\inf}\,J_{\ell}\;\left(\frac{R}{r}\right)^{\ell}\;P_{\ell}(\theta)\right]$
(1)

thus acquires equatorially symmetric contributions of even degree $\ell=2n$
with $n=1,2,3,...$. Here $G$ is the gravity constant, $M$ the planetary mass,
$R$ the planetary radius, $\theta$ the colatitude, and $P_{\ell}$ a Schmitt-
normalized Legendre Polynomial of degree $\ell$. The gravity harmonics
$J_{\ell}$ are given by the volume integral

$J_{\ell}=\frac{2\pi}{MR^{\ell}}\;\int\,d\,V\;{r}^{\ell}\;\rho(r,\theta)\;P_{\ell}(\theta)$
(2)

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=\Omega^{2}/(G\rho)$, where $\Omega$ 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

$2\Omega\;\frac{\partial\overline{\rho}\,U_{\phi}}{\partial
z}=\hat{\mathbf{\phi}}\cdot{\mathbf{\nabla}}\times\;\left({\rho}^{\prime}\,{\mathbf{\nabla}}\overline{\varPsi}_{e}+{\varPsi}^{\prime}\,{\mathbf{\nabla}}\overline{\rho}\;\right)\;\;,$
(3)

where $z$ is the distance to the equatorial plane, $\overline{\varPsi}_{e}$
the effective background potential, $\overline{\rho}$ the background density,
${\rho}^{\prime}$ the density perturbation and ${\varPsi}^{\prime}$ 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 ${\varPsi}^{\prime}$ 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 ${\rho}^{\prime}$ (Zhang et al., 2015; Kong et al., 2016, 2017,
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.

##  2 From Navier-Stokes Equation to
Inhomogeneous Helmholtz Equation

The link between the dynamics and gravity is provided by the Navier-Stokes
equation

$\rho\left(\frac{\partial}{\partial
t}+\mathbf{u}\cdot{\mathbf{\nabla}}\right)\,\mathbf{u}+2\varOmega\rho\;\hat{\mathbf{z}}\times\mathbf{u}=-{\mathbf{\nabla}}p+\rho\,\mathbf{g}_{e}+\mathbf{j}\times\mathbf{B}+\nu\,{\mathbf{\nabla}}\cdot\mathcal{S}\;\;,$
(4)

where $\mathbf{u}$ is velocity, $\hat{\mathbf{z}}$ the unit vector in the
direction of the rotation axis, $p$ the pressure, $\mathbf{j}$ the electric
current, $\mathbf{B}$ the magnetic field, $\nu$ the kinematic viscosity, and
$\mathcal{S}$ the traceless rate-of-strain tensor for constant kinematic
viscosity:

$\mathcal{S}=\rho\left(\frac{\partial u_{i}}{\partial x_{j}}+\frac{\partial
u_{j}}{\partial
x_{i}}-\frac{2}{3}\delta_{ij}{\mathbf{\nabla}}\cdot\mathbf{u}\right)\;\;.$ (5)

The effective gravity $\mathbf{g}_{e}$ can be expressed by an effective
gravity potential,

$\mathbf{g}_{e}=-{\mathbf{\nabla}}\varPsi_{e}=-{\mathbf{\nabla}}\left(\varPsi+\varPsi_{\Omega}\right)\;\;,$
(6)

which is the sum of the gravity potential obeying the Poisson equation

$\nabla^{2}\varPsi=4\pi G\;\rho$ (7)

and the centrifugal potential

$\varPsi_{\Omega}=-\frac{1}{2}\;\Omega^{2}s^{2}\;\;,$ (8)

with $s=r\sin{\theta}$ being the distance to the rotation axis.

The zero order force balance is given by the hydrostatic equilibrium with
vanishing flow and magnetic field:

${\mathbf{\nabla}}\overline{p}=-\overline{\rho}\;{\mathbf{\nabla}}\overline{\varPsi}_{e}\;\;,$
(9) $\nabla^{2}\overline{\varPsi}_{e}=\left(4\pi
G\;\overline{\rho}-\Omega^{2}\right)\;\;.$ (10)

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

$\overline{\rho}\left(\frac{\partial}{\partial
t}+\mathbf{u}\cdot{\mathbf{\nabla}}\right)\,\mathbf{u}+2\varOmega\overline{\rho}\;\hat{\mathbf{z}}\times\mathbf{u}=-{\mathbf{\nabla}}{p}^{\prime}-\overline{\rho}\,{\mathbf{\nabla}}{\varPsi}^{\prime}_{e}{\rho}^{\prime}\,{\mathbf{\nabla}}\overline{\varPsi}+\mathbf{j}\times\mathbf{B}+\nu\,{\mathbf{\nabla}}\cdot\mathcal{S}\;\;,$
(11) $\nabla^{2}{\varPsi}^{\prime}=4\pi G\;{\rho}^{\prime}\;\;.$ (12)

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:

$\overline{\rho}{\mathbf{\nabla}}{\varPsi}^{\prime}={\mathbf{\nabla}}(\overline{\rho}\,{\varPsi}^{\prime})\;-\;{\varPsi}^{\prime}\,{\mathbf{\nabla}}\overline{\rho}\;\;.$
(13)

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:

$2\Omega\;\frac{\partial}{\partial
z}\left(\overline{\rho}\,U_{\phi}\right)=\hat{\mathbf{\phi}}\cdot\left({\mathbf{\nabla}}\times\left[\;{\rho}^{\prime}\,{\mathbf{\nabla}}\overline{\varPsi}_{e}-{\varPsi}^{\prime}\,{\mathbf{\nabla}}\overline{\rho}\;\right]\right)\;\;.$
(14)

The next step is to assume that $\psi_{\Omega}$ can be neglected in comparison
to the background gravity contribution $\overline{\varPsi}$, as discussed in
the introduction. The background state then becomes spherically symmetric and
equation (14) simplifies to

$2\Omega\;\frac{\partial}{\partial
z}\left(\overline{\rho}\,U_{\phi}\right)=\frac{1}{r}\,\frac{\partial}{\partial\theta}\left(\;{\rho}^{\prime}\,\frac{\partial}{\partial
r}\overline{\varPsi}-{\varPsi}^{\prime}\,\frac{\partial}{\partial
r}\overline{\rho}\;\right)\;\;.$ (15)

This is the thermo-gravitational wind equation (TGWE) solved for a given
$U_{\phi}$ 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, $\overline{\rho}{\mathbf{\nabla}}{\varPsi}^{\prime}$, or
more precisely its non-conservative contribution.

Integrating equation (15) in latitude, dividing by background gravity
$\overline{\mathbf{g}}=-\partial\overline{\varPsi}/\partial 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:

$\left(\nabla^{2}+\mu\right){\varPsi}^{\prime}=4\pi\,G\;\rho^{U}\;\;,$ (16)

with

$\mu(r)=4\pi\,G\,{\mathbf{\nabla}}\overline{\rho}\,\big{/}\,\overline{g}\;\;,$
(17)

and the dynamic density perturbation

$\rho^{U}(r,\theta)=\frac{2\Omega
r}{\overline{g}}\;\int_{0}^{\theta}\,d\hat{\theta}\;\frac{\partial}{\partial
z}\,\left(\overline{\rho}\,U_{\phi}\right)\;\;,$ (18)

as a source term. Note that $\rho^{U}$ is an auxiliary variable different from
${\rho}^{\prime}$. We will refer to $\mu(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
${\varPsi}^{\prime}$ vanishes at $r=0$. Outside of the source, the solutions
must obey

$\nabla^{2}{\varPsi}^{\prime}=0\;\;.$ (19)

A respective matching condition at the outer radius $R$ yields the second
boundary condition that we provide further below.

Because $\rho^{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
$\rho^{U}={\rho}^{\prime}$ 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 ${\rho}^{\prime}$ derived by Zhang et al. (2015) and the Possion
equation for ${\varPsi}^{\prime}$ is then solved in a second step. Their
solution is cumbersome and numerically time-consuming. We avoid this
complication by directly solving the inhomogeneous Helmholtz-type equation
(16) to obtain $\psi^{\prime}$. The true density perturbation can be recovered
by

${\rho}^{\prime}=\rho^{U}-\frac{\mu}{4\pi G}\;{\varPsi}^{\prime}\;\;,$ (20)

which is obtained from equation (12) and equation (16). We note that
$\rho^{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 $\mu$? 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

$\overline{\rho}(r)=\overline{\rho}_{c}\,\frac{\sin{\chi}}{\chi}\;\;,$ (21)

where $\rho_{c}$ is the density at $r=0$, and $\chi$ a rescaled radius:

$\chi=\pi\;\frac{r}{R}\,\frac{\rho_{c}-\rho(R)}{\rho_{c}}\;\;.$ (22)

The gravity profile is then

$\overline{g}(r)=-4\pi\,G\;\frac{\rho_{c}^{2}}{\rho_{c}-\rho(R)}\;\frac{R}{\pi}\;\frac{\chi\cos{\chi}-\sin{\chi}}{\chi^{2}}$
(23)

and the DSG coefficient becomes constant:

$\mu(r)=\frac{\pi^{2}}{R^{2}}\;\left(\frac{\rho_{c}-\rho(R)}{\rho_{c}}\right)^{2}\approx\frac{\pi^{2}}{R^{2}}\;\;.$
(24)

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:

${\mathbf{\nabla}}\overline{\rho}=\beta_{S}\,\overline{\rho}\;{\mathbf{\nabla}}\overline{p}\;\;,$
(25)

with

$\beta_{S}=\frac{1}{\overline{\rho}}\left(\frac{\partial\rho}{\partial
p}\right)_{S}$ (26)

being the compressibility at constant entropy. Combining equation (25) and
equation (9) shows that the gradient in the background density is given by

$\frac{\partial}{\partial
r}\overline{\rho}=\beta_{S}\,\overline{\rho}^{2}\,\overline{g}\;\;.$ (27)

The DSG coefficient is thus given by

$\mu(r)=4\pi G\;\beta_{S}\,\overline{\rho}^{2}\;\;.$ (28)

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 $\mu(r)$
variations remain remarkable small. In the lower layer $r<0.25R$, $\mu(r)$ is
nearly constant and close to $\pi^{2}/R^{2}$. In the outer envelope
$r>0.85\,R$, $\mu$ becomes more variable, reaching amplitudes $40$% larger
than $\pi^{2}/R^{2}$. A constant $\mu$ 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.

Figure 1: Panel a) shows pressure versus density (solid line) for the Jupiter
model by Nettelmann et al. (2012) and French et al. (2012) and a polytrope of
index unity (dashed line). The double logarithmic plot highlights that this
polytrope, i.e. $p\sim\rho^{2}$, provides a decent approximation. The Jupiter
model by Nettelmann et al. (2012) and French et al. (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 $\mu(r)\,R^{2}$ (solid line) suggested by the ab-initio data
points by French et al. (2012) (circles) with the constant value $\pi^{2}$
expected for the polytrope (dashed line).

##  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 $\mu$. The involved methods represent textbook knowledge, but
their application to the specific gravity problem is new however, we
nevertheless discuss them in some detail.

###  3.1 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 $\varGamma$ is defined by

$\nabla^{2}\varGamma(\mathbf{r},\tilde{\mathbf{r}})=\delta(\mathbf{r}-\tilde{\mathbf{r}})\;\;,$
(29)

where vectors $\mathbf{r}$ and $\tilde{\mathbf{r}}$ 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

$\varPsi(\mathbf{r})=4\pi\,G\;\int\,d\,\tilde{V}\;\varGamma(\mathbf{r},\tilde{\mathbf{r}})\,\rho(\tilde{\mathbf{r}})\;\;,$
(30)

where

$\int\,d\,\tilde{V}=\int_{0}^{R}\,d\,\tilde{r}\;{\tilde{r}}^{2}\;\int_{0}^{\pi}\,d\,\tilde{\phi}\;\int_{0}^{2\pi}\,d\,\tilde{\theta}\sin{\tilde{\theta}}$
(31)

denotes the integration over the spherical volume.

The classical Green’s function for the Poisson problem is given by

$\varGamma(\mathbf{r},\tilde{\mathbf{r}})=-1\big{/}\,\left(4\pi\left|\mathbf{r}-\tilde{\mathbf{r}}\right|\right)\;\;,$
(32)

but of more practical use is the representation where $\varGamma$ 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:

$\nabla^{2}\;f(r)\;P_{\ell}(\theta)\;=\;\left(\;\frac{\partial^{2}}{\partial
r^{2}}\,+\,\frac{2}{r}\frac{\partial}{\partial
r}\,-\,\frac{\ell(\ell+1)}{r^{2}}\;\right)\;f(r)\;P_{\ell}(\theta)\;\;.$ (33)

The Schmitt normalization assumed here means that

$\int_{0}^{\pi}\,d\theta\sin\theta\;P_{\ell}(\theta)\;P_{\ell^{\prime}}(\theta)\;=\;\frac{2}{2\ell+1}\;\delta_{\ell\ell^{\prime}}\;\;.$
(34)

The two possibilities for the radial function are $f_{\ell}(r)=r^{\ell}$ and
$f_{\ell}(r)=r^{-(\ell+1)}$. The expanded Green’s function then reads

$\varGamma(\mathbf{r},\tilde{\mathbf{r}})=-\frac{1}{4\pi}\sum_{\ell=0}^{\infty}\;\frac{{r_{<}}^{\ell}}{r_{>}^{\ell+1}}\;P_{\ell}(\theta)\;P_{\ell}(\tilde{\theta})\;\;,$
(35)

where $r_{>}$ ($r_{<}$) denotes that larger (smaller) of the two radii $r$ and
$\tilde{r}$. The matching condition to the field for $r>R$ reduces to the
mixed boundary condition

$\frac{\partial}{\partial
r}\;f_{\ell}(r)=-\frac{(\ell+1)}{R}\;f_{\ell}(r)\;\;,$ (36)

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

$\varPsi(\mathbf{r})=\sum_{\ell=0}^{\infty}\;\varPsi_{\ell}\;\left(\frac{R}{r}\right)^{\ell+1}\;P_{\ell}(\theta)\;\;,$
(37)

with the expansion coefficients

$\varPsi_{\ell}=-\frac{G}{4\pi
R}\;\int\,d\tilde{V}\;\left(\frac{\tilde{r}}{R}\right)^{\ell}\;\rho(\tilde{\mathbf{r}})\,P_{\ell}(\tilde{\theta})\;\;.$
(38)

This is equivalent to the differently normalized classical expansion equation
(1) and equation (2).

The same solution applies to ${\varPsi}^{\prime}$ when replacing $\rho$ by
${\rho}^{\prime}$. Should the impact of DSG $\mu$ be negligible, we could
simply use ${\rho}^{\prime}\approx\rho^{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).

###  3.2 Solving the Inhomogeneous Helmholtz equation

For constant $\mu(r)=K^{2}$, the modified potential field equation becomes an
inhomogeneous Helmholtz equation

$\left(\nabla^{2}+K^{2}\right)\;{\varPsi}^{\prime}=4\pi\,G\;\rho^{U}\;\;.$
(39)

The respective Green’s function is now defined by

$\left(\,\nabla^{2}+K^{2}\,\right)\varGamma=\delta(\mathbf{r}-\tilde{\mathbf{r}})\;\;.$
(40)

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_{\ell}(r)P_{\ell}(\theta)$ of the Laplace operator where the $f_{\ell}(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

$\left(\frac{\partial^{2}}{\partial r^{2}}+\frac{2}{r}\frac{\partial}{\partial
r}-\frac{\ell(\ell+1)}{r^{2}}\;+\;1\right)j_{\ell}(r)=0\;\;.$ (41)

We only use the spherical Bessel functions of the first kind, $j_{\ell}$, with
$\ell>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:

$\nabla^{2}\;j_{\ell}(k_{\ell n}r)\;P_{\ell}(\theta)=\lambda\;j_{\ell}(k_{\ell
n}r)\;P_{\ell}(\phi)\;\;,$ (42)

with eigenvalues

$\lambda=-k_{\ell n}^{2}\;\;.$ (43)

The different $k_{\ell n}$ are chosen so that $j_{\ell}(k_{\ell n}R)$ fulfills
the boundary condition (36). Because of recurrence relation (88) (see App. C),
this condition reduces to

$j_{\ell-1}(k_{\ell n}R)=0\;\;,$ (44)

which means that the $k_{\ell n}$ are the roots of $j_{\ell-1}(x)$ divided by
the outer boundary radius $R$. We start the numbering at the smallest root
larger than zero so that $0<k_{\ell 1}<k_{\ell 2}<k_{\ell 3}<...$. Panel (a)
of Fig. 2 illustrates the spherical Bessel functions $j_{\ell}$ for different
degrees $\ell$. Table 1 list the first five roots for $\ell\leq 5$.

$\ell$ / $n$ | 1 | 2 | 3 | 4 | 5
---|---|---|---|---|---
1 | 1 | 2 | 3 | 4 | 5
2 | 1.4303 | 2.4590 | 3.4709 | 4.4774 | 5.4815
3 | 1.8346 | 2.8950 | 3.9225 | 4.9384 | 5.9489
4 | 2.2243 | 3.3159 | 4.3602 | 5.3870 | 6.4050
5 | 2.6046 | 3.7258 | 4.7873 | 5.8255 | 6.8518
Table 1: List of $k_{\ell n}R/\pi$. The $k_{\ell n}R$ are the roots of
$j_{\ell-1}$.

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_{\ell n}\;N_{\ell n^{\prime}}\;\int\,dr\;r^{2}\;j_{\ell}(k_{\ell
n}r)\;j_{\ell}(k_{\ell n^{\prime}}r)=\delta_{n,n^{\prime}}\;\;,$ (45)

where the $N_{\ell n}$ are normalization constants derived analytically in
Sect. B:

$N_{\ell n}=\left(\frac{2}{R^{3}j_{\ell}^{2}(k_{\ell n}R)}\right)^{1/2}\;\;.$
(46)

Panel (b) of Fig. 2 shows the first five normalized functions,

$j^{\star}_{\ell n}(r)=N_{\ell n}\;j_{\ell}(k_{\ell n}r)\;\;,$ (47)

for $\ell=2$.

Figure 2:  Panel a) shows the first five spherical Bessel functions of the
first kind. Panel b) shows the first orthonormal normalized functions
$j^{\star}_{\ell n}$ for degree $\ell=2$.

We can now expand the potential field perturbation in Legendre polynomials and
the new orthonormal radial functions:

${\varPsi}^{\prime}(\mathbf{r})=\sum_{n=1}^{\infty}\,\sum_{\ell=1}^{\infty}\;{\varPsi}^{\prime}_{\ell
n}\;j^{\star}_{\ell n}(r)\,P_{\ell}(\theta)\;\;.$ (48)

Using this expansion in equation (39), multiplying with the ansatz functions
$j^{\star}_{\ell n}(\tilde{r})P_{\ell}(\tilde{\theta})$ and integrating over
the volume yields a spectral equation for the expansion coefficients:

$\frac{4\pi\,\left(k_{\ell
n}^{2}-K^{2}\right)}{(2\ell+1)}\;{\varPsi}^{\prime}_{\ell
n}=-4\pi\;G\;\int\,d\,V\;\rho^{U}(\mathbf{r})\;j^{\star}_{\ell
n}\;(\tilde{r})\,P_{\ell}(\tilde{\theta})\;\;.$ (49)

The coefficients are thus simply given by

${\varPsi}^{\prime}_{\ell n}=-\frac{G\,(2\ell+1)}{k_{\ell
n}^{2}-K^{2}}\;\int\,d\,\tilde{V}\;\rho^{U}(\tilde{\mathbf{r}})\;j^{\star}_{\ell
n}(\tilde{r})\;P_{\ell}(\tilde{\theta})\;\;.$ (50)

A comparison with equation (35) shows that the Green’s function for the
inhomogeneous Helmholz equation is then

$\varGamma(\mathbf{r},\tilde{\mathbf{r}})=-\frac{1}{4\pi}\sum_{n=1}^{\infty}\,\sum_{\ell=1}^{\infty}\;\frac{(2\ell+1)}{k_{\ell
n}^{2}-K^{2}}\;j^{\star}_{\ell n}(r)\;j^{\star}_{\ell
n}(\tilde{r})\;P_{\ell}(\theta)P_{\ell}(\tilde{\theta})\;\;.$ (51)

The potential field for $r>R$ has to decay like $\left(R/r\right)^{\ell+1}$.
The respective solution is thus given by

${\varPsi}^{\prime}(\mathbf{r})=\sum_{\ell=1}^{\infty}\;{\varPsi}^{\prime}_{\ell}(R)\;\left(\frac{R}{r}\right)^{\ell+1}\;P_{\ell}(\theta)\;\;,$
(52)

with

${\varPsi}^{\prime}_{\ell}(R)=\sum_{n=1}^{\infty}\;{\varPsi}^{\prime}_{\ell
n}\;j^{\star}_{\ell n}(R)\;\;.$ (53)

As expected, this solution is identical to the classical result (37) for
$K^{2}=0$. We show this analytically in Sect. D.

##  4 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:

$\rho^{U}=j^{\star}_{\ell n}(r)\;P_{\ell}(\theta)\;\;.$ (54)

Only the respective potential field coefficient thus has to be considered and
the solution for $r<R$ is

${\varPsi}^{\prime}(\mathbf{r})=-\frac{G\,(2\ell+1)}{k_{\ell
n}^{2}-K^{2}}\;j^{\star}_{\ell n}(r)\;P_{\ell}(\theta)\;\;.$ (55)

Solving for a more general source thus boils down to the question: How well
can $\rho^{U}$ be expanded in the ansatz functions?

A special situation arises when $K^{2}=k_{\ell n}^{2}$. For the polytropic
density distribution with polytropic index unity, this happens for $\ell=1$
and $n=1$ where $K=k_{1,1}=\pi/R$. The two non-conservative buoyancy terms
then cancel exactly,

${\rho}^{\prime}{\mathbf{\nabla}}\overline{\varPsi}-{\varPsi}^{\prime}{\mathbf{\nabla}}\overline{\rho}=0\;\;,$
(56)

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 $\rho^{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 ${\varPsi}^{\prime}_{1,1}$ contributions. Since the $\ell=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

$\rho^{U}=\frac{2\Omega}{\overline{g}}\;\left(\overline{\rho}\,\sin{\theta}\,U_{\phi}+r\frac{\partial\overline{\rho}}{\partial
r}\;\int_{0}^{\pi}\,d\theta\;\cos\theta\,U_{\phi}\right.+\\\
\left.r\,\overline{\rho}\;\int_{0}^{\pi}\,d\theta\;\cos\theta\,\frac{\partial
U_{\phi}}{\partial r}\;\right)\;\;.$ (57)

While latitude-dependence is this purely determined by the zonal flow,
$\overline{\rho}$, $U_{\phi}$ and their radial derivatives influence the
radial profile of $\rho^{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 $\rho^{U}$:

$\rho^{U\\!N}_{\ell}(r)=\sum_{n=1}^{N}\;\rho^{U}_{\ell n}\;j^{\star}_{\ell
n}(r)\;\;,$ (58)

with

$\rho^{U}_{\ell
n}=\int_{0}^{R}\,d\,r\;r^{2}\;\rho^{U}_{\ell}(r)\;j^{\star}_{\ell n}(r)\;\;,$
(59)

and

$\rho^{U}_{\ell}(r)=\int_{0}^{\pi}\,d\,\theta\;\sin{\theta}\;\rho^{U}(r,\theta)\;P_{\ell}(\theta)\;\;.$
(60)

The quality of the representation is quantified by the misfit

$D(N)=\frac{\int_{0}^{r_{o}}\,d\,r\;r^{2}\;\left[\,\rho^{U\\!N}_{\ell}(r)-\rho^{U}_{\ell}(r)\,\right]^{2}}{\int_{0}^{r_{o}}\,d\,r\;r^{2}\;{\rho^{U}_{\ell}}^{2}(r)}\;\;.$
(61)

N | h=0.143 | h=1.143
---|---|---
| TWE | TGWE | TWE | TWGE
10 | $3.133\mbox{$\times 10^{-5}$}$ | $4.961\mbox{$\times 10^{-5}$}$ | $0.8419\mbox{$\times 10^{-4}$}$ | $1.489\mbox{$\times 10^{-4}$}$
20 | $3.165\mbox{$\times 10^{-5}$}$ | $4.992\mbox{$\times 10^{-5}$}$ | $0.8433\mbox{$\times 10^{-4}$}$ | $1.491\mbox{$\times 10^{-4}$}$
40 | $3.169\mbox{$\times 10^{-5}$}$ | $4.997\mbox{$\times 10^{-5}$}$ | $0.8435\mbox{$\times 10^{-4}$}$ | $1.491\mbox{$\times 10^{-4}$}$
60 | $3.170\mbox{$\times 10^{-5}$}$ | $4.998\mbox{$\times 10^{-5}$}$ | $0.8435\mbox{$\times 10^{-4}$}$ | $1.491\mbox{$\times 10^{-4}$}$
100 | $3.170\mbox{$\times 10^{-5}$}$ | $4.998\mbox{$\times 10^{-5}$}$ | $0.8435\mbox{$\times 10^{-4}$}$ | $1.491\mbox{$\times 10^{-4}$}$
Z2015 | $3.17\mbox{$\times 10^{-5}$}$ | $5.00\mbox{$\times 10^{-5}$}$ | $0.874\mbox{$\times 10^{-4}$}$ | $1.553\mbox{$\times 10^{-4}$}$
Table 2: Gravity harmonic $J_{2}$ for the equatorially symmetric test case
suggested by 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).

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

$U_{\phi}=U_{0}\;f_{1}(r)\;\sin^{2}{\theta}$ (62)

with amplitude $U_{0}=R\Omega/100$ and radial dependence

$f_{1}(r)=\left(\frac{r}{R}\right)^{2}\;\exp{\left(-\frac{1}{h}\frac{R-r}{R}\right)}\;\;.$
(63)

Jupiter values used to define flow and gravity are $R=6.9894\mbox{$\times
10^{7}$}\,$m, $\Omega=1.759\mbox{$\times 10^{-4}$}\,$s-1, and
$M=1.898\mbox{$\times 10^{27}$}\,$kg. Two relative decay scale heights
$h=0.143$ and $h=1.143$ are explored. The flow yields $\ell=0$ and $\ell=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).

Figure 3: Expansion of the function $f_{1}(r)$ with $h=1.143$ into the
$j^{\star}_{\ell n}$ for $\ell=2$. Panel a) compares the normalized function
with representations for truncations $N=10$, $40$, and $100$. Panel b) shows
the same in a logarithmic plot. Panel c) shows the spectrum for $N=101$ and
panel d) the misfit $D(N)$.

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^{\star}_{\ell 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\rightarrow 0$.

As a second example we explore the function

$f_{2}(r)=r^{\ell}$ (64)

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 $\ell=3$. Panels a)
and b) once more illustrate the difficulties of representing the function at
the boundaries.

Figure 4: Expansion of the function $f_{2}=r^{3}$ for $R=10$ into the
$j^{\star}_{\ell n}$ for $\ell=3$. Panel a) compares the normalized function
with representations for truncations $N=10$, $40$, and $100$. Panel b) shows
the same in a logarithmic plot. Panel c) shows the spectrum for $N=101$ and
panel d) the misfit $D(N)$.

The last example is the radial function

$f_{3}(r)=\frac{\overline{\rho}\,r}{\overline{g}}\;\frac{\partial
U_{\phi}(r)}{\partial r}$ (65)

that determines the radial dependence of one term in $\rho^{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:

$U_{\phi}(r)=\left\\{\begin{array}[]{ll}\exp{\left(\frac{1}{h}\,\frac{d^{2}}{D^{2}-d^{2}}\right)}&\mbox{for}\;d\leq
D\\\ 0&\mbox{for}\;d>D\end{array}\right.\;\;.$ (66)

where $d=R-r$ is the depth, $D=0.15\,R$ is the maximum depth of $\rho^{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^{\star}_{\ell 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^{\star}_{\ell n}$.

Figure 5: Same as Fig. 4 but for function $f_{3}(r)$. The $j^{\star}_{\ell n}$
for $\ell=3$ have been used.

##  5 Relative Importance of Dynamic Self Gravity

The analytical solution shows that the impact of the DSG simply depends on the
ratio $k_{\ell n}^{2}/K^{2}$. The relative importance of $K^{2}$ in the
inhomogenous Helmholtz equation for a given spherical harmonic degree $\ell$
and radial index $n$ can be quantified by

$S_{\ell n}=\frac{\left(k_{\ell n}^{2}-K^{2}\right)^{-1}}{k_{\ell
n}^{-2}}\;-\;1=\frac{1}{k_{\ell n}^{2}/K^{2}-1}\;\;.$ (67)

Table 3 lists $S_{\ell n}$ for spherical harmonic degrees up to $\ell=30$ and
$n$ up to $5$, assuming $K=\pi$. The values indicate that the DSG should be
considered a first order effect for $\ell\leq 4$, reaches the $10$% level at
$\ell=5$ or $\ell=6$ and amounts to only about $1$% for $\ell\geq 20$.

$\ell$ / $n$ | 1 | 2 | 3 | 4 | 5
---|---|---|---|---|---
1 | — | $3.3\mbox{$\times 10^{-1}$}$ | $1.2\mbox{$\times 10^{-1}$}$ | $6.7\mbox{$\times 10^{-2}$}$ | $4.2\mbox{$\times 10^{-2}$}$
2 | $9.6\mbox{$\times 10^{-1}$}$ | $2.0\mbox{$\times 10^{-1}$}$ | $9.1\mbox{$\times 10^{-2}$}$ | $5.3\mbox{$\times 10^{-2}$}$ | $3.4\mbox{$\times 10^{-2}$}$
3 | $4.2\mbox{$\times 10^{-1}$}$ | $1.4\mbox{$\times 10^{-1}$}$ | $7.0\mbox{$\times 10^{-2}$}$ | $4.3\mbox{$\times 10^{-2}$}$ | $2.9\mbox{$\times 10^{-2}$}$
4 | $2.5\mbox{$\times 10^{-1}$}$ | $1.0\mbox{$\times 10^{-2}$}$ | $5.6\mbox{$\times 10^{-2}$}$ | $3.6\mbox{$\times 10^{-2}$}$ | $2.5\mbox{$\times 10^{-2}$}$
5 | $1.7\mbox{$\times 10^{-1}$}$ | $7.8\mbox{$\times 10^{-2}$}$ | $4.6\mbox{$\times 10^{-2}$}$ | $3.0\mbox{$\times 10^{-2}$}$ | $2.2\mbox{$\times 10^{-2}$}$
6 | $1.3\mbox{$\times 10^{-1}$}$ | $6.2\mbox{$\times 10^{-2}$}$ | $3.8\mbox{$\times 10^{-2}$}$ | $2.6\mbox{$\times 10^{-2}$}$ | $1.9\mbox{$\times 10^{-2}$}$
8 | $7.8\mbox{$\times 10^{-2}$}$ | $4.3\mbox{$\times 10^{-2}$}$ | $2.8\mbox{$\times 10^{-2}$}$ | $2.0\mbox{$\times 10^{-2}$}$ | $1.5\mbox{$\times 10^{-2}$}$
10 | $5.3\mbox{$\times 10^{-2}$}$ | $3.2\mbox{$\times 10^{-2}$}$ | $2.2\mbox{$\times 10^{-2}$}$ | $1.6\mbox{$\times 10^{-2}$}$ | $1.3\mbox{$\times 10^{-2}$}$
14 | $3.0\mbox{$\times 10^{-2}$}$ | $2.0\mbox{$\times 10^{-2}$}$ | $1.4\mbox{$\times 10^{-2}$}$ | $1.1\mbox{$\times 10^{-2}$}$ | $8.9\mbox{$\times 10^{-3}$}$
20 | $1.6\mbox{$\times 10^{-2}$}$ | $1.2\mbox{$\times 10^{-2}$}$ | $8.9\mbox{$\times 10^{-3}$}$ | $7.2\mbox{$\times 10^{-3}$}$ | $6.0\mbox{$\times 10^{-3}$}$
30 | $7.9\mbox{$\times 10^{-3}$}$ | $6.0\mbox{$\times 10^{-3}$}$ | $4.9\mbox{$\times 10^{-3}$}$ | $4.1\mbox{$\times 10^{-3}$}$ | $3.6\mbox{$\times 10^{-3}$}$
Table 3: Relative importance of DSG measured by $S_{\ell n}$ for spherical
harmonic degrees up to $\ell=30$ and $n$ up to $5$.

When specifying a source term $\rho^{U}$, we can quantify the relative
importance of the DSG at each spherical harmonic degree by

$S_{\ell}(N)=\frac{\sum_{n=1}^{N}\,j^{\star}_{\ell n}(R)\;\rho_{\ell
n}^{U}\big{/}\left(k_{\ell
n}^{2}-K^{2}\right)}{\sum_{n=1}^{N}\,j^{\star}_{\ell n}(R)\;\rho_{\ell
n}^{U}\big{/}k_{\ell n}^{2}}\;-\;1\;\;.$ (68)

Figure 6: Measure $S_{\ell}(N)$ quantifying the relative importance of self
gravity at different spherical harmonic degrees $\ell$. Line types indicate
the different radial profiles used for $\rho^{U}$: $f_{1}=r^{5}$ (solid),
$f_{2}$ (dotted), and $f_{3}$ (dashed).

Fig. 6 compares $S_{\ell}$ for the three radial $\rho^{U}$ profiles explored
in Sect. 4. In order to be on the safe side, we have used $N=200$. Selected
values of $S_{\ell}(200)$ are listed in Table 4. All cases show a similar
decay with $\ell$, reaching $10$% relative importance between $\ell=5$ and
$\ell=7$ and $1$% between $\ell=22$ and $\ell=30$. At least for degrees
$\ell>20$, the specific radial profile hardly seems to matter. Because $n=1$
contributions are always significant, the respective ratio (67) listed in
Table 3 already provides a decent estimate of the relative importance for the
DSG.

$\ell$ | $f_{2}=r^{5}$ | $f_{3}$
---|---|---
2 | $7.6\mbox{$\times 10^{-1}$}$ | $6.8\mbox{$\times 10^{-1}$}$
3 | $3.4\mbox{$\times 10^{-1}$}$ | $3.0\mbox{$\times 10^{-1}$}$
4 | $2.0\mbox{$\times 10^{-1}$}$ | $1.9\mbox{$\times 10^{-1}$}$
5 | $1.4\mbox{$\times 10^{-1}$}$ | $1.3\mbox{$\times 10^{-1}$}$
6 | $1.0\mbox{$\times 10^{-1}$}$ | $1.0\mbox{$\times 10^{-1}$}$
8 | $6.4\mbox{$\times 10^{-2}$}$ | $6.4\mbox{$\times 10^{-2}$}$
10 | $4.4\mbox{$\times 10^{-2}$}$ | $4.7\mbox{$\times 10^{-2}$}$
14 | $2.5\mbox{$\times 10^{-2}$}$ | $2.9\mbox{$\times 10^{-2}$}$
20 | $1.3\mbox{$\times 10^{-2}$}$ | $1.8\mbox{$\times 10^{-2}$}$
30 | $6.5\mbox{$\times 10^{-3}$}$ | $1.0\mbox{$\times 10^{-2}$}$
Table 4: Relative importance of the DSG measured by $S_{\ell}$ for two
different radial functions. A radial truncation of $N=200$ has been used.

##  6 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 ${\varPsi}^{\prime}$ is
defined by the inhomogeneous differential equation

$\left(\nabla^{2}+\mu\right)\;{\varPsi}^{\prime}=4\pi G\;\rho^{U}$ (69)

where $\mu$ is the DSG factor and $\rho^{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 $\rho^{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 $\mu=\pi^{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 $\ell=1$. For the background density with
polytropic index unity, the $\ell=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 $\ell=1$ contributions vanish.

Table 5 compares the relative DSG impact with the precision of newest gravity
harmonics of Jupiter and Saturn. The even harmonics $J_{2}$ to $J_{6}$ are not
listed since they are dominated by the rotational deformation of the planet.
For 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.

$\ell$ | Jupiter | Saturn | $S_{\ell}$
---|---|---|---
3 | $0.24$ | $0.39$ | $0.30$
5 | $0.11$ | $0.24$ | $0.13$
7 | $0.14$ | $1.13$ | $0.08$
9 | $0.42$ | $0.70$ | $0.05$
10 | $0.40$ | $0.09$ | $0.05$
11 | $3.39$ | $1.44$ | $0.04$
12 | $3.78$ | $0.67$ | $0.04$
Table 5: Relative error of gravity harmonics for Jupiter (Iess et al., 2018)
(second column) and Saturn (Iess et al., 2019) (third column). The fourth
column shows $S_{\ell}$, the relative impact of the DSG for radial profile
$f_{3}$ also listed in Table 4.

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.

## References

  * Abramowitz and Stegun (1984) Abramowitz, M., Stegun, I., 1984. Pocketbook of mathematical functions. Verlag Harry Deutsch, Thun.
  * Braginsky and Roberts (1995) Braginsky, S. I., Roberts, P. H., 1995. Equations governing convection in earth’s core and the geodynamo. Geophys. Astrophys.l Fluid Dyn. 79, 1–97.
  * Cao and Stevenson (2017) Cao, H., Stevenson, D. J., Apr 2017. Gravity and zonal flows of giant planets: From the Euler equation to the thermal wind equation. J. Geosphys. Res. (Planets) 122 (4), 686–700.
  * Debras and Chabrier (2019) Debras, F., Chabrier, G., feb 2019. New models of jupiter in the context of juno and galileo. APJ 872 (1), 100.
  * French et al. (2012) French, M., Becker, A., Lorenzen, W., Nettelmann, N., Bethkenhagen, M., Wicht, J., Redmer, R., Sep. 2012. Ab Initio Simulations for Material Properties along the Jupiter Adiabat. Astrophys. J. Supp. 202, 5.
  * Galanti et al. (2019) Galanti, E., Kaspi, Y., Miguel, Y., Guillot, T., Durante, D., Racioppa, P., Iess, L., Jan. 2019. Saturn’s Deep Atmospheric Flows Revealed by the Cassini Grand Finale Gravity Measurements. Geosphy. Res. Lett. 46, 616–624.
  * Galanti et al. (2017) Galanti, E., Kaspi, Y., Tziperman, E., Jan 2017. A full, self-consistent treatment of thermal wind balance on oblate fluid planets. J. Fluid Mech. 810, 175–195.
  * Guillot et al. (2018) Guillot, T., Miguel, Y., Militzer, B., Hubbard, W. B., Kaspi, Y., Galanti, E., Cao, H., Helled, R., Wahl, S. M., Iess, L., Folkner, W. M., Stevenson, D. J., Lunine, J. I., Reese, D. R., Biekman, A., Parisi, M., Durante, D., Connerney, J. E. P., Levin, S. M., Bolton, S. J., Mar. 2018. A suppression of differential rotation in Jupiter’s deep interior. Nature 555, 227–230.
  * Hubbard (1982) Hubbard, W. B., Dec 1982. Effects of differential rotation on the gravitational figures of Jupiter and Saturn. Icarus 52 (3), 509–515.
  * Hubbard (2013) Hubbard, W. B., May 2013. Concentric Maclaurin Spheroid Models of Rotating Liquid Planets. APJ 768 (1), 43.
  * Iess et al. (2018) Iess, L., Folkner, W. M., Durante, D., Parisi, M., Kaspi, Y., Galanti, E., Guillot, T., Hubbard, W. B., Stevenson, D. J., Anderson, J. D., Buccino, D. R., Casajus, L. G., Milani, A., Park, R., Racioppa, P., Serra, D., Tortora, P., Zannoni, M., Cao, H., Helled, R., Lunine, J. I., Miguel, Y., Militzer, B., Wahl, S., Connerney, J. E. P., Levin, S. M., Bolton, S. J., Mar. 2018. Measurement of Jupiter’s asymmetric gravity field. Nature 555, 220–222.
  * Iess et al. (2019) Iess, L., Militzer, B., Kaspi, Y., Nicholson, P., Durante, D., Racioppa, P., Anabtawi, A., Galanti, E., Hubbard, W., Mariani, M. J., Tortora, P., Wahl, S., Zannoni, M., Jun 2019. Measurement and implications of Saturn’s gravity field and ring mass. Science 364 (6445), aat2965.
  * Kaspi et al. (2016) Kaspi, Y., Davighi, J. E., Galanti, E., Hubbard, W. B., Sep 2016. The gravitational signature of internal flows in giant planets: Comparing the thermal wind approach with barotropic potential-surface methods. Icarus 276, 170–181.
  * Kaspi et al. (2018) Kaspi, Y., Galanti, E., Hubbard, W. B., Stevenson, D. J., Bolton, S. J., Iess, L., Guillot, T., Bloxham, J., Connerney, J. E. P., Cao, H., Durante, D., Folkner, W. M., Helled, R., Ingersoll, A. P., Levin, S. M., Lunine, J. I., Miguel, Y., Militzer, B., Parisi, M., Wahl, S. M., Mar. 2018. Jupiter’s atmospheric jet streams extend thousands of kilometres deep. Nature 555, 223–226.
  * Kong et al. (2016) Kong, D., Zhang, K., Schubert, G., Oct 2016. Odd gravitational harmonics of Jupiter: Effects of spherical versus nonspherical geometry and mathematical smoothing of the equatorially antisymmetric zonal winds across the equatorial plane. Icarus 277, 416–423.
  * Kong et al. (2017) Kong, D., Zhang, K., Schubert, G., Jul 2017. On the interpretation of the equatorially antisymmetric Jovian gravitational field. MNRAS 469 (1), 716–720.
  * Kong et al. (2018) Kong, D., Zhang, K., Schubert, G., Anderson, J. D., May 2018. Origin of Jupiter’s cloud-level zonal winds remains a puzzle even after Juno. PNAS 115 (34), 8499–8504.
  * Nettelmann (2017) Nettelmann, N., Oct 2017. Low- and high-order gravitational harmonics of rigidly rotating Jupiter. AAP 606, A139.
  * Nettelmann et al. (2012) Nettelmann, N., Becker, A., Holst, B., Redmer, R., May 2012. Jupiter Models with Improved Ab Initio Hydrogen Equation of State (H-REOS.2). Astrophys. J. 750, 52.
  * Wisdom (1996) Wisdom, J., 1996. Non-perturbative hydrostatic equilibrium, available at http://web.mit.edu/wisdom/www/ interior.pdf.
  * Wisdom and Hubbard (2016) Wisdom, J., Hubbard, W. B., Mar 2016. Differential rotation in Jupiter: A comparison of methods. Icarus 267, 315–322.
  * Zhang et al. (2015) Zhang, K., Kong, D., Schubert, G., Jun. 2015. Thermal-gravitational Wind Equation for the Wind-induced Gravitational Signature of Giant Gaseous Planets: Mathematical Derivation, Numerical Method, and Illustrative Solutions. Astrophys. J. 806, 270.
  * Zharkov and Trubitsyn (1978) Zharkov, V. N., Trubitsyn, V. P., 1978. Physics of Planetary Interiors. Pachart, Tucson Ariz.

##  Appendix A Orthogonality

In this section we show that the spherical Bessel functions for different
$k_{\ell n}$ are orthogonal and that $k_{\ell n}^{2}$ is real. We start by
recalling the properties of a self-adjoint or Hemitian linear operator $L$.
Let $f$ and $g$ be eigenvectors (functions) of $L$ with eigenvalues $\lambda$
and $\mu$:

$L\;f=\lambda\;f\;\;,\;\;L\;g=\mu\;g\;\;.$ (70)

For a self-adjoint operator we have

$\langle g,Lf\rangle=\langle Lg,f\rangle\;\;.$ (71)

It follows that

$\lambda\;\langle g,f\rangle=\mu^{\star}\;\langle g,f\rangle$ (72)

and thus $\lambda=\mu^{\star}$. The eigenvalue is thus real and for
$\lambda\neq\mu$ we must have

$\langle f,g\rangle=0\;\;.$ (73)

Here the angular brackets denote the integration over the interval of
interest, in our case

$\langle f,g\rangle=\int_{0}^{R}\,dr\;r^{2}f^{\star}\;g\;\;.$ (74)

To show under which conditions an operator is Hermitian, we chose a somewhat
more general textbook example:

$L=a(r)\frac{\partial^{2}}{\partial r^{2}}\;+\;b(r)\frac{\partial}{\partial
r}+c(r)\;\;.$ (75)

Partial integration yields

$\langle f,L\,g\rangle=\left.\left[r^{2}af\frac{\partial g}{\partial
r}+r^{2}bfrg-g\frac{\partial(r^{2}af)}{\partial r}\right]\right|_{0}^{R}\\\
+\int_{0}^{R}\,dr\,g\,\left[\frac{\partial^{2}(r^{2}af^{\star})}{\partial
r^{2}}-\frac{\partial(r^{2}bf^{\star})}{\partial r}+r^{2}f^{\star}c\right]$
(76)

Rewriting part of the last integral in terms of the operator $L$ leads to

$\langle f,L\,g\rangle=\langle L\,f,g\rangle\\\
+\left.\left[r^{2}af\frac{\partial g}{\partial
r}+r^{2}bfg-r^{2}ag\frac{\partial f}{\partial
r}-fg\frac{\partial(r^{2}a)}{\partial r}\right]\right|_{0}^{r_{o}}\\\
+\int_{0}^{r_{o}}\,dr\ g\,\left[f^{\star}\frac{\partial^{2}(r^{2}a)}{\partial
r^{2}}+2\frac{\partial(r^{2}a)}{\partial r}\frac{\partial f^{\star}}{\partial
r}-f^{\star}\frac{\partial(r^{2}b)}{\partial r}-2r^{2}b\frac{\partial
f^{\star}}{\partial r}\right]$ (77)

The remaining integral vanishes when

$\frac{\partial(r^{2}a)}{\partial r}=r^{2}b\;\;,$ (78)

which is certainly the case for the Laplace operator.

The surface contributions only vanish for particular boundary conditions. When
using equation (78), the surface contributions vanish for:

$f\frac{\partial g}{\partial r}-g\frac{\partial f}{\partial r}=0\;\;.$ (79)

There are the three classical options:

  1. 1.

Dirichlet boundary conditions $f=0$

  2. 2.

Neumann boundary conditions $\partial f/\partial r=0$

  3. 3.

mixed boundary conditions $\partial f/\partial r+df=0$, where $d$ is a
constant.

The third option is used for the gravity problem.

We have thus shown that the different eigenfunctions defined for each
spherical Bessel function $j_{\ell}(k_{\ell n}r)$ (or the second kind
$y_{\ell}(k_{\ell n}r)$) must be orthogonal as long as the functions fulfill
the boundary conditions.

##  Appendix B Normalization

Using

$\langle f,Lg\rangle-\langle Lf,g\rangle=(\mu-\lambda)\,\langle
f,g\rangle\;\;,$ (80)

we can define the integral $\langle f,f\rangle$ as the limit

$\langle f,f\rangle=\lim_{\lambda\rightarrow\mu}\frac{\langle
f,Lg\rangle-\langle Lf,g\rangle}{\mu-\lambda}$ (81)

Using equation (77) shows that

$\langle
f,f\rangle=\lim_{\lambda\rightarrow\mu}\frac{\left.\left[r^{2}af\;\partial
g\big{/}\partial r\;-\;r^{2}ag\;\partial f\big{/}\partial
r\right]\right|_{r_{i}}^{r_{o}}}{\mu-\lambda}$ (82)

This limit can be evaluated using l’Hospital’s rule.

For the spherical Bessel functions and the Laplace operator we are interested
in, equation (82) reads

$\int_{0}^{r_{o}}\,dr\;r^{2}\;j_{\ell}^{2}(kr)=\\\ \lim_{k^{\prime}\rightarrow
k}\frac{r_{o}^{2}\;\left[j_{\ell}(kr_{o})\;\partial
j_{\ell}(k^{\prime}r_{o})\big{/}\partial
r\;-\;j_{\ell}(k^{\prime}r_{o})\;\partial j_{\ell}(kr_{o})\big{/}\partial
r\right]}{{k}^{2}-{k^{\prime}}^{2}}\;\;,$ (83)

where we have used $k=k_{\ell n}$ for brevity.

The result depends on the boundary conditions. For the mixed condition the
limit becomes

$\int_{0}^{r_{o}}\,dr\;r^{2}\;j_{\ell}^{2}(kr)=\\\ \lim_{k^{\prime}\rightarrow
k}\frac{r_{o}^{2}j_{\ell}(kr_{o})\;\left[\partial
j_{\ell}(k^{\prime}r_{o})\big{/}\partial
r\;+\;(\ell+1)\big{/}r_{o}\;j_{\ell}(k^{\prime}r_{o})\right]}{{k}^{2}-{k^{\prime}}^{2}}\;\;.$
(84)

Using recurrence relation (88) and L’Hopital’s rule yields

$\int_{0}^{R}\,dr\;r^{2}\;j_{\ell}^{2}(kr)=-\frac{R^{3}j_{\ell}(kR)\;\partial
j_{\ell-1}(kR)\big{/}\partial r}{2k^{2}}\;\;.$ (85)

Finally, using recurrence relations (89) leads to

$\int_{0}^{R}\,dr\;r^{2}\;j_{\ell}^{2}(kr)=\frac{R^{3}}{2}\;j_{\ell}^{2}(kR)$
(86)

and thus the normalization constant

$N_{\ell n}=\frac{2^{1/2}}{r_{o}^{3/2}\,j_{\ell}(k_{\ell n}r_{o})}\;\;.$ (87)

##  Appendix C Recurrence relations

Some recurrence relations for determining derivatives of spherical Bessel
functions come in handy. Standard relations (Abramowitz and Stegun, 1984, e.
g. ) are

$\partial j_{\ell}(x)\big{/}\partial
x=j_{\ell-1}(x)\;-\;(\ell+1)\big{/}x\;j_{\ell}(x)\;\;,$ (88)

and

$\partial j_{\ell}(x)\big{/}\partial
x=-j_{\ell+1}(x)\;+\;\ell\big{/}x\;j_{\ell}(x)\;\;.$ (89)

Combining both allows us to express the second derivative as

$\partial^{2}j_{\ell}(x)\big{/}\partial
x^{2}=-2\big{/}x\;j_{\ell-1}(x)\;-\left[1-(\ell+1)(\ell+2)\big{/}x^{2}\right]\;j_{\ell}(x)\;\;.$
(90)

##  Appendix D Equivalence of new and classical solution

For $K^{2}=0$, both the classical solution equation (37) and the new expansion
(48)/(50) in spherical Bessel functions should be identical. A comparison
shows that this would require

$\int_{0}^{r_{o}}\,d\,\tilde{r}\;{\tilde{r}}^{\ell+2}\;{{\rho}^{\prime}}(\tilde{r})\;\overset{?}{=}\;\\\
(2\ell+1)\,r_{o}^{\ell+1}\;\sum_{n=1}^{\infty}\;\frac{j^{\star}_{\ell
n}(r)}{k_{\ell
n}^{2}}\;\;\int_{0}^{r_{o}}\,d\,\tilde{r}\;{\tilde{r}}^{2}j^{\star}_{\ell
n}(\tilde{r})\;{{\rho}^{\prime}}(\tilde{r})\;\;,$ (91)

where $j^{\star}_{\ell n}=N_{\ell n}j_{\ell}(k_{\ell 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:

$\tilde{r}^{\ell}=\sum_{n=1}^{\infty}\;j^{\star}_{\ell
n}(\tilde{r})\;\int_{0}^{R}\,d\,r\;r^{\ell+2}\;j^{\star}_{\ell n}(r)\;\;.$
(92)

Partial integration and using the boundary conditions (36) yields

$\int_{0}^{R}\,d\,r\;r^{\ell+2}\;j^{\star}_{\ell
n}(r)=\frac{(2\ell+1)}{k_{\ell
n}}\;\int_{0}^{R}\,d\,r\;r^{\ell+1}\;j^{\star}_{\ell-1n}(r)\;\;.$ (93)

Using recurrence relation (88) and performing another partial integration
finally gives

$\int_{0}^{R}\,d\,r\;r^{\ell+2}\;j^{\star}_{\ell
n}(r)=\frac{(2\ell+1)}{k_{\ell n}^{2}}\;r_{o}^{\ell+1}\;j^{\star}_{\ell
n}(R)\;\;.$ (94)

Plugging this into equation (92) and then the result into the left hand side
of equation (91) finally proves equation (91).