1 Introduction
The observed cloud-level flow on Jupiter and Saturn is dominated by strong east–west (zonal) flows. The depth to which these flows extend is unknown, and has been a topic of great debate over the past few decades (see the reviews by Vasavada & Showman Reference Vasavada and Showman2005 and Showman et al. Reference Showman, Kaspi, Achterberg and Ingersoll2016). One of the prime goals of the Juno mission to Jupiter and the Cassini Grande Finale at Saturn is to estimate the depth of these flows via precise gravitational measurements. If the flows are indeed deep, and therefore perturb significant mass, then they can produce a gravity signal that will be measurable (Hubbard Reference Hubbard1999; Kaspi et al. Reference Kaspi, Hubbard, Showman and Flierl2010). Constraining the depth of these flows will help explore the mechanisms driving the jets (e.g. Busse Reference Busse1976; Williams Reference Williams1978; Ingersoll & Pollard Reference Ingersoll and Pollard1982; Cho & Polvani Reference Cho and Polvani1996; Showman, Gierasch & Lian Reference Showman, Gierasch and Lian2006; Scott & Polvani Reference Scott and Polvani2007; Kaspi & Flierl Reference Kaspi and Flierl2007; Lian & Showman Reference Lian and Showman2010; Liu & Schneider Reference Liu and Schneider2010; Heimpel, Gastine & Wicht Reference Heimpel, Gastine and Wicht2016), and give better constraints on interior structure models (e.g. Guillot Reference Guillot2005; Militzer et al. Reference Militzer, Hubbard, Vorberger, Tamblyn and Bonev2008; Nettelmann et al. Reference Nettelmann, Becker, Holst and Redmer2012; Helled & Guillot Reference Helled and Guillot2013).
Several studies over the past few decades have examined the effects of interior flow on the gravitational moments. The gravity moment spectrum mostly results from the planet’s oblate shape due to its rotation, and from the corresponding interior density distribution. However, density perturbations due to atmospheric dynamics and internal flows can affect the measured gravity moments especially if the flows extend deep enough into the planets. Hubbard (Reference Hubbard1982) and Hubbard (Reference Hubbard1999) used potential theory to calculate the gravity moments due to internal flows, by extending the observed cloud-level winds along cylinders throughout the planet as suggested by Busse (Reference Busse1976). This approach takes into account the oblateness of the planet, yet is only possible for the barotropic case, meaning that the flow is constant along lines parallel to the axis of rotation. This occurs if the baroclinicity vector
$\unicode[STIX]{x1D735}\unicode[STIX]{x1D70C}\times \unicode[STIX]{x1D735}p$
vanishes (e.g. if density is a function of pressure only) at small Rossby number and negligible dissipation. More recently, Hubbard introduced more accurate calculations, for the gravitational signature of the flows (Hubbard Reference Hubbard2012; Kong, Zhang & Schubert Reference Kong, Zhang and Schubert2012; Hubbard Reference Hubbard2013), using concentric Maclaurin spheroids (CMS), but these are also only limited to the fully barotropic case.
A different approach, assuming the large scale flow is dominated by the rotation of the planet, used thermal wind (TW) balance to calculate the gravity moments due to the wind field (Kaspi et al.
Reference Kaspi, Hubbard, Showman and Flierl2010; Kaspi Reference Kaspi2013; Kaspi et al.
Reference Kaspi, Showman, Hubbard, Aharonson and Helled2013; Liu, Schneider & Kaspi Reference Liu, Schneider and Kaspi2013; Liu, Schneider & Fletcher Reference Liu, Schneider and Fletcher2014). The TW approach is not limited to the barotropic case (it can account for any wind field), and in the barotropic limit has been shown to be equivalent to the potential theory and CMS methods (Kaspi et al.
Reference Kaspi, Davighi, Galanti and Hubbard2016). In addition, this approach allows for the calculation of the odd gravitational moments, which can emerge from north–south hemispherical asymmetries in the wind structure (Kaspi Reference Kaspi2013). However, this TW approach was originally limited to spherical symmetry, resulting in an inability to calculate the effects of the planet oblateness on the gravity signature of the winds. Cao & Stevenson (Reference Cao and Stevenson2016) added the effects of oblateness on the background state density and gravity and concluded that it should be considered when estimating the effects of the winds on the gravity moments using TW. Similarly, Zhang, Kong & Schubert (Reference Zhang, Kong and Schubert2015) included another effect, of the gravity anomalies due to density perturbations associated with the winds, and found an effect on the second gravity moment
$J_{2}$
, terming their approach the thermal-gravity wind (TGW) method. However, while both of these recent studies found some effects of the terms they added, their choice of added physics did not result from a systematic and self-consistent approach.
The purpose of this study is to develop a full, self-consistent thermal wind (FTW) perturbation approach for the treatment of the general TW balance on a fluid planet. This will allow to calculate density anomalies and gravity moments due to prescribed winds, omitting the traditional sphericity assumptions which have been adopted from dynamics on terrestrial planets. Our approach includes the effect of oblateness, as well as that of gravity anomalies due to the dynamical density perturbations themselves. We show these two effects to be but two of several different factors that should be considered in a self-consistent calculation. Our approach is based on a systematic perturbation expansion, which both allows us to consider all effects, and also points the way to improving the estimated gravity moments using a higher-order perturbation that can be considered by future studies.
By including all relevant terms in the general TW balance, we are able to evaluate the relative contribution of different terms. We find that the simplified TW approach captures most of the relation between the wind shear and density gradients. The term added in TGW is found to be one of several smaller terms that all need to be added together for consistency in order to improve the estimates of the simplified TW balance. Furthermore, previous applications of the TW balance encountered an unknown integration constant that was a function of radius only and could not be solved for. We show that this integration constant may have an effect, although small, on the gravity moments, and develop a method for calculating it.
The following section describes the perturbation approach, the resulting equations and how they are solved. Next, in § 3, we first verify this approach by comparing it with the results of the CMS method in the barotropic limit, and then compare our results to the less complete approaches of simplified TW and TGW. We then also apply the self-consistent solution to a case with baroclinic winds where CMS cannot be used. We discuss the results and conclude in § 4.
2 Methods: perturbation expansion of the momentum equations
We begin by taking the standard form of the momentum equations on a planet rotating at an angular velocity
$\unicode[STIX]{x1D6FA}$
,
$$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{u}}{\unicode[STIX]{x2202}t}+(\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{u}+2\unicode[STIX]{x1D734}\times \boldsymbol{u}+\unicode[STIX]{x1D734}\times \unicode[STIX]{x1D734}\times \boldsymbol{r}=-\frac{1}{\unicode[STIX]{x1D70C}}\unicode[STIX]{x1D735}p+\unicode[STIX]{x1D735}\unicode[STIX]{x1D6F7}, & & \displaystyle\end{eqnarray}$$
where
$\boldsymbol{u}$
is the three-dimensional wind vector,
$\unicode[STIX]{x1D70C}$
is density,
$p$
is pressure,
$\unicode[STIX]{x1D734}$
is the planetary rotation rate and
$\unicode[STIX]{x1D6F7}$
is the body force potential. The first term on the left-hand side is the local acceleration of the flow, the second is the Eulerian advection, the third is the Coriolis acceleration and the fourth is the centrifugal acceleration. On the right-hand side appear the pressure gradient term and the body force (gravity in this case, so that
$\unicode[STIX]{x1D735}\unicode[STIX]{x1D6F7}=-\boldsymbol{g}$
). Note that by gravity we refer here to the force due to the Newtonian potential, not to the modified gravity which is commonly used in geostrophic studies and includes the centrifugal potential as well. Typical values for a Jupiter-like planet are
$U=O(100)~\text{m}~\text{s}^{-1}$
,
$\unicode[STIX]{x1D6FA}=O(10^{-4})~\text{s}^{-1}$
,
$a=O(7\times 10^{7})$
m, where
$a$
is the planet radius. The resulting Rossby number (
$Ro$
) is therefore much smaller than one (
$Ro\equiv U/\unicode[STIX]{x1D6FA}a\approx 10^{-2}$
), and the first two terms in (2.1) can be neglected so that the resulting balance is
$$\begin{eqnarray}\displaystyle 2\unicode[STIX]{x1D734}\times (\unicode[STIX]{x1D70C}\boldsymbol{u})=-\unicode[STIX]{x1D735}p-\unicode[STIX]{x1D70C}\boldsymbol{g}-\unicode[STIX]{x1D70C}\unicode[STIX]{x1D734}\times \unicode[STIX]{x1D734}\times \boldsymbol{r}. & & \displaystyle\end{eqnarray}$$
Next, we denote the static solution (
$\boldsymbol{u}=0$
) as
$\unicode[STIX]{x1D70C}_{0},p_{0},g_{0}$
, and the perturbations due to the non-zero wind (dynamical part of the solution) as
$\unicode[STIX]{x1D70C}^{\prime },p^{\prime },g^{\prime }$
, such that
$$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\unicode[STIX]{x1D70C}=\unicode[STIX]{x1D70C}_{0}(r,\unicode[STIX]{x1D703})+\unicode[STIX]{x1D70C}^{\prime }(r,\unicode[STIX]{x1D703}),\\ p=p_{0}(r,\unicode[STIX]{x1D703})+p^{\prime }(r,\unicode[STIX]{x1D703}),\\ \boldsymbol{g}=\boldsymbol{g}_{0}(r,\unicode[STIX]{x1D703})+\boldsymbol{g}^{\prime }(r,\unicode[STIX]{x1D703}).\end{array}\right\} & & \displaystyle\end{eqnarray}$$
Note that both the static and the dynamic solutions are functions of latitude and radius, and that the gravity is directly related to the density via a relation shown below.
The equation obtained by setting the small Rossby number to zero as a first approximation is effectively static and does not include the velocity field,
$$\begin{eqnarray}\displaystyle 0=-\unicode[STIX]{x1D735}p_{0}-\unicode[STIX]{x1D70C}_{0}\boldsymbol{g}_{0}-\unicode[STIX]{x1D70C}_{0}\unicode[STIX]{x1D734}\times \unicode[STIX]{x1D734}\times \boldsymbol{r}, & & \displaystyle\end{eqnarray}$$
and the dynamical perturbations therefore satisfy
$$\begin{eqnarray}\displaystyle 2\unicode[STIX]{x1D734}\times (\unicode[STIX]{x1D746}_{0}\boldsymbol{u}) & = & \displaystyle -\unicode[STIX]{x1D735}p^{\prime }-\unicode[STIX]{x1D70C}_{0}\boldsymbol{g}^{\prime }-\unicode[STIX]{x1D70C}^{\prime }\boldsymbol{g}_{0}-\unicode[STIX]{x1D70C}^{\prime }\unicode[STIX]{x1D734}\times \unicode[STIX]{x1D734}\times \boldsymbol{r}.\end{eqnarray}$$
The solution procedure outlined here involves first finding the static solution and then solving (2.5) for the dynamical perturbations to the density due to the effects of the prescribed winds. Taking the curl of (2.5) yields a single equation in the azimuthal direction
$$\begin{eqnarray}\displaystyle -2\unicode[STIX]{x1D6FA}r\unicode[STIX]{x2202}_{z}(\unicode[STIX]{x1D70C}_{0}u) & = & \displaystyle -rg_{0}^{(\unicode[STIX]{x1D703})}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}^{\prime }}{\unicode[STIX]{x2202}r}-g_{0}^{(r)}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}^{\prime }}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}+r\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{0}}{\unicode[STIX]{x2202}r}{g^{\prime }}^{(\unicode[STIX]{x1D703})}\nonumber\\ \displaystyle & & \displaystyle -\,g^{\prime (r)}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}\unicode[STIX]{x1D6FA}^{2}r\left[\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}^{\prime }}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}\cos ^{2}\unicode[STIX]{x1D703}+\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}^{\prime }}{\unicode[STIX]{x2202}r}r\cos \unicode[STIX]{x1D703}\sin \unicode[STIX]{x1D703}\right],\end{eqnarray}$$
where the notation
$\unicode[STIX]{x2202}_{z}=\cos \unicode[STIX]{x1D703}(\unicode[STIX]{x2202}/\unicode[STIX]{x2202}r)-\sin \unicode[STIX]{x1D703}(\unicode[STIX]{x2202}/\unicode[STIX]{x2202}\unicode[STIX]{x1D703})$
denotes the derivative along the direction of the axis of rotation (
$z$
), and gravity is expressed as function of the density as
$$\begin{eqnarray}\boldsymbol{g}(r,\unicode[STIX]{x1D703})=2\unicode[STIX]{x03C0}G\left[\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r},\frac{\unicode[STIX]{x2202}}{r\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}\right]\iint \left\langle \frac{1}{|\boldsymbol{r}-\boldsymbol{r}^{\prime }|}\right\rangle \unicode[STIX]{x1D70C}(r^{\prime },\unicode[STIX]{x1D703}^{\prime })r^{\prime 2}\cos \unicode[STIX]{x1D703}^{\prime }\,\text{d}\unicode[STIX]{x1D703}^{\prime }\,\text{d}r^{\prime },\end{eqnarray}$$
where the gravity can be either
$\boldsymbol{g}_{\mathbf{0}}$
or
$\boldsymbol{g}^{\prime }$
, calculated from
$\unicode[STIX]{x1D70C}_{0}$
or
$\unicode[STIX]{x1D70C}^{\prime }$
, respectively, and
$G=6.67\times 10^{-11}~\text{m}^{3}~\text{kg}^{-1}~\text{s}^{-2}$
is the gravitational constant. Note that (2.6), together with (2.7), form an integro-differential equation whose solution requires the calculation of integration constants as is done below and demonstrated in a simpler context in appendix B.
The static model solution, as function of radius and latitude, from the CMS model for the density and gravity (a–c), and their deviation from the latitudinal mean (d–f).
The above equations for the perturbation density are the first-order perturbation equations, which are solved in this study. It is important to note, though, that because this is a self-consistent treatment of the density perturbations, it also allows improving on the approximation by proceeding to the next orders. As a demonstration, we write the second-order perturbation equations in appendix A.
2.1 The background static solution
The static density
$\unicode[STIX]{x1D70C}_{0}$
and gravity
$\boldsymbol{g}_{\mathbf{0}}$
are taken from the solution of the CMS model (Hubbard Reference Hubbard2012, Reference Hubbard2013). The model is based on a numerical method for solving the equilibrium shape of a rotating planet, for which an analytic solution exists in the form of a Maclaurin spheroid. The continuously varying density and pressure structures are represented by a discrete set of layers, in which the density and pressure are constant. This onion-like structure is then decomposed into a set of CMS, and a solution is sought by requiring that the sum of the gravitational potential and the rotational potential are constant on the surface of the planet (Hubbard Reference Hubbard2012, Reference Hubbard2013). Solutions using this method give similar results to other methods (Wisdom & Hubbard Reference Wisdom and Hubbard2016; Kaspi et al.
Reference Kaspi, Davighi, Galanti and Hubbard2016). The static fields resulting from the CMS solution are shown in figure 1. Both the density and gravity fields show a structure that is mostly radial. While the density (figure 1
a) ranges from zero to about
$4000~\text{kg}~\text{m}^{-3}$
, its latitude-dependent component (figure 1
d) ranges only between
$-$
150 and
$100~\text{kg}~\text{m}^{-3}$
. Similarly, the radial gravity is also dominated by its radial component (figure 1
b,e), with the peak gravity being
${\sim}32~\text{m}~\text{s}^{-2}$
at about
$0.7a$
. The latitudinal component of the gravity is much weaker than the radial component (figure 1
c,f). Comparing figure 1(c) and (f) shows that the latitudinal mean of the latitudinal component of the gravity is zero, and has a much smaller magnitude than the radial component. Note that in figure 1(c,f) a positive value means gravity pointing northward.
2.2 Solving for the dynamic density perturbations
We solve (2.6) by writing the equations in matrix form (e.g. Zhang et al.
Reference Zhang, Kong and Schubert2015). The two-dimensional problem (radius and latitude) is discretized in both directions, where
$N_{r}$
and
$N_{\unicode[STIX]{x1D703}}$
are the number of grid points in radius and latitude, respectively. Here
$N=N_{r}\times N_{\unicode[STIX]{x1D703}}$
is the total number of grid points. Equation (2.6) is then written in matrix form,
$$\begin{eqnarray}\displaystyle \boldsymbol{b}=\unicode[STIX]{x1D63C}\unicode[STIX]{x1D746}^{\prime }, & & \displaystyle\end{eqnarray}$$
where
$\unicode[STIX]{x1D63C}$
is a
$N\times N$
matrix with contributions from all terms in the right-hand side of the equation,
$\boldsymbol{b}$
is the known left-hand side of the equation, due to the prescribed wind field, written as an
$N\times 1$
vector, and the unknown density perturbation
$\unicode[STIX]{x1D746}^{\prime }$
is written as an
$N\times 1$
vector. All partial derivatives are written in centre finite difference form, aside from near the boundaries where the derivatives are evaluated between grid points and weighted together with the adjacent derivative. Note that the gravity is fully calculated in the matrix (2.8). It is done explicitly for each grid point.
Solving (2.8) involves the inversion of the matrix
$\unicode[STIX]{x1D63C}$
, which is not possible because the matrix is singular and thus has a null space, which leads to a part of the solution that cannot be determined from the equation. The physical source for the singularity (and, hence, the existence of the null space) can be better understood in the simpler case where (2.6) is reduced to the TW balance with a spherical base state, involving only the left-hand side and the second term on the right-hand side. In that case, the equation can be integrated in latitude, leaving an unknown integration constant that is a function of radius alone (e.g. Kaspi et al.
Reference Kaspi, Hubbard, Showman and Flierl2010). Had we solved that equation numerically instead of integrating it, we would thus find a null space, corresponding to the unknown integration constant, whose dimension is the number of the radial grid points
$N_{r}$
. In the more general case solved here, while the base state is a function of both radius and latitude, the null space still has a dimension equal to the number of radial grid points
$N_{r}$
, similar to the unknown function of radius only, but which has latitudinal dependence as well. We show in appendix B how the integration constant is calculated analytically in the simpler case, and discuss in § 2.3 the solution in the more general case. In order to solve for the non-null part of the solution, we use singular value decomposition (Strang Reference Strang2006) as follows. Let
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D63C}=\unicode[STIX]{x1D650}\unicode[STIX]{x1D72E}\unicode[STIX]{x1D651}^{\text{T}}, & & \displaystyle\end{eqnarray}$$
where
$\unicode[STIX]{x1D650}$
and
$\unicode[STIX]{x1D651}$
are unitary matrixes, and
$\unicode[STIX]{x1D72E}$
is a rectangular diagonal matrix (whose dimensions are those of
$\unicode[STIX]{x1D63C}$
and whose entries are all zeros except along the diagonal), then the pseudo-inverse of
$\unicode[STIX]{x1D63C}$
is given by
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D63C}^{\dagger }=\unicode[STIX]{x1D651}\unicode[STIX]{x1D72E}^{\dagger }\unicode[STIX]{x1D650}^{\text{T}}, & & \displaystyle\end{eqnarray}$$
where
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D72E}=\text{diag}(\unicode[STIX]{x1D70E}_{1},\ldots ,\unicode[STIX]{x1D70E}_{N-N_{r}},0,\ldots ,0), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D72E}^{\dagger }=\text{diag}(\unicode[STIX]{x1D70E}_{1}^{-1},\ldots ,\unicode[STIX]{x1D70E}_{N-N_{r}}^{-1},0,\ldots ,0). & \displaystyle\end{eqnarray}$$
$\unicode[STIX]{x1D70E}_{i}$
, are the square root of the eigenvalues of
$AA^{\text{T}}$
or
$A^{\text{T}}A$
. The solution is now written as
$\unicode[STIX]{x1D70C}^{\prime }=\hat{\unicode[STIX]{x1D70C}}^{\prime }+\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}^{\prime }$
, such that
$\hat{\unicode[STIX]{x1D70C}}^{\prime }$
is the part obtained from the pseudo-inverse, 
$$\begin{eqnarray}\displaystyle \hat{\unicode[STIX]{x1D746}}^{\prime }=\unicode[STIX]{x1D63C}^{\dagger }\boldsymbol{b}, & & \displaystyle\end{eqnarray}$$
and the additional component
$\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}^{\prime }$
is in the null space of
$A$
, so that
$\unicode[STIX]{x1D63C}\unicode[STIX]{x1D746}^{\prime }=\unicode[STIX]{x1D63C}(\hat{\unicode[STIX]{x1D746}}^{\prime }+\unicode[STIX]{x1D739}\unicode[STIX]{x1D746}^{\prime })=\unicode[STIX]{x1D63C}\hat{\unicode[STIX]{x1D746}}^{\prime }=\boldsymbol{b}$
. That is, denoting the eigenvectors of
$\unicode[STIX]{x1D63C}$
by
$\boldsymbol{e}_{i}$
, such that
$\unicode[STIX]{x1D63C}\boldsymbol{e}_{i}=\unicode[STIX]{x1D706}_{i}\boldsymbol{e}_{i}$
, then the null space of the solution corresponds to any linear combination of eigenvectors corresponding to the zero eigenvalues, which may be added to the solution
$\hat{\unicode[STIX]{x1D746}^{\prime }}$
while still satisfying
$\unicode[STIX]{x1D63C}\unicode[STIX]{x1D746}^{\prime }=\boldsymbol{b}$
. We next discuss the calculation of the null space and its contribution to the solution for
$\unicode[STIX]{x1D746}^{\prime }$
.
2.3 Calculating the null-space solution
As mentioned above, the number of null eigenvalues is always found to be
$N_{n}=N_{r}$
, i.e. equal to the number of grid points in the radial direction, hinting to the possibility that the modes are predominantly radially dependent. Indeed, looking at the zero eigenvectors does show this characteristic (see the example in figure 2
a). Nevertheless, since the null eigenvectors do depend on latitude (figure 2
b), and this dependence is concentrated close to the planet upper levels, it could have a substantial effect on the gravity moments. The contribution of the null space to the density may be written as
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D739}\unicode[STIX]{x1D746}^{\prime }=\mathop{\sum }_{i=N-N_{r}+1}^{N}c_{i}\boldsymbol{e}_{i}=\unicode[STIX]{x1D640}\boldsymbol{c}, & & \displaystyle\end{eqnarray}$$
where
$\unicode[STIX]{x1D640}$
is a
$N\times N_{n}$
matrix whose columns are the
$N_{n}$
null eigenvectors, and
$\boldsymbol{c}$
is a
$N_{n}\times 1$
vector of the unknown amplitudes of the null eigenvectors. The null space discussed above is now represented by the vector
$\boldsymbol{c}$
, and as explained above it cannot be determined using the TW balance alone, and we must therefore introduce additional physics. This additional physics is the relation between the pressure and the density, as well as a constraint on the total mass. We specifically use the polytropic relation (see the further discussion below). A demonstration of a full analytical solution using the polytropic equation is shown in appendix B for a simpler case, and here we present the numerical procedure for solving the general problem for
$\boldsymbol{c}$
. Note that the polytropic relation is only used for the solution of the null space, but the rest of the solution is independent of it.
An example of a null mode
$\boldsymbol{V}$
: (a) the full null model; (b) the latitudinal dependent part of the null mode (
$\boldsymbol{V}_{a}=\boldsymbol{V}-\overline{\boldsymbol{V}}$
). Note that in most of the domain, the structure of the mode depends on radius only, aside from the region close to the surface where a large-scale latitude-dependent structure exists, but whose value is around half the magnitude of the full structure.
The polytropic equation with index unity,
$p=K\unicode[STIX]{x1D70C}^{2}$
, linearized around the static solution, gives
$$\begin{eqnarray}p^{\prime }=K2\unicode[STIX]{x1D70C}_{0}\unicode[STIX]{x1D70C}^{\prime },\end{eqnarray}$$
which may be used to define
$\hat{p}^{\prime }$
as
$$\begin{eqnarray}\hat{p}^{\prime }\equiv Kn2\unicode[STIX]{x1D70C}\hat{\unicode[STIX]{x1D70C}}^{\prime }.\end{eqnarray}$$
Previously we took the curl of the momentum equation (2.5), to solve for
$\hat{\unicode[STIX]{x1D70C}}^{\prime }(r,\unicode[STIX]{x1D703})$
, and it satisfies the original equation up to a gradient of some scalar function
$\unicode[STIX]{x1D713}$
. We therefore may replace the pressure
$\hat{p}^{\prime }$
with
$\hat{p}^{\prime }+\unicode[STIX]{x1D713}$
. Using this augmented function, the perturbation momentum equation may be written as
$$\begin{eqnarray}2\unicode[STIX]{x1D734}\times \boldsymbol{u}\unicode[STIX]{x1D70C}_{0}=-\unicode[STIX]{x1D735}(\hat{p}^{\prime }+\unicode[STIX]{x1D713})-\unicode[STIX]{x1D70C}_{0}\hat{\boldsymbol{g}}^{\prime }-\hat{\unicode[STIX]{x1D70C}}^{\prime }\boldsymbol{g}_{0}-\hat{\unicode[STIX]{x1D70C}}^{\prime }\unicode[STIX]{x1D734}\times \unicode[STIX]{x1D734}\times \boldsymbol{r}.\end{eqnarray}$$
Next, taking the difference between the momentum equation for
$\unicode[STIX]{x1D70C}^{\prime }$
and for
$\hat{\unicode[STIX]{x1D70C}}^{\prime }$
, we find
$$\begin{eqnarray}0=-\unicode[STIX]{x1D735}(p^{\prime }-\hat{p}^{\prime }-\unicode[STIX]{x1D713})-\unicode[STIX]{x1D70C}_{0}\unicode[STIX]{x1D739}\boldsymbol{g}^{\prime }-\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}^{\prime }\boldsymbol{g}_{0}-\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}^{\prime }\unicode[STIX]{x1D734}\times \unicode[STIX]{x1D734}\times \boldsymbol{r}.\end{eqnarray}$$
The function
$\unicode[STIX]{x1D713}$
appears as a correction to the perturbation pressure, and we therefore define the perturbation pressure to be
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FF}p^{\prime }\equiv p^{\prime }-\hat{p}^{\prime }-\unicode[STIX]{x1D713}. & & \displaystyle\end{eqnarray}$$
Taking the difference between the polytropic equation for
$\unicode[STIX]{x1D70C}$
and
$\hat{\unicode[STIX]{x1D70C}}^{\prime }$
, we have
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FF}p^{\prime }=K2\unicode[STIX]{x1D70C}_{0}\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}^{\prime }, & & \displaystyle\end{eqnarray}$$
which leads to an equation for the unknown perturbation density,
$$\begin{eqnarray}0=-\unicode[STIX]{x1D735}(Kn\unicode[STIX]{x1D70C}_{0}\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}^{\prime })-\unicode[STIX]{x1D70C}_{0}\unicode[STIX]{x1D739}\boldsymbol{g}^{\prime }-\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}^{\prime }\boldsymbol{g}_{0}-\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}^{\prime }\unicode[STIX]{x1D734}\times \unicode[STIX]{x1D734}\times \boldsymbol{r}\end{eqnarray}$$
and, explicitly,
$$\begin{eqnarray}\displaystyle & \displaystyle F_{r}:\quad 0=-2K\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}(\unicode[STIX]{x1D70C}_{0}\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}^{\prime })-\unicode[STIX]{x1D70C}_{0}\unicode[STIX]{x1D6FF}g^{\prime (r)}-\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}^{\prime }g_{0}^{(r)}+\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}^{\prime }\unicode[STIX]{x1D6FA}^{2}r\cos ^{2}\unicode[STIX]{x1D703}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle F_{\unicode[STIX]{x1D703}}:\quad 0=-2K\frac{\unicode[STIX]{x2202}}{r\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}(\unicode[STIX]{x1D70C}_{0}\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}^{\prime })-\unicode[STIX]{x1D70C}_{0}\unicode[STIX]{x1D6FF}g^{\prime (\unicode[STIX]{x1D703})}-\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}^{\prime }g_{0}^{(\unicode[STIX]{x1D703})}-\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}^{\prime }\unicode[STIX]{x1D6FA}^{2}r\cos \unicode[STIX]{x1D703}\sin \unicode[STIX]{x1D703}, & \displaystyle\end{eqnarray}$$
$F_{r}$
and
$F_{\unicode[STIX]{x1D703}}$
are the equations in the radial and latitudinal directions, respectively. This represents the radial and latitudinal components of a vector equation, and we next take the divergence to obtain a single equation, 
$$\begin{eqnarray}\displaystyle 0=\frac{1}{r^{2}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}(r^{2}F_{r})+\frac{1}{r\sin \unicode[STIX]{x1D703}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}(\sin \unicode[STIX]{x1D703}F_{\unicode[STIX]{x1D703}}), & & \displaystyle\end{eqnarray}$$
where everything but
$\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}^{\prime }$
is known. Numerically, this equation can be written as a set of linear equations,
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D63D}\unicode[STIX]{x1D739}\unicode[STIX]{x1D746}^{\prime }=0, & & \displaystyle\end{eqnarray}$$
where
$\unicode[STIX]{x1D63D}$
is an
$N\times N$
matrix, and the right-hand side is an
$N\times 1$
vector with zeros in all entries. An additional constraint is that the total mass of the planet must not change due to the existence of the wind, so that
$\int \unicode[STIX]{x1D70C}^{\prime }d^{3}\boldsymbol{r}=0$
and, therefore,
$$\begin{eqnarray}\int \unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}^{\prime }\,\text{d}^{3}\boldsymbol{r}=-\int \hat{\unicode[STIX]{x1D70C}}^{\prime }\,\text{d}^{3}\boldsymbol{r}=-\unicode[STIX]{x1D6FF}M.\end{eqnarray}$$
The mass
$\unicode[STIX]{x1D6FF}M$
, due to the non-null space solution
$\hat{\unicode[STIX]{x1D70C}}^{\prime }$
, is known from the solution to (2.6). Adding this constraint to (2.24) results in augmenting the matrix
$\unicode[STIX]{x1D63D}$
with an additional row, to form a matrix
$\tilde{\unicode[STIX]{x1D63D}}$
whose size is now
$(N+1)\times N$
and the right-hand side, now denoted by
$\boldsymbol{m}$
, is now a vector of length
$N+1$
with a non-zero value only in the last entry,
$m_{N+1}=-\unicode[STIX]{x1D6FF}M$
. Using the definition for the null space part of the solution (2.14) we obtain
$$\begin{eqnarray}\displaystyle \tilde{\unicode[STIX]{x1D63D}}\unicode[STIX]{x1D640}\boldsymbol{c}=\boldsymbol{m}. & & \displaystyle\end{eqnarray}$$
This is a formally overdetermined problem, which is solved for
$\boldsymbol{c}$
using least squares (Strang Reference Strang2006), allowing us to then calculate
$\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}^{\prime }$
. Note that even though
$\tilde{\unicode[STIX]{x1D63D}},\unicode[STIX]{x1D640}$
and
$\boldsymbol{m}$
are all complex,
$\unicode[STIX]{x1D739}\unicode[STIX]{x1D746}^{\prime }$
is found to be real, as expected.
The use of a polytropic relation between the pressure and density implies that the baroclinic vector
$\unicode[STIX]{x1D735}p\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D70C}$
vanishes, and therefore that the velocity field is necessarily barotropic at small Rossby numbers. This is in line with most of the cases discussed here that are indeed barotropic, aside from the last case that is baroclinic and is analysed in § 3.2.1. Note that a different pressure–density relation is also possible. Nonetheless, in all cases discussed here the contribution of the null space solution to the overall solution is negligible. Our purpose here is to show how adding information regarding the equation of state (in this case
$p(\unicode[STIX]{x1D70C})$
) can be used to calculate the unknown integration constant arising in the TW formulation (e.g. Kaspi et al.
Reference Kaspi, Hubbard, Showman and Flierl2010; Zhang et al.
Reference Zhang, Kong and Schubert2015). However, in future application, one would need to use a more realistic equation of state that allows for determining the null space for a baroclinic wind field as well.
2.4 Prescribed winds
The wind profile used in this study is based on the measured cloud-tracking wind during the Cassini flyby (Porco et al. Reference Porco, West, McEwen, Del Genio, Ingersoll, Thomas, Squyres, Dones, Murray and Johnson2003). Since we compare our results to the CMS model solution as a reference for the full oblate solution, and the CMS wind profile must be truncated for numerical convergence (see Kaspi et al. Reference Kaspi, Davighi, Galanti and Hubbard2016 for details), we use a 24th-degree expansion of its differential potential. Kaspi et al. (Reference Kaspi, Davighi, Galanti and Hubbard2016) shows a comparison between the resulting gravity moments using the truncated and untruncated wind profiles. The choice of the specific wind profile does not affect the results. In order to do a proper comparison with the CMS model, which is limited to barotropic winds, the wind profile is extended along cylinders parallel to the direction of the axis of rotation. For the baroclinic case discussed in (§ 3.3), the wind profile is extended toward the centre of the planet using an exponential decay function (e.g. as in Galanti & Kaspi Reference Galanti and Kaspi2016a ) with a decay scale height of 1000 km.
2.5 Calculating gravity moments
In all cases discussed in the following, in addition to examining the solution for the density perturbations, we calculate the resulting gravitational moments given by
$$\begin{eqnarray}\unicode[STIX]{x0394}J_{n}=-\frac{2\unicode[STIX]{x03C0}}{Ma^{n}}\int _{0}^{a}r^{\prime n+2}\,\text{d}r^{\prime }\int _{-1}^{1}P_{n}(\unicode[STIX]{x1D707}^{\prime })\unicode[STIX]{x1D70C}^{\prime }(r^{\prime },\unicode[STIX]{x1D707}^{\prime })\,\text{d}\unicode[STIX]{x1D707}^{\prime },\end{eqnarray}$$
where
$M$
is the mass of the planet,
$P_{n}$
are the Legendre polynomials and
$\unicode[STIX]{x1D707}=\cos (\unicode[STIX]{x1D703})$
. Note that any part of
$\unicode[STIX]{x1D70C}^{\prime }$
that is a function of radius only does not contribute to the gravity moments. For instance, using the latitudinal average of the density,
$\unicode[STIX]{x1D70C}_{m}^{\prime }(r)=\overline{\unicode[STIX]{x1D70C}^{\prime }(r,\unicode[STIX]{x1D703})}$
in (2.28) will give
$$\begin{eqnarray}\unicode[STIX]{x0394}J_{n}=-\frac{2\unicode[STIX]{x03C0}}{Ma^{n}}\int _{0}^{a}r^{\prime n+2}\unicode[STIX]{x1D70C}_{m}^{\prime }(r^{\prime })\,\text{d}r^{\prime }\int _{-1}^{1}P_{n}(\unicode[STIX]{x1D707}^{\prime })\,\text{d}\unicode[STIX]{x1D707}^{\prime }=0,\end{eqnarray}$$
which vanishes due to the Legendre polynomials having a zero latitudinal mean for any value of
$n$
. Therefore, any solution for
$\unicode[STIX]{x1D70C}^{\prime }$
needs to be examined with respect to its latitudinal-dependent part.
3 Results for wind-induced density and gravity moments
We now consider the solution for the density field and gravitational moments in several cases. First, we examine the case of barotropic winds where the results of the perturbation approach can be compared with the CMS solution (§ 3.1). Second, we compare our approach to earlier methods and approximations (§ 3.2). Finally, we analyse an example of the more general case of baroclinic winds, where a CMS solution is not possible (§ 3.3).
3.1 Verification of perturbation method via a comparison with CMS
Solving numerically (2.6), with all six terms on the right-hand side included and adding the null-space solution, we obtain the anomalous density field
$\unicode[STIX]{x1D70C}^{\prime }=\hat{\unicode[STIX]{x1D70C}}+\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}^{\prime }$
from which we calculate the gravitational moments shown in figure 3 (blue line), together with the reference CMS solution (red line). The perturbation analysis captures most of the signal of the moments. The dashed line shows the contribution of the null-mode solution
$\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}^{\prime }$
that is much smaller than the total solution.
The wind-induced gravitational moments solution of (2.28) when all terms are included (blue), compared with the CMS solution (red). Also shown are the contribution from
$\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}^{\prime }$
(black dashed line) and the solid-body-induced gravitational moments (black).
Next, consider each term in (2.6) as function of radius and latitude (figure 4), where figure 4(a) shows the left-hand side, figure 4(b) the total of the right-hand side and figure 4(c–h) show the individual contribution from the six different terms on the right-hand side. The dominant term on the right-hand side balancing the left-hand side is the second term (figure 4
d), i.e. the TW term, whose magnitude is about 10 times larger than any of the other terms. In § 3.2, we examine less-complete solutions each including only some of the terms in (2.6), including the simplified TW approach of Kaspi et al. (Reference Kaspi, Hubbard, Showman and Flierl2010) and the TGW approximation of Zhang et al. (Reference Zhang, Kong and Schubert2015). It is already clear, though, that the TGW term (figure 4
e) is of the same magnitude as several others, so a self-consistent approximation can either neglect all terms, but the most dominant one as done in Kaspi et al. (Reference Kaspi, Hubbard, Showman and Flierl2010), or include all other terms as well as done here in the FTW approach. The perturbation density solution,
$\unicode[STIX]{x1D70C}^{\prime }$
, is concentrated near the surface to a large degree (figure 5
a,b), and therefore terms that depend on its vertical derivative are also concentrated near the surface. Terms that depend on gradients of the zeroth-order density are characterized by a larger-scale structure (figure 4
e,f).
Solution of (2.6) when all terms are kept: (a) left-hand side term; (b) total right-hand side; (c–h) the six terms on the right-hand side. Note the different scales in the different panels.
Figure 5(c,d) shows the solution to the null-space part of the density,
$\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}^{\prime }$
. It is negative everywhere, in order to compensate for
$\hat{\unicode[STIX]{x1D70C}}^{\prime }$
which is generally positive so that mass is conserved (§ 2.3). The null-space solution
$\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}^{\prime }$
is smaller than the full solution (figure 5
a) by an order of magnitude. Furthermore, the latitude-dependent part of
$\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}^{\prime }$
(figure 5
d), which is the only part contributing to the gravity moments, is an order of magnitude smaller than
$\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}^{\prime }$
itself (figure 5
a). This explains why the contribution of
$\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}^{\prime }$
to the gravitational moments (figure 3, dashed line) is at least two orders of magnitude smaller than that of the non-null-space part of the solution. Overall, this analysis shows that this perturbation approach gives, to leading order, results that are very close to those of the CMS.
(a,b) The solution for
$\unicode[STIX]{x1D70C}^{\prime }$
(a) and the latitude-dependent part of the solution (b). (c,d) The solution for
$\unicode[STIX]{x1D6FF}\hat{\unicode[STIX]{x1D70C}}^{\prime }$
(c) and its latitude-dependent part (d). Unlike the eigenvectors (figure 2), the solution for the null space is large scale in both radius and latitude. It is negative everywhere (as a result of the need to compensate for
$\unicode[STIX]{x1D6FF}\hat{\unicode[STIX]{x1D70C}}^{\prime }$
, which is positive everywhere). A large-scale latitude-dependent structure exists, with the highest values close to the planet surface. All values are in
$\text{kg}~\text{m}^{-3}$
.
3.2 Analysis of solution and comparison with previous approximations
We now assess the contribution of each term on the right-hand side of the equation to the density solution (figure 4), and to the gravitational moments in particular. Since the equation is linear with respect to
$\unicode[STIX]{x1D70C}^{\prime }$
the analysis can be done by solving the equation when different terms are excluded. Following is a discussion of the TW approximation of (Kaspi et al.
Reference Kaspi, Hubbard, Showman and Flierl2010), and of the TGW solution of Zhang et al. (Reference Zhang, Kong and Schubert2015).
3.2.1 Spherically symmetric TW approximation
The simplest solution to (2.6), the TW approximation, is obtained when assuming that the static solution is spherically symmetric (Kaspi et al.
Reference Kaspi, Hubbard, Showman and Flierl2010), density and gravity as in figure 1(a,c) satisfying
$\unicode[STIX]{x1D70C}_{0}=\unicode[STIX]{x1D70C}_{0}^{m}$
,
$g_{0}^{(r)}=g_{0}^{(r)m}$
, and
$g_{0}^{(\unicode[STIX]{x1D703})}=0$
, and neglecting the gravity anomaly
$\boldsymbol{g}^{\prime }$
, so that
$$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\unicode[STIX]{x1D70C}=\unicode[STIX]{x1D70C}_{0}(r)+\unicode[STIX]{x1D70C}^{\prime }(r,\unicode[STIX]{x1D703}),\\ \boldsymbol{g}=g_{0}^{(r)}(r).\end{array}\right\} & & \displaystyle\end{eqnarray}$$
These assumptions reduce (2.6) to
$$\begin{eqnarray}\displaystyle 2\unicode[STIX]{x1D6FA}r\unicode[STIX]{x2202}_{z}(\unicode[STIX]{x1D70C}_{0}u)=g_{0}^{(r)}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}^{\prime }}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}, & & \displaystyle\end{eqnarray}$$
where the centrifugal terms drop under the background sphericity assumptions (see the discussion in § 4). The solution for the anomalous density can be simply found by integrating the right-hand side of (3.2) so that
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70C}^{\prime }(r,\unicode[STIX]{x1D703})=\int ^{\unicode[STIX]{x1D703}}\frac{2\unicode[STIX]{x1D6FA}r}{g_{0}^{(r)}(r)}\unicode[STIX]{x2202}_{z}(\unicode[STIX]{x1D70C}_{0}(r)u(r,\unicode[STIX]{x1D703}^{\prime }))\,\text{d}\unicode[STIX]{x1D703}^{\prime }+\tilde{\unicode[STIX]{x1D70C}}^{\prime }(r), & & \displaystyle\end{eqnarray}$$
where
$\tilde{\unicode[STIX]{x1D70C}}^{\prime }$
is an unknown integration coefficient that does not contribute to the gravitational moments (see (2.29)). Note that (3.2) is not the standard form of the TW equation, e.g. Vallis (Reference Vallis2006), since it includes
$\unicode[STIX]{x1D70C}_{0}$
on the left-hand side, and the right-hand side is not a purely baroclinic term; nonetheless, the two forms are equivalent (see the details in Kaspi et al.
Reference Kaspi, Davighi, Galanti and Hubbard2016). Solving (3.3) and calculating the gravity moments using (2.28) we can compare the solution with the CMS method (figure 6
a). The TW solution follows the full CMS solution with the ratio between the calculated moments being
$0.91,1.44,1.6,1.53,1.56,0.88$
for
$J_{2},J_{4},\ldots ,J_{12}$
, respectively. Note that in order to maintain the same framework, the equation was solved numerically using the same methodology as in the full case. Solving the equation using the method of Kaspi et al. (Reference Kaspi, Hubbard, Showman and Flierl2010), which is much more efficient numerically, gives the same results up to numerical roundoff.
The wind-induced gravitational moments from several limits of (2.28). (a) For barotropic winds, showing the FTW solution (blue), standard TW (green) and the TGW (grey). The CMS solution is shown in red for comparison. All solutions are quite similar, indicating that the simple TW approximation produces essentially the correct solution shown by the fuller approximations. (b) For a baroclinic case (with a wind-decay scale of 1000 km). Shown are the FTW solution (blue), standard TW (green) and TGW (grey) as well as a solution where the terms depending on
$\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{1}/\unicode[STIX]{x2202}r$
are included (see the text for a discussion of the TW+ case).
A variation on this case (Cao & Stevenson Reference Cao and Stevenson2016), is to allow the background density
$\unicode[STIX]{x1D70C}_{0}$
, as well as the gravity in the radial direction,
$g_{0}^{(r)}$
, to vary with latitude (figure 1
b,d). The gravity in the latitudinal direction is kept zero. The resulting equation is the same as (3.3), but with the background fields being a function of both radius and latitude,
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70C}^{\prime }(r,\unicode[STIX]{x1D703})=\int ^{\unicode[STIX]{x1D703}}\frac{2\unicode[STIX]{x1D6FA}r}{g_{0}^{(r)}(r,\unicode[STIX]{x1D703})}\unicode[STIX]{x2202}_{z}(\unicode[STIX]{x1D70C}_{0}(r,\unicode[STIX]{x1D703})u(r,\unicode[STIX]{x1D703}^{\prime }))\,\text{d}\unicode[STIX]{x1D703}^{\prime }+\tilde{\unicode[STIX]{x1D70C}}^{\prime }(r). & & \displaystyle\end{eqnarray}$$
The solution to this approximation is very similar to the above simplified TW balance (indistinguishable from the green line in figure 6), aside from some differences in the higher gravity moments, especially
$J_{12}$
as was also found by Cao & Stevenson (Reference Cao and Stevenson2016).
3.2.2 The TGW approximation
Next, we examine the contribution of the anomalous gravity to the solution, termed by Zhang et al. (Reference Zhang, Kong and Schubert2015) the TGW equation. They suggested that since the density perturbations
$\unicode[STIX]{x1D70C}^{\prime }$
result also in perturbations to the gravity field
$\boldsymbol{g}^{\prime }$
, these in turn affect the solution, and therefore need to be included in the balance. As in Zhang et al. (Reference Zhang, Kong and Schubert2015) we assume the background state to vary with radius only, so that
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70C}=\unicode[STIX]{x1D70C}_{0}(r)+\unicode[STIX]{x1D70C}^{\prime }(r,\unicode[STIX]{x1D703}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{g}=g_{0}^{(r)}(r)+\boldsymbol{g}^{\prime }(r,\unicode[STIX]{x1D703}). & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle -2\unicode[STIX]{x1D6FA}r\unicode[STIX]{x2202}_{z}(\unicode[STIX]{x1D70C}_{0}u)=-g_{0}^{(r)}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}^{\prime }}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}+r\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{0}}{\unicode[STIX]{x2202}r}{g^{\prime }}^{(\unicode[STIX]{x1D703})}. & & \displaystyle\end{eqnarray}$$
This equation cannot be easily integrated in
$\unicode[STIX]{x1D703}$
, and needs to be solved numerically (Zhang et al.
Reference Zhang, Kong and Schubert2015). The resulting gravity moments are shown in figure 6(a) (grey), together with the TW solution and the full perturbation method solution. The overall effect of the term added in (3.7) relative to the simplified TW is small. It is mostly apparent in
$J_{2}$
which increases by 54 %. The effect on higher moments, not calculated by Zhang et al. (Reference Zhang, Kong and Schubert2015), is much smaller. The small effect of the additional term in TGW approximation is already clear from the magnitude of the relevant term figure 4(e) (repeating figure 4 with a radially dependent background state gives a similar structure and magnitude to that shown in figure 4
a,c,e). Note that solving the equation with background fields that are both radially and latitudinally dependent shows similar results in the gravity moments.
3.3 The perturbation method in the more general case of baroclinic winds
So far we have focused only on barotropic cases since the CMS solution, which we used as our reference, may only be obtained for barotropic winds. The FTW can be used also to analyse baroclinic winds, which are considered in this section. Using the baroclinic winds described in § 2.4 (decay depth of wind is 1000 km), we repeat the above calculations of the density field and gravity moments. The gravitational moments for this case are shown in figure 4(b), and the individual terms in the FTW equation are shown in figure 7.
The solution for the gravitational moments shows that the TW (Kaspi et al.
Reference Kaspi, Hubbard, Showman and Flierl2010) and TGW (Zhang et al.
Reference Zhang, Kong and Schubert2015) solutions are again remarkably similar, apart from
$J_{2}$
. The solutions for both of these approximations are similar to the fuller FTW approximation, except for the moments
$J_{2}$
,
$J_{6}$
and
$J_{8}$
. In particular, the fuller solution for
$J_{8}$
is an order of magnitude smaller than both cruder approximations, underlining the importance of considering the additional physical effects included in this paper. In this case, there is a significant difference between TW (green) and the similar TGW (grey) on the one hand, and FTW (blue) on the other. The main reason for this difference are the two terms involving
$\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}^{\prime }/\unicode[STIX]{x2202}r$
. This is shown by the dash black curve denoted TW+, where we have used the TW solution plus the terms shown in figure 7(c,h), which both involve the radial derivative of the perturbation density. The importance of these terms is a direct consequence of the structure of the perturbation density solution, which tends to be strongly concentrated near the upper surface. This surface enhancement is not surprising given that the wind forcing decays rapidly away from the surface in this baroclinic case. This also implies that all terms in the equation tend to be more concentrated near the surface than in the barotropic case (compare figure 4 and figure 7). For this baroclinic case, the TGW term (figure 7
e) is negligible relative to most other terms considered here. Note also that the solution of the null space (§ 2.3) relies on the barotropically based polytropic equation, therefore some inconsistencies might arise due to that. This, however, should not affect much the solutions since the null-space contribution to the solution is small.
Solution of (2.6) for a baroclinic case (depth of winds is 1000 km) when all terms are included: (a) left-hand side term; (b) total right-hand side; (c–h) the six terms on the right-hand side. Values of
$\unicode[STIX]{x1D70C}^{\prime }$
are in
$\text{kg}~\text{m}^{-3}$
.
4 Discussion and conclusion
In the traditional approximation for terrestrial planets the centrifugal term is often merged with gravity in the momentum equation, by choosing the vertical direction to be that perpendicular to the planet’s geopotential surface and defining an effective gravity. This is then traditionally followed by approximating the planet as a sphere, so that the vertical direction coincides with the radial direction, and thus effectively neglecting the horizontal component of the centrifugal term. It is important to note that this centrifugal term is not smaller than the Coriolis term even for the Earth case, but because of the nearly spherical shape of Earth, this approximation allows trading a large dynamical component in the momentum balance with a small geometric error (Vallis Reference Vallis2006, § 2.2). This approximation has proven to hold well for Earth and other terrestrial planets. On the giant planets, the oblateness is not small (6.5 % and 9.8 % on Jupiter and Saturn, respectively, compared with 0.3 % on Earth). Therefore, the contribution of the centrifugal (fifth and sixth terms on the right-hand side of (2.6)) and self-gravitation terms (third and fourth terms on the right-hand side of (2.6)) can potentially lead to significant contributions to the momentum balance, and therefore may alter thermal-wind balance as well. The goal of this study is to assess the importance of these terms in a fluid planet.
By solving numerically the full second-order momentum equation in our perturbation approach, which includes the original thermal-wind balance terms (left-hand side term and second term on right-hand side in (2.6)), self-gravity terms, centrifugal terms, and other non-spherical contributions (first term on right-hand side of (2.6)), we have shown that the original TW balance is still the leading order. In the barotropic limit, the TW results, with the various higher-order contributions, are systematically compared with results from the oblate CMS model. In the context of recent studies that argue that additional terms are important in the balance for calculating the gravitational moments (Zhang et al.
Reference Zhang, Kong and Schubert2015; Cao & Stevenson Reference Cao and Stevenson2016), we show that to leading-order these terms are negligible, and have a small contribution to the gravity moments. Consistently with Zhang et al. (Reference Zhang, Kong and Schubert2015), we find that the self-gravity term (TGW) increases the value of
$J_{2}$
, though it does not bring the TW
$J_{2}$
closer to the CMS result. This terms has a negligible contribution to all higher harmonics, which were not discussed in Zhang et al. (Reference Zhang, Kong and Schubert2015).
We conclude therefore that while more complete solutions are possible, as we do in this study, the traditional thermal-wind gives a very good approximation to the balance between the wind shear and the density gradients, and integrating it gives a very good approximation to the gravity moments. In particular, taking into account the accuracy of the Juno and Cassini measurements, this gives an excellent approximation. Quantitatively, for the barotropic cases, its results differ by at most a factor of 1.6 compared with the full solution. This difference is small considering the other uncertainties of the interior flow. For the baroclinic case, where wind structures decay rapidly near the surface, terms involving the radial derivative of the perturbation density become more important for calculating the gravity moments.
The main advantage of using the TW model compared with FTW is numerical. While the TW equation (3.3) allows for local calculation of the density from the wind, the FTW equation (2.6) is an integro-differential equation that needs to be solved globally. It is mostly complicated from the need to integrate the dynamical density (
$\unicode[STIX]{x1D70C}^{\prime }$
) globally to calculate the dynamical self-gravity (
$\boldsymbol{g}^{\prime }$
). The TW approximation allows therefore using much higher resolution, which is necessary for resolving the high-order moments. As a consequence of the simplicity of the TW model, more sophisticated and numerically demanding methods can be applied in order to find the best-matching wind field given the gravity measurements (Galanti & Kaspi Reference Galanti and Kaspi2016a
; Galanti & Kaspi Reference Galanti and Kaspi2016b
). The invertibility of the solution using the TW model is a major advantage for the upcoming analysis of the Juno and Cassini data. Given the extremely small contribution of the null space to the overall solution, we expect that the more complete FTW model would also be invertible, still the computational challenge involved is much greater.
In summary, deciphering the effect of the atmospheric and internal flows from the measured gravity spectrum of Jupiter and Saturn provides a major challenge. The methods suggested to date have been either limited to barotropic cases (e.g. Hubbard Reference Hubbard1982, Reference Hubbard1999, Reference Hubbard2012; Kong et al. Reference Kong, Zhang and Schubert2012; Hubbard et al. Reference Hubbard, Schubert, Kong and Zhang2014), or approximations limited to spherical symmetry or partial solutions (e.g. Kaspi et al. Reference Kaspi, Hubbard, Showman and Flierl2010; Zhang et al. Reference Zhang, Kong and Schubert2015; Cao & Stevenson Reference Cao and Stevenson2016). Here, we have developed a self-consistent perturbation approach to the TW balance that incorporates all physical effects, including the effects of oblateness on the dynamics and the gravity perturbation induced by the flow itself. The full self-consistent perturbation approach to the TW balance considered here allows us to objectively examine the role of different physical processes, allows obtaining and even more accurate approximation by proceeding to higher-order perturbation corrections (appendix B), and allows the interpretation of the expected Juno and Cassini observations in a more complete way than was possible in previous approaches, thus maximizing the benefits of these observations.
Acknowledgements
We thank the Juno Science Team Interiors Working Group for valuable discussions. E.T. is funded by the NSF Physical Oceanography program, grant OCE-1535800, and thanks the Weizmann Institute of Science (WIS) for its hospitality during parts of this work. Y.K. and E.G. acknowledge support from the Israeli Ministry of Science (grant 45-851-641), the Minerva foundation with funding from the Federal German Ministry of Education and Research, and from the WIS Helen Kimmel Center for Planetary Science.
Appendix A. Higher-order perturbation equations
We write here the perturbation equations to the second order to demonstrate how the results of our approach can be made more accurate if needed. The momentum equation (2.2) is
$$\begin{eqnarray}\displaystyle 2\unicode[STIX]{x1D734}\times (\unicode[STIX]{x1D70C}\boldsymbol{u})=-\unicode[STIX]{x1D735}p-\unicode[STIX]{x1D70C}\boldsymbol{g}-\unicode[STIX]{x1D70C}\unicode[STIX]{x1D734}\times \unicode[STIX]{x1D734}\times \boldsymbol{r}. & & \displaystyle\end{eqnarray}$$
Writing the density as
$\unicode[STIX]{x1D70C}=\unicode[STIX]{x1D70C}_{0}+\unicode[STIX]{x1D70C}_{1}+\unicode[STIX]{x1D70C}_{2}$
, and substituting into the above equation, treating
$\unicode[STIX]{x1D70C}_{1}$
as an order
$\unicode[STIX]{x1D716}$
correction and
$\unicode[STIX]{x1D70C}_{2}$
as an order
$\unicode[STIX]{x1D716}^{2}$
correction, we can separate the different orders to find equations for each correction order. Note that in the previous sections we denote
$\unicode[STIX]{x1D70C}_{1}$
as
$\unicode[STIX]{x1D70C}^{\prime }$
. We view, as defined earlier, the zeroth-order balance as the balance without the effects of the winds, so that our zeroth-order equation is
$$\begin{eqnarray}\displaystyle 0=-\unicode[STIX]{x1D735}p_{0}-\unicode[STIX]{x1D70C}_{o}\boldsymbol{g}_{0}-\unicode[STIX]{x1D70C}_{0}\unicode[STIX]{x1D734}\times \unicode[STIX]{x1D734}\times \boldsymbol{r}. & & \displaystyle\end{eqnarray}$$
The winds enter at the first order, where the equation is
$$\begin{eqnarray}\displaystyle 2\unicode[STIX]{x1D734}\times (\unicode[STIX]{x1D70C}_{0}\boldsymbol{u})=-\unicode[STIX]{x1D735}p_{1}-\unicode[STIX]{x1D70C}_{o}\boldsymbol{g}_{1}-\unicode[STIX]{x1D70C}_{1}\boldsymbol{g}_{0}-\unicode[STIX]{x1D70C}_{1}\unicode[STIX]{x1D734}\times \unicode[STIX]{x1D734}\times \boldsymbol{r}, & & \displaystyle\end{eqnarray}$$
and the second-order correction is then obtained by solving
$$\begin{eqnarray}\displaystyle 2\unicode[STIX]{x1D734}\times (\unicode[STIX]{x1D70C}_{1}\boldsymbol{u})+\text{}\unicode[STIX]{x1D70C}_{1}\boldsymbol{g}_{1}\text{}=-\unicode[STIX]{x1D735}p_{2}-\unicode[STIX]{x1D70C}_{o}\boldsymbol{g}_{2}-\unicode[STIX]{x1D70C}_{2}\boldsymbol{g}_{0}-\unicode[STIX]{x1D70C}_{2}\unicode[STIX]{x1D734}\times \unicode[STIX]{x1D734}\times \boldsymbol{r}. & & \displaystyle\end{eqnarray}$$
In the second-order perturbation equation we moved all terms that are known from previous orders to the left-hand side. This equation is again solved by taking a curl and then taking care of the integration constant (null space solution) if needed. Because the second-order equation is generally similar to the first-order equation, its numerical solution follows the same approach and should not pose significant additional difficulties. While this procedure should be formally done in a non-dimensional form, we present here the dimensional equations for clarity. The small parameter in this expansion when it is done in non-dimensional form is the Rossby number
$Ro=U/(\unicode[STIX]{x1D6FA}L)$
, where
$U$
is a scale for the velocity and
$L$
the horizontal length scale of the relevant motions.
Appendix B. Integration constant and null space
The objective of this appendix is to show how the integration constant encountered in the TW approach (e.g. Kaspi et al. Reference Kaspi, Hubbard, Showman and Flierl2010; Zhang et al. Reference Zhang, Kong and Schubert2015) can be calculated. This is meant to aid the understanding of our approach to solving for the null space of the more general solution considered in the main text. The main message of this appendix is that the integration constant may be determined by adding the missing physics of the polytropic equation, and by requiring the total mass of the density perturbation to vanish. The momentum equation for the simple example is
$$\begin{eqnarray}0=-\unicode[STIX]{x1D735}p-\boldsymbol{g}\unicode[STIX]{x1D70C}+\unicode[STIX]{x1D6FA}^{2}\unicode[STIX]{x1D70C}r\cos \unicode[STIX]{x1D703}\hat{\boldsymbol{r}}_{\bot }.\end{eqnarray}$$
We treat the effects of rotation as a perturbation and the leading-order balance is then hydrostatic, with the corresponding fields being a function of radius alone (
$\unicode[STIX]{x1D70C}_{0}=\unicode[STIX]{x1D70C}_{0}(r)$
),
$$\begin{eqnarray}0=-\unicode[STIX]{x1D735}p_{0}-\boldsymbol{g}_{0}\unicode[STIX]{x1D70C}_{0}.\end{eqnarray}$$
The next order contains the deviations in density and gravity due to rotation,
$$\begin{eqnarray}0=-\unicode[STIX]{x1D735}p_{1}-\boldsymbol{g}_{1}\unicode[STIX]{x1D70C}_{0}-\boldsymbol{g}_{0}\unicode[STIX]{x1D70C}_{1}+\unicode[STIX]{x1D6FA}^{2}\unicode[STIX]{x1D70C}_{0}r\cos \unicode[STIX]{x1D703}\hat{\boldsymbol{r}}_{\bot }.\end{eqnarray}$$
Taking the curl,
$$\begin{eqnarray}\displaystyle 0=-g_{1}^{(\unicode[STIX]{x1D703})}r\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{0}}{\unicode[STIX]{x2202}r}+g_{0}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{1}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}-\unicode[STIX]{x1D6FA}^{2}r^{2}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{0}}{\unicode[STIX]{x2202}r}\cos \unicode[STIX]{x1D703}\sin \unicode[STIX]{x1D703}, & & \displaystyle\end{eqnarray}$$
substituting the expression for the gravity fields, and integrating over
$\unicode[STIX]{x1D703}$
,
$$\begin{eqnarray}\displaystyle -\frac{1}{2}\unicode[STIX]{x1D6FA}^{2}r^{2}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{0}}{\unicode[STIX]{x2202}r}\cos ^{2}\unicode[STIX]{x1D703}=g_{0}\unicode[STIX]{x1D70C}_{1}+2\unicode[STIX]{x03C0}G\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}_{0}}{\unicode[STIX]{x2202}r}\iint \frac{\unicode[STIX]{x1D70C}_{\mathbf{1}}(r^{\prime },\unicode[STIX]{x1D703}^{\prime })r^{\prime \mathbf{2}}\cos (\unicode[STIX]{x1D703}^{\prime })}{\langle |\boldsymbol{r}-\boldsymbol{r}^{\prime }|\rangle }\,\text{d}r^{\prime }\,\text{d}\unicode[STIX]{x1D703}^{\prime }+C(r),\quad & & \displaystyle\end{eqnarray}$$
where
$C(r)$
is the unknown integration constant to be solved for, and the left-hand side represents rotation effects due to the already known zeroth-order solution. In the following, we solve for the perturbation density by expressing the density as a sum of two terms. First,
$\hat{\unicode[STIX]{x1D70C}}_{1}(r,\unicode[STIX]{x1D703})$
, that solves the above equation with the integration constant set to zero and, second,
$\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}_{1}(r)$
that satisfies the same equation with a zero on the left-hand side and with the integration constant. Together,
$\unicode[STIX]{x1D70C}_{1}=\hat{\unicode[STIX]{x1D70C}}_{1}(r,\unicode[STIX]{x1D703})+\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}_{1}(r)$
solves the full equation (B 5).
Given that we took the curl of the momentum equation,
$\hat{\unicode[STIX]{x1D70C}}_{1}$
which solves the above equation without
$C(r)$
satisfies the original momentum equation up to a gradient of some function which we may write as
$\hat{p}_{1}+\unicode[STIX]{x1D6FF}p_{1}$
,
$$\begin{eqnarray}\displaystyle 0=-\unicode[STIX]{x1D735}(p_{1}+\unicode[STIX]{x1D6FF}p_{1})+\unicode[STIX]{x1D70C}_{0}\hat{\boldsymbol{g}}_{1}+\hat{\unicode[STIX]{x1D70C}}_{1}\boldsymbol{g}_{0}+\unicode[STIX]{x1D70C}_{0}\unicode[STIX]{x1D6FA}^{2}r\cos \unicode[STIX]{x1D703}\hat{\boldsymbol{r}}_{\bot }. & & \displaystyle\end{eqnarray}$$
Next, take the difference between the momentum equation for
$\unicode[STIX]{x1D70C}_{1}$
and for
$\hat{\unicode[STIX]{x1D70C}}_{1}$
, we find
$$\begin{eqnarray}0=-\unicode[STIX]{x1D735}\unicode[STIX]{x1D6FF}p_{1}+\unicode[STIX]{x1D70C}_{0}\unicode[STIX]{x1D739}\boldsymbol{g}_{1}+\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}_{1}\boldsymbol{g}_{0}.\end{eqnarray}$$
The second and third terms in (B 7) are functions of radius
$r$
only. Therefore, we expect the pressure term
$-\unicode[STIX]{x1D735}\unicode[STIX]{x1D6FF}p_{1}$
to also be a function of the radius only. Furthermore, assuming now that the polytropic relation holds for the perturbation pressure and density, we have
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FF}p_{1}=K2\unicode[STIX]{x1D70C}_{0}\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}_{1}, & & \displaystyle\end{eqnarray}$$
which leads to an equation for the unknown perturbation density,
$$\begin{eqnarray}\displaystyle 0=-\unicode[STIX]{x1D735}(K2\unicode[STIX]{x1D70C}_{0}\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}_{1})+\unicode[STIX]{x1D70C}_{0}\unicode[STIX]{x1D739}\boldsymbol{g}_{1}+\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}_{1}\boldsymbol{g}_{0}. & & \displaystyle\end{eqnarray}$$
Substituting the expression for the gravity,
$$\begin{eqnarray}\displaystyle 0=-\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}(Kn2\unicode[STIX]{x1D70C}\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}_{1})+\unicode[STIX]{x1D70C}_{0}2\unicode[STIX]{x03C0}G\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}\iint \left\langle \frac{1}{|\boldsymbol{r}-\boldsymbol{r}^{\prime }|}\right\rangle \unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}_{1}(r^{\prime })r^{\prime \mathbf{2}}\cos (\unicode[STIX]{x1D703}^{\prime })\,\text{d}r^{\prime }\,\text{d}\unicode[STIX]{x1D703}^{\prime }+\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}_{1}g_{0}, & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$
and substituting the zeroth-order solutions for
$\unicode[STIX]{x1D70C}_{0},$
$p_{0}$
and
$g_{0}$
(Zhang et al.
Reference Zhang, Kong and Schubert2015, see notation there), and after some more rearrangement and integration in
$\unicode[STIX]{x1D709}$
we obtain
$$\begin{eqnarray}\displaystyle 0=-\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}_{1}+\frac{1}{2}\iint \left\langle \frac{1}{|\unicode[STIX]{x1D709}-\unicode[STIX]{x1D709}^{\prime }|}\right\rangle \unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}_{1}(\unicode[STIX]{x1D709}^{\prime })\unicode[STIX]{x1D709}^{\prime 2}\cos (\unicode[STIX]{x1D703})\,\text{d}\unicode[STIX]{x1D709}^{\prime }\,\text{d}\unicode[STIX]{x1D703}+D, & & \displaystyle\end{eqnarray}$$
where
$D$
is the constant of integration. Use the expansion for the average over longitude (Zhang et al.
Reference Zhang, Kong and Schubert2015), and the fact that the integral over all Legendre polynomials but the first vanish,
$$\begin{eqnarray}\displaystyle 0=-\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}_{1}(\unicode[STIX]{x1D709})+\int _{0}^{\unicode[STIX]{x03C0}}f_{0}(\unicode[STIX]{x1D709},\unicode[STIX]{x1D709}^{\prime })\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}_{1}(\unicode[STIX]{x1D709}^{\prime })\unicode[STIX]{x1D709}^{\prime 2}\,\text{d}\unicode[STIX]{x1D709}^{\prime }+D, & & \displaystyle\end{eqnarray}$$
where
$$\begin{eqnarray}\displaystyle f_{0}(\unicode[STIX]{x1D709},\unicode[STIX]{x1D709}^{\prime })=\left\{\begin{array}{@{}ll@{}}\displaystyle \frac{1}{\unicode[STIX]{x1D709}}\quad & \unicode[STIX]{x1D709}^{\prime }\leqslant \unicode[STIX]{x1D709}\\ \displaystyle \frac{1}{\unicode[STIX]{x1D709}^{\prime }}\quad & \unicode[STIX]{x1D709}^{\prime }>\unicode[STIX]{x1D709}.\end{array}\right. & & \displaystyle\end{eqnarray}$$
A solution may be found by assuming a Frobenius-form solution,
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}_{1}(\unicode[STIX]{x1D709})=\mathop{\sum }_{m=0}^{\infty }a_{m}\unicode[STIX]{x1D709}^{m+r}. & & \displaystyle\end{eqnarray}$$
Now find what is
$a_{m}$
. Write (B 12) explicitly
$$\begin{eqnarray}\displaystyle 0=-\mathop{\sum }_{m=0}^{\infty }a_{m}\unicode[STIX]{x1D709}^{m+r}+\mathop{\sum }_{m=0}^{\infty }\int _{0}^{\unicode[STIX]{x1D709}}\frac{1}{\unicode[STIX]{x1D709}}a_{m}\unicode[STIX]{x1D709}^{\prime m+r}\unicode[STIX]{x1D709}^{\prime \mathbf{2}}\,\text{d}\unicode[STIX]{x1D709}^{\prime }+\mathop{\sum }_{m=0}^{\infty }\int _{\unicode[STIX]{x1D709}}^{\unicode[STIX]{x03C0}}\frac{1}{\unicode[STIX]{x1D709}^{\prime }}a_{m}\unicode[STIX]{x1D709}^{\prime m+r}\unicode[STIX]{x1D709}^{\prime \mathbf{2}}\,\text{d}\unicode[STIX]{x1D709}^{\prime }+D,\quad & & \displaystyle\end{eqnarray}$$
and performing the integral and collecting powers, while also defining
$a_{-1}=a_{-2}=0$
,
$$\begin{eqnarray}\displaystyle 0=\mathop{\sum }_{m=0}^{\infty }\left(-a_{m}+\frac{1}{(m+1+r)}a_{m-2}-\frac{1}{m+r}a_{m-2}\right)\unicode[STIX]{x1D709}^{m+r}+\mathop{\sum }_{m=0}^{\infty }\frac{1}{m+2+r}a_{m}\unicode[STIX]{x03C0}^{m+2+r}+D. & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$
The coefficient of
$\unicode[STIX]{x1D709}^{m+r}$
should vanish for all
$m$
, giving us the recursion relation. Multiplying by
$(m+r)(m+r+1)$
and considering the
$m=0$
and
$m=1$
cases, we find that there are two possible solutions for the Frobenius power,
$r=0$
or
$r=-1$
. The first leads to the recursion relation,
$$\begin{eqnarray}\displaystyle a_{m}=-\frac{1}{m(m+1)}a_{m-2}. & & \displaystyle\end{eqnarray}$$
We term this solution
$F_{1}(\unicode[STIX]{x1D709}).$
The second solution, with
$r=-1$
, leads to
$$\begin{eqnarray}\displaystyle a_{m}=-\frac{1}{(m-1)m}a_{m-2}, & & \displaystyle\end{eqnarray}$$
which is the recursion relation for cosine. Therefore,
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}_{1}=C_{1}F_{1}+C_{2}\frac{\cos (\unicode[STIX]{x1D709})}{\unicode[STIX]{x1D709}}. & & \displaystyle\end{eqnarray}$$
The
$\cos (\unicode[STIX]{x1D709})/\unicode[STIX]{x1D709}$
solution is not physical because it diverges at
$\unicode[STIX]{x1D709}=0$
and we conclude that
$C_{2}=0$
. Considering the coefficients of
$\unicode[STIX]{x1D709}^{0}$
in (B 16) we find
$a_{0}$
in terms of the unknown integration constant
$D$
. We may then use
$D$
to satisfy the constraint that volume integral over the total perturbation
$\unicode[STIX]{x1D70C}_{1}(\unicode[STIX]{x1D709},\unicode[STIX]{x1D703})+\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D70C}_{1}(\unicode[STIX]{x1D709})$
should vanish, as done in the manuscript for the fuller problem considered there.
