A Linear Approximation for the Effect of Cylindrical Differential Rotation on Gravitational Moments: Application to the Non-Unique Interpretation of Saturn's Gravity
Yayaati Chachan David J. Stevenson

Evolution from B_2g nematics to B_1g nematics in heavily
hole-doped iron-based superconductors
    Roser Valentí
    December 30, 2023
==============================================================================================




The higher order gravitational moments of a differentially rotating planet are greatly affected or even dominated by the near surface differential flow. Unlike the contribution from rigid body rotation, this part of the gravity field can be well estimated by linear theory when the differential rotation is on cylinders and the corresponding gravity field arises from the higher order moments directly introduced by that flow. In the context of an n = 1 polytrope, we derive approximate analytical formulas for the gravity moments. We find that Δ J_2n is typically at most a few times (-1)^n+1  a   q   d^5/2, where a is the amplitude of the differential rotation (as a fraction of the background rigid body rotation), q = Ω^2 R^3 / G M is the usual dimensionless measure of rotation for the planet (mass M, radius R), and d << 1 is the characteristic depth of the flow as a fraction of the planetary radius. Applied to Saturn, with a set by the observed surface wind amplitude, we find first that the observed signs of the Δ J_2n are a trivial consequence of the definition of the corresponding Legendre polynomials, but the Cassini observations can not be explained by a simple exponentially decaying flow and instead require a substantial retrograde flow at depth and a larger d than the simple scaling suggests. This is consistent with the results reported by <cit.>. However, there is no fluid dynamical requirement that the flows observed in the atmosphere are a guide to the flows thousands of km deeper. We explore a wide range of flow depths and amplitudes which yield values for Δ J_n that are acceptable within the error estimates and thus highlight the inherent non-uniqueness of inferences made from the higher order gravity moments.




§ INTRODUCTION
 

The external gravity field of an isolated, uniformly rotating, hydrostatic planet is expressible in the form


    V_ext (r, θ) = G M/r[ 1 - ∑_n=1^∞(a/r)^2n  J_2n P_2n(cos θ) ]


where J_n, the gravitational moments are given by:


    J_n = - 1/M a^n∫_R^3ρ (r) r^n P_n (cos θ) d^3r


where M is planet mass, a is planet radius (equatorial or volume averaged), P_n is the order n Legendre polynomial, ρ is the density, G is the gravitational constant, and θ is the colatitude. The dimensionless coefficients of this expansion are called gravitational moments and denoted J_2n where n is a positive integer. They exist because the planet rotates. The expansion contains only even moments because the rotational perturbation is quadratic in rotation rate and therefore does not distinguish between Northern and Southern hemispheres. In linear response theory, J_2 is proportional to q = Ω^2 R^3 / G M and all higher J's are zero, where Ω is the rotation rate, R is the planet's radius, M is the planet's mass, and G is the gravitational constant. This assumes rigid body rotation. However, q is not very small for giant planets and so the linear response is far from adequate, especially given the high accuracy with which the gravity field can be measured by orbiting spacecraft such as Juno and Cassini. The higher J's arise explicitly from non-linearity; that is, the centrifugal potential is exactly of degree 2, but the change in shape and resulting self gravity introduces the higher order terms so that J_2n is approximately proportional to q^n <cit.>. Moreover, even low order moments are affected by higher order corrections. Although the J's are linearly related to the density anomalies of the same harmonic degree, the underlying non-linearity makes the theory complicated, even for simple models of planetary structure, and it can be difficult to develop an intuition for the dependence of J's on features of that structure, which is, after all, the whole point of making the precise measurements. However, a great deal of successful effort has been put into forward modeling of the J's from physical models of the planet's interior <cit.>. This is necessarily non-unique (you cannot invert the J's to determine planet structure except in the context of models with only a few degrees of freedom) but it has provided us with a great deal of understanding of fluid planets.

However, fluid bodies do not rotate as rigid bodies. This is self-evident in the zonal flows observed in their atmospheres. In general, this flow need not be purely cylindrical, a special case of the Taylor-Proudman theorem (i.e. with no variation of the East-West flow along an axis parallel to the rotation axis). Deviations from the Taylor-Proudman theorem in a fluid with large density variations but no MHD effects arise primarily from latitudinal entropy gradients but also potentially from the Reynolds stress <cit.>. The advantage of assuming cylindrical differential rotation is that it allows for the definition of a centrifugal potential, the gradient of which provides the acceleration and thus the wind strength at each cylindrical radius. At first sight, it might seem that even this simple kind of differential rotation only exacerbates the difficulty of developing an intuitive or approximate understanding of the J's. However, there is a sense in which the effect of differential rotation may actually be easier to understand than the part due to rigid body rotation. The reason is that the centrifugal potential is no longer just degree 2, but contains higher order moments, much higher if the differential rotation is largely confined to a thin shell. One can then compute the gravitational response to that part alone using linear theory. Moreover, we can then approximate the planet as spherical. This will be an adequate approximation (with a fractional error of order q) provided the differential rotation is sufficiently small relative to the total rotation. This is the approach presented in this paper. It assumes that one can isolate the contribution to the J's that arises from the differential rotation from the often larger or comparable part that arises from the rigid part of the rotation. This can often be done because the higher order J's arising from rigid body rotation are usually confined to a narrow range of possible values based on interior models that fit the lower J's [This is predicated on the standard view that we actually know the structure of the outermost ten percent or so of the planet very well: an adiabatic hydrogen-dominated region for Jupiter and Saturn.] <cit.>. Moreover, the emphasis on great precision that usually accompanies analysis of the interior models need no longer apply – one might well be content with getting the differential rotation of a planet to just tens of percent, especially if the inference is necessarily non-unique, but also because the observational data can have substantial errors at the higher order moments.

This paper is motivated by a desire to develop a better intuition about what it is that determines the differential rotation contribution to the J's, more specifically how the amplitude and patterns of their contribution to the J's is set by the amplitude, depth, and shape of the differential flow profile. It focuses on the high order J's where the effect of differential rotation is more readily discerned in the data.

In keeping with our emphasis on transparency and simplicity, we use an n = 1 polytropic model (pressure everywhere proportional to the square of the density). This is not essential to the idea of our approach but it leads to a linear equation for hydrostatic equilibrium and enables the derivation of analytical formulas. As we shall explain, it is not necessary that the polytrope apply everywhere; models with a separate but small core can also be described in this way. The biggest error in this assumption lies in the breakdown of that equation of state as one approaches the surface, but this is still a smaller error than the error arising from assuming sphericity.



In Section 2, we describe our model and derive the main results for an arbitrary choice of functional form for the cylindrical differential rotation. In Section 3, we show how these results can yield simple scaling laws and simple analytical formulas for the case of differential flow that is exponentially declining with depth. We also show how to estimate the corrections that arise from planetary oblateness (i.e. corrections of order q). In Section 4, we apply our procedure to the Cassini observations of Saturn and show how we can recover flows similar to those suggested by more precise calculation, but also how they can be quite non-unique. We conclude with a discussion of what this approach tells us about the general features of gravitational effects caused by differential rotation.



§ THEORETICAL MODEL
 

If the angular velocity within the planet depends only on the cylindrical radius r_c = r sinθ, then the sole dependence of flow properties on r_c allows us to write their effect on the hydrostatic equilibrium of the planet as a gradient of a potential V_rot:


    1/ρ∇ P = ∇ (V + V_rot)


where ρ is the density, P is the pressure, and V is the gravitational potential. Figure <ref> shows an example of the differential flows considered in this work. All references to the depth of a flow apply to the cylindrical radius. We need an equation of state to relate the pressure to the density. This hydrostatic equation is difficult to solve for an arbitrary equation of state but it has a particularly simple form for one with polytropic index n = 1, i.e., P = K ρ^2. Here, K is a constant of proportionality that is mainly determined by the properties of the degenerate electron gas but is also somewhat affected by the planet's entropy. Using this polytropic equation, the hydrostatic equilibrium simplifies to: 


    ∇ (2 K ρ -V - V_rot) = 0


This can only be satisfied if 2 K ρ = V + V_rot. We assume that the contribution of differential rotation to the rotational potential is small compared to the background solid-body rotation potential. We normalize the potential to the monopole, G M / R, so we can write the rotational potential as:


    V_rot = q  ∫_0^s (1 + ϵ(s'))^2 s' ds' 
    ≈ q  ∫_0^s (1 + 2ϵ(s')) s' ds'

where s = r_c / R is the cylindrical radius normalised to the planet's radius, ϵ (s) is the differential flow profile as a function of s, the integrand is Taylor-expanded to first order under the assumption that ϵ << 1, and q is given by


    q = Ω^2 R^3 / G M

Here, G is the gravitational constant, M is the planet's mass, R its radius, and Ω is the background solid body rotation. The first term is simply the solid body rotation term. The second term, linear in ϵ, is henceforth denoted δ V_rot. The Poisson equation relates the Laplacian of the potential to the density and using the above relationship between density and potential, we obtain:


    ∇^2 δ V = - 4 π G ρ = - k^2 (δ V + δ V_rot)

where k^2 = (2π G) / K, δ V is only the self gravity potential arising from δ V_rot. We can solve this by writing δ V = V_h + V_p as a sum of a general (homogeneous) and a particular solution. The general solution is the solution to the homogeneous Helmholtz equation and is given by:


    V_h = ∑_Even   n  A_n   j_n(kr)   P_n(cos θ)

where j_n are the spherical Bessel's functions, P_n are Legendre Polynomials of order n, and A_n are the calculable coefficients. To calculate the particular solution V_p, we make an additional assumption that the differential rotation is significant near the surface of the planet and that it decays away rapidly inside the planet. This is a reasonable assumption for Jupiter, perhaps less so for Saturn. In this approximation, the curvature terms in the cylindrical Laplacian can be ignored:


    d^2 V_p/d r_c^2 = - k^2 δ V_rot


Solving for V_p, we get an expression which is a function of cylindrical radius r_c. However, the homogeneous solution (V_h) and the external potential (V_ext obtained from solving ∇^2 V_ext =0) are normally expressed as a function of the spherical coordinates. We therefore need to project V_p from cylindrical to spherical coordinates to match the boundary condition for the potential at the planet's surface. This can be done using the Legendre polynomials which form a complete basis set for axisymmetric functions of the spherical coordinates:


    V_p (r sinθ) = ∑_Even   n  V_pn(r)   P_n (cosθ) 
    
    V_pn(r) = 2n+1/2∫_-1^1 V_p(r sinθ)   P_n(cosθ)  d(cosθ)


We want to calculate the deviation in gravitational moments due to a small shallow flow and we can do that by applying boundary conditions to the potential at the planetary surface. We refer to these as Δ J_n where the J's are defined by Equation <ref>. We assume that this surface can be approximated by a sphere (we deal in Section <ref> with the error introduced by this). The potential and its derivative must be continuous at the surface which gives us two equations that relate the quantities V_pn, A_n j_n (kr), and Δ J_n at r = R and their derivatives. 


    V_pn (R) + A_n j_n (kR) = - Δ J_n 
    
    R V_pn' + k R A_n j_n'(kR) = (n+1) Δ J_n

Here, the prime ' refers to differentiation with respect to r (so j_n'(kR) =  d j_n(kr)/ dr evaluated at r = R). Eliminating A_n, we obtain the following expression for Δ J_n:


    Δ J_n = R V'_pn(R) j_n(kR) - kR V_pn(R) j'_n(kR)/(n+1)j_n(kR) + kR j'_n(kR)


For a planet with no core, k R = π. For a planet with a high density (rock or ice) core, k R = π (1-δ) where δ is the ratio of core to planet mass. However, this correction is small for the models we consider for the Δ J_n's. In the next section, we demonstrate the utility and power of this model that arises out of its simplicity.



§ ANALYTICAL SOLUTIONS
 
We use the theoretical framework established in the previous section to do analytic calculations for the gravitational moments due to differential flow in the limit of very shallow flow. Although, we do not use this analytic solution directly to fit Saturn's higher order gravitational moments, it is still instructive to perform this exercise to gain an understanding of how the deviations in gravitational moments relate to the properties of differential flow.



 §.§ Limit of Shallow Flow

Let us assume that the differential flow profile can be represented by a steeply declining exponential which is only a function of the cylindrical radius of a fluid planet:


    ϵ (s) = a   e^-b   (1-s)

Here, a and b are constants that characterize the amplitude and the depth of the flow, and s = r_c / R is again the cylindrical radius as a fraction of the equatorial radius. We can integrate this function analytically to obtain V_rot and calculate V_p from the differential equation <ref>, with the boundary condition that V_p and its derivative vanish at s = 0. To proceed further with the analytic calculation, we need to project V_p onto spherical coordinates using the Legendre Polynomials (equation <ref>). For very shallow flow, we can use a Taylor series expansion of the Legendre Polynomials up to a certain order (that we choose) and calculate the integrals analytically. We choose to expand the even Legendre Polynomials to second order in cos θ around θ = π/2:


    P_n (cos θ) = n  n/2(-1)^n/2/2^n[ 1 - n   (n+1)/2 (cosθ)^2 ]

where n  n/2 is the binomial coefficient. With this approximation, we can project V_p onto a basis set of Legendre polynomials with s = r  sinθ / R. We then have an expression for V_pn from equation <ref>. Given that we calculate this integral in the approximation of shallow flow, we can neglect all terms that do not contain the large exponential term in the expression for V_p. To calculate V_pn, we resort to the saddle point method (Appendix A) which works well for strongly peaked functions such as the exponential in the expression for V_p. Evaluating V_pn and its derivative at r = R = 1 (fractional radius in our calculations) and substituting it in our expression for Δ J_n, we obtain (see Appendix B for a detailed derivation):


    Δ J_n = - 2n + 1/2n  n/22 a    q  π^2/b^3(-1)^n/2/2^n√(2 π/b) ×
    [ ( -5/2 + b + 3/2b) + n (n+1)/4b( 7/2 -b - 9/2b) ] j_n(π) - π[ (1 - 3/b) ( 1 - n(n+1)/4b) ] j'_n(π)/(n+1)j_n(π) + π j'_n(π)




This result uses the saddle point method expansion (Appendix A) to first order in O(1/b) (constant term + first order term) as the contributions from higher orders are negligible. If we consider only the lowest order term then we find that


    Δ J_n = C_n (-1)^n/2+1  a   q   b^-5/2


where C_6 = 8.16, C_8 = 6.98, and C_10 = 6.22. However, these asymptotic (very small d = 1/b) limiting forms are not sufficiently accurate for the case of Saturn, as we shall discuss. The calculations from this analytical model are compared with the numerical calculations for a shallow flow that decays exponentially in Figure <ref>. The agreement between the analytical model and the numerical calculation improves with the inclusion of the O (1/b) term as well as for high b (increasingly shallow flow). The flow depth at which agreement becomes good is also a function of the order of gravitational moment, as expected. This is simply because cosθ changes a lot as sinθ changes by a small amount, when θ is near π/2. Consequently, the Legendre polynomials of high order n vary rapidly in the region of interest where the flow is non-negligible. The relationship between the flow properties and gravitational moments also becomes apparent in the plot. The deviations in gravitational moments decrease as the flow becomes shallower and as n becomes larger.




 §.§ Effect of Planet's Oblateness
 
We now address one of the key assumptions that we made in our model, i.e., the sphericity of the planet. Using our framework and perturbation theory, we can estimate the expected alteration in the calculated Δ J_n due to the non-spherical shape of the planet. To first approximation, the small deviation from sphericity can be captured with the second order Legendre polynomial:


    r → r (1 + ε P_2 ( cosθ) )


We Taylor-series expand the expression for gravitational moment J_n (equation <ref>), spherical Bessel function j_n (equation <ref>), and V_pn (equation <ref>) in the small quantity ε, which scales with q. The product of P_2 with P_n can be written as a sum of three Legendre Polynomials, P_n-2, P_n, and P_n+2, where each of the polynomials is weighted by the corresponding Clebsch-Gordon coefficient <cit.>. Expanding in small ε and keeping the first order terms (see Appendix C), we get:


    Δ J_n,oblate≈Δ J_n,sph  - ε[ 1/4 (n+1) Δ J_n,sph + 3/8 (n-1) Δ J_n-2,sph  +  3/8 (n+3) Δ J_n+2,sph]

where the subscript sph stands for spherical i.e. values of Δ J_n calculated under the spherical planet approximation. The oblateness of a planet couples gravitational moments of different order.

Importantly, the effect of oblateness is small and of the order ε∼ 0.1 for Saturn. Another thing to note is that Δ J_n-2 and Δ J_n+2 are opposite in sign to Δ J_n for the flow profiles considered here[This does not necessarily have to be the case. The sign of Δ J_n depends on the functional form of differential flow ϵ and the decay depth adopted and can change when the depth is changed.] and the expression in the bracket is small. The effect of oblateness is roughly proportional to n, which implies that it has a greater effect on higher order Δ J_n. This makes intuitive sense because higher order Δ J_n probe density in regions very close to the surface and therefore are more sensitive to the shape and oblateness of the planet.

Numerous studies have investigated the impact of assuming a spherical background state on Δ J_n's and compared their values from the (spherical) Thermal Wind equation and methods such as Concentric Maclaurin Spheroids (CMS) or Euler equation that account for the non-spherical shape of the planet <cit.>. The differences in the Δ J_n's values from these two approaches do not exhibit the simple proportionality to the gravity moment degree n that we propose above. We suspect this is because our perturbative approach breaks down for flows of moderate depth.



§ CHARACTERISATION OF SATURN'S DIFFERENTIAL ROTATION



 §.§ Numerical Calculations
 

Cassini's measurements have revealed the deviations of higher order gravitational moments from solid body rotation and we adopt rough values of these in our work <cit.>:


    ΔJ_6 ∼ (5 ± 0.5) × 10^-6  ,  ΔJ_8 ∼ (-5.5 ± 0.5) × 10^-6  ,  ΔJ_10∼ (3.5 ± 0.5) × 10^-6


What is remarkable about these values is their similarity in magnitude even though n increases from 6 to 10. It is immediately clear that an exponentially decaying flow cannot match these observations for any reasonable value of a, the amplitude of differential flow (Figure <ref>). We therefore must resort to another form for the flow profile.

The next natural step in formulating a flow profile is to consider the observed cloud-top motion as a surface manifestation of differential rotation in the planet's interior. There is no reason to suppose that winds and flows in the deep interior match the flow observed at the cloud-top level <cit.>. However, a sinc-like cylindrical flow profile that matches the surface flows of Saturn roughly reproduces the higher order gravitational moments surprisingly well <cit.> and also happens to have nice analytical properties (it is coincidentally the zeroth order Bessel function). In reality, it is not the lack of a differential flow profile that fits the data that has been the problem for Saturn. Degeneracy is the reason we cannot make conclusive statements about Saturn's flow properties using these measurements <cit.>.

Cloud-top motion extended deep into the planet as flows on cylinders <cit.> is plausible provided the cylinders are external to the electrically conducting region. It is a useful working hypothesis due to the simplicity of the centrifugal potential that arises as a result of this flow. A more realistic model would include a radial depth dependent decay term <cit.> but we cannot accommodate a radial dependence in our simple model because it violates Equation <ref>. We therefore adopt a general sinc-like flow profile in our study.

For Saturn, differential rotation, as manifested on the surface and in the zonal flows, is known to be the order of a few percent of the solid body rotation rate (∼ 4 %). This is a rough estimate and the exact value is not known. [This is due to the difficulty posed by the inference of the solid body rotation rate for Saturn <cit.>, which still has some uncertainty and cannot be inferred from the highly axisymmetric magnetic field <cit.>. This uncertainty is much smaller than the zonal flows however.] Given these uncertainties and the fact that flow in the interior may be different from cloud-top motion, we allow the amplitude and the depth of the flow to vary significantly from the values inferred from cloud-top motion.



Three broad functional forms based on the sinc function are chosen in the study (Figure <ref>). The first is a sinc function truncated at its first zero with no flow interior to it.


    ϵ (s) = 
    
    A  sinc( B (1-s ) ),     1 -  π / B < s < 1 
    
    0 ,     0 < s <  1 - π / B


The second is a sinc squared functional form which only has prograde flow and is truncated at the position of second zero.


    ϵ (s) = 
    
    A_1  sinc^2( B_1 (1-s ) ),      1 - π / B_1 < s < 1 
    
    A_2  sinc^2( B_2 (1-s ) ),      1 - π / B_1 - π / B_2 < s < 1 - π / B_1  
    
    0 ,     0 < s < 1 - π / B_1 - π / B_2


The third form has the equatorial prograde flow and mid-latitude retrograde flow and the sinc function is truncated at the position of the second zero.


    ϵ (s) = 
    
    A_p  sinc( B_p (1-s ) ),      1 - π / B_p < s < 1 
    
    A_r  sinc( B_r (1-s ) )     1 - π / B_p - π / B_r < s < 1-  π / B_p 
    
    0 ,     0 < s < 1 - π / B_p - π / B_r


For the second and third functional forms that have mid-latitude flow (prograde or retrograde), we allow the equatorial and mid-latitude flow to have independent amplitudes and depths. Although analytic forms are chosen for ϵ (s), subsequent calculations inadvertently involve numerical solutions and approximations. This is because analytical calculations for the sinc profile are time consuming. The entire calculation can be performed rapidly for exponential flow profiles. Therefore, we choose to represent the functional forms presented above as a linear combination of exponentials. A set of 15 exponentials (giving the desired accuracy) with a wide range of decay depths is chosen for this purpose:


    ϵ (s) = ∑_i=1^15α_i e^-2i (1-s)


The details of the numerical calculations are given in Appendix D. We now use this simple model to calculate gravitational moment deviations due to differential flow in Saturn and statistically characterize the dependence of these deviations on flow properties. This allows us to quantify the degeneracy in flow properties: a task that cannot be accomplished using the available suite of complex models due to their computational cost. Note that our approach does not work for Jupiter because it has significant hemispheric asymmetry, as indicated by the measurement of the odd moments <cit.>.





 §.§ The Necessity of Mid-Latitude Retrograde Flow


Efforts to fit the observed gravitational moments using highly accurate and sophisticated forward models have indicated that mid-latitude retrograde flow seems to be necessary to match the data <cit.>. We use our simple linear model to find out whether we reach the same conclusion or not. This serves as a test of both our model and the conclusion reached by other studies.



The deviations in gravitational moments due to prograde flows (equatorial and/or mid-latitude; defined in Equation <ref> and <ref>) are shown in Figure <ref> and <ref>. The different flow profiles are plotted in the left panels of the figures and the corresponding values of Δ J_n are shown in the right panels. The dashed lines in the plot mark the measured values of the gravitational moments. None of the flow profiles agree with the observations. The calculated Δ J_n for prograde flows exhibit behaviour that is similar to those obtained for exponentially decaying profiles. That is, the prograde flow are unable to produce Δ J_n of similar magnitude for n = 6, 8, and 10. We therefore conclude that for cylindrical differential flow profiles, the absence of a retrograde flow at intermediate latitudes makes it difficult to fit the observed values of Δ J_n. This is in agreement with other studies that use more sophisticated models.



 §.§ Non-Uniqueness of Flow Properties

Given that a prograde flow is not able to reproduce the observed Δ J_n values, we focus on our third general flow profile: a sinc function truncated at its second zero (Equation <ref>). This function is characterised by four parameters: the prograde flow depth and amplitude (D_p and A_p) and the retrograde flow depth and amplitude (D_r and A_r)[The depth and amplitude are normalised by Saturn's radius and solid body rotation amplitude respectively]. To understand the degeneracy in flow properties, we use affine-invariant Markov Chain Monte Carlo (MCMC) simulations to constrain the parameters that characterize Saturn's differential flow. This is made possible due to our choice to write all flow profiles as a linear combination of the same set of exponentials (Equation <ref>). The gravitational moments due to each of these exponentially decaying flows is calculated in advance and for each flow profile, we simply sum the weighted contribution of these moments (see Appendix D). The widely used Python package emcee is used to perform the MCMC simulations <cit.>. We assume uniform priors for the four model parameters that characterize the data flow and plot the corresponding posterior probability distributions from the MCMC analysis in order to determine the information provided by these new data.

As mentioned before, the chosen range of parameters is inspired by the observed surface zonal wind profile. However, we allow for variation in flow properties to account for the possibility that surface flows do not trace the flows in the interior. We use 20 `walkers' and 35,000 steps in each chain, which gives us a total of 700,000 samples of the posterior probability distribution. We initialize our walkers at randomly chosen positions clustered near a solution that is in reasonably good agreement with the data. We then run the MCMC for an initial 1000 step burn-in phase in order to ensure that all walkers have reached the preferred region of parameter space. The integrated auto-correlation time (τ_f) for each of the model parameters is ∼ 100 - 150. This implies that we draw 700,000/τ_f∼ 4500 independent samples from the posterior distribution, indicating that we have adequately sampled the relevant parameter space <cit.>.



The posterior probability distributions for the prograde and retrograde flow amplitude and depth are shown in Figure <ref>. This corner plot gives us the posteriors for each of these parameters as well as the correlations for all pairs of parameters. The surface (cloud-top) flow amplitude and the extent of the prograde jet on Saturn are marked by red lines in these plots. One thing that becomes immediately evident is that the data do not constrain the parameters very tightly. The posteriors for the parameters taper off at the edges of the parameter space but are not that different from the assumed flat prior for the MCMC simulations. This implies that a wide range of flow parameters can reasonably explain the observed gravitational moments, given the errors in our measurements.

Correlations between flow parameters allow us to understand the non-uniqueness of the possible solutions and place greater constraints than the 1D histograms. We find that the flow depth and amplitude for both the prograde and retrograde flow are strongly correlated. A shallow flow with a large amplitude and a deep flow with a small amplitude produce similar values for the deviations in gravitational moments. Moreover, we note that correlations between properties of the prograde and retrograde flow are weak. The amplitude of the prograde flow and the depth of the retrograde flow are almost uncorrelated. As for the other parameter pairs, we find that the depth and the amplitude of prograde and retrograde flow are roughly anti-correlated with each other (the actual posterior distribution is not simple). That is, the retrograde flow tends to be shallow if the prograde flow is deep (and vice versa) and the amplitude of the retrograde flow tends to be smaller if the prograde flow is strong (and vice versa).

We plot a sample of flow profiles, colored according to their log-likelihood (- χ^2), for which the MCMC calculations are performed in Figure <ref>a. These profiles are a part of the MCMC `chain' as the simulation explores the parameter space to fit the observed gravitational moments. The profiles with the highest log-likelihood (equivalent to the smallest χ^2) best fit the observations. Smaller amplitude flows that are deep tend to fit the observations as well as larger amplitude shallow flows.



Inclusion of retrograde flow clearly improves the agreement between the data and the model (Figure <ref>b). It is therefore worth noting how the properties of retrograde flow influence the gravitational moments. This will help us understand why including retrograde flow seems necessary to fit the data. Essentially, models with only prograde flow that match Δ J_6 and Δ J_8 underestimate the value of Δ J_10 (see Figure <ref> and <ref>). Including retrograde flow and increasing its amplitude or depth increases Δ J_10 while leaving the other two gravitational moments relatively unaltered (Figure <ref>c). This is what allows us to obtain a reasonable match between the model and the observations.



§ DISCUSSION AND CONCLUSIONS

The ideas and theory presented here aid in the development of an intuitive understanding of how the properties of the differential rotation are related to the deviations (from solid body rotation) in gravitational moments. Gravitational moments are given by the integral of density and a Legendre polynomial of certain order. In our formalism, density and gravitational potential are related simply: 2 K ρ = V + V_rot. The contribution due to differential rotation, δ V_rot, therefore directly relates to the density, and hence to the gravitational moments, which are now given by the integral of gravitational potential with the Legendre polynomial. Additionally, δ V_rot = q ∫_0^s 2 ϵ(s') s' ds' is related to cylindrical flow properties contained in ϵ. This establishes a somewhat more comprehensible relation between differential rotation properties and the gravitational moments.

Using this model, we have shown that deviations in gravitational moments alone do not suffice to constrain the differential flow properties of Saturn very accurately due to the errors in the measurement and the modeling. Even using a single sinc-like cylindrical flow profile, which already severely restricts the parameter space, we are able to match the observed gravitational moments for a wide range of flow properties. There is added non-uniqueness in the functional form of the flow profile as well and the surface flows, although coupled to the interior dynamics, may be susbtantially different from those in the interior <cit.>. To deal with this non-uniqueness, we must seek some other ways of constraining the flow depth or amplitude which will narrow down the range of possible flow properties. This requires invoking additional physics that might place constraints on the differential flow of Saturn.

One promising way of constraining flow depth is by studying the interaction of interior differential flow and planetary magnetic fields. Electric conductivity rises steadily as one ventures deep into the planet's interior and any differential flow would lead to the generation of currents and Ohmic dissipation. In the deep interior, the conductivity is very high and any differential rotation is quickly damped out by MHD drag. However, in the semi-conducting region, it is possible to have small amplitude flows with speeds of ∼ 1 cm/s to 1 m/s. Electrical conductivity and Ohmic dissipation along with the measured interior heat flux of the planet place this upper limit on the amplitude of differential flow because strong flows would generate a large amount of heat, which would lead to a discrepancy between the expected and the measured heat flux. <cit.> studied this interaction in the semi-conducting region of Jupiter and Saturn to determine if flow induced magnetic field variations can be measured with the Juno mission or the Cassini grand finale orbits. Their work indicates that for an intermediate electrical conductivity σ≈ 10^-2 S/m, it is not possible to have large amplitude (100 m/s) flows. However, 1 m/s flows are possible. The electrical conductivity reaches this value at ∼ 0.85 R_Saturn, which implies that the differential flow must decay from 100s of m/s at the surface to a negligible 1 m/s at this radius. This could be used as a constraint on the flow depth in our efforts to determine the differential flow properties. Indeed, the measured gravity moments of Saturn seem to imply an upper limit on the flow depth that is consistent with this constraint placed by MHD <cit.>. Notably, Ohmic dissipation places a constraint on the radial depth of the flow, not just the cylindrical depth r_c considered in this work. It is therefore possible to have differential rotation at high latitudes as long as it becomes negligible at this transition radius.

Using such a constraint and our results (in particular, Figure <ref> and <ref>), we note that deep flows with small amplitudes are excluded and shallow flows with large amplitudes are preferable because they agree with the gravity data as well as the theoretical expectations from magnetohydrodynamics and Ohmic dissipation. More specifically, the flow depth needs to meet the following conditions: 1/D_p + 1/D_r ⪅ 25 with both 1/D_p ≥ 6 and 1/D_r ≥ 6. The required amplitude of the flow would depend on the flow depths chosen and it can vary from ∼ 4 - 8%.

In conclusion, we have presented a linear model that relates differential flow to gravitational moments in this paper. We demonstrated its utility in the development of an intuitive understanding of the underlying phenomenon and the rapidity of the calculations that enable a statistical analysis of the flow parameters. Our analytical calculations for a shallow flow are useful for understanding the relationship between flow properties and resulting gravitational moments. In addition, we have derived an expression that quantifies the effect of a planet's oblateness on its gravitational moments, providing first order corrections to the assumption of sphericity that we made to make our model linear. This expression will prove useful in comparing the effect of oblateness with other small corrections to the gravitational moments. Moreover, the linearity of the model allows us to make a quick forward calculation for the gravitational moments given a (cylindrical) flow profile. We used this feature of our model to constrain the properties of Saturn's differential flow based on the higher order gravitational moments measured by Cassini by performing MCMC simulations. This calculation allowed us to quantify the widely suspected non-uniqueness of Saturn's flow parameters inferred from the gravity data. We found that retrograde flow seems necessary to explain the observed gravitational moments as flows that are purely prograde tend to underestimate Δ J_10. A wide range of flow parameters yield gravitational moments that agree with the observations and the flow depth and amplitude are anti-correlated. Given the non-uniqueness of the flow properties, additional physics needs to be invoked to place tighter constraints on them. Matching the gravitational moments along with the theoretical expectations from Ohmic dissipation due to flow in the semi-conducting region of Saturn and its heat flux might provide an additional way of diminishing the allowed parameter space of differential flow properties.



§ ACKNOWLEDGEMENTS

Y.C. thanks Hao Cao for stimulating discussions, Steve Markham for help with Mathematica ©, and Heather Knutson for discussions on MCMC. We are grateful to the reviewers for their thoughtful comments and suggestions which improved this manuscript. This work uses matplotlib <cit.>, scipy <cit.>, and sympy <cit.>.

apalike




§ SADDLE POINT METHOD

The saddle point method is used for estimation of integrals of functions with very sharp peaks and it was first published by <cit.>. For an integral of the form below for real functions f(x) and g(x) and for large A > 0


    I (A) = ∫_x1^x2 f(x) e^Ag(x)dx


Since the function is strongly peaked at some x = x_o for large A, we can approximate the integral as:


    I (A) ≈ f(x_o) e^Ag(x_o)√(2 π/-A g”(x_o))( 1 + ∑_n=1^∞C_2n/A^n)


The C coefficients depend on the derivatives of f (x) and g (x). We include C_2 in our model and it is given by: C_2 = -f”/4fg” + f' g”' / 8 f (g”)^2 + g”” / 32 (g”)^2 - 5 (g”')^2 / 192 (g”)^3.



§ SHALLOW FLOW: DETAILED CALCULATION

Let us assume that the differential flow profile can be represented by a steeply declining exponential which is only a function of the cylindrical radius of a fluid planet:


    ϵ (s) = a   e^-b   (1-s)


Here, a and b are constants and s is the cylindrical radius (as a fraction of the equatorial radius). We can integrate this function analytically to obtain δ V_rot:


    δ V_rot (s) = 2 q ∫_0^sϵ (s')   s'  ds' = a   q   e^-b  1 + e^bs (b s - 1)/b^2


Solving the equation for V_p (Equation 6) and using the condition that V_p and its derivative are zero at s = 0 to set the constants:


    V_p (s) = -  2 a    q   e^-b k^2/b^3[   b s^2/2  + e^b s( s - 3/b) + 2s + 3/b]


To proceed further with the analytic calculation, we need to project V_p onto spherical coordinates using the Legendre Polynomials (Equation 7 and 8). For very shallow flow, we can use a Taylor series expansion of the Legendre Polynomials up to a certain order (that we choose) and calculate the integrals analytically. We choose to expand the even Legendre Polynomials to second order in cos θ around θ = π/2:


    P_n (cos θ) = n  n/2(-1)^n/2/2^n[ 1 - n   (n+1)/2 (cosθ)^2 ]


With this approximation, we can project V_p onto a basis set of Legendre Polynomials with s = r  sinθ. We then have an expression for V_pn from Equation 8:


    V_pn (r) = - 2n + 1/2n  n/22 a    q   e^-b k^2/b^3(-1)^n/2/2^n ×
    ∫_-1^1[   b (r  sinθ)^2/2  + e^b r  sinθ( r  sinθ - 3/b)  + 3/b + 2 r  sinθ] ×
    [ 1 - n   (n+1)/2 (cosθ)^2 ] dcosθ


We calculate this integral in the approximation of shallow flow such the exponential approaches a delta function. In this limit, we can neglect all terms except the exponential in the expression for V_p and the expression for V_pn simplifies to:


    V_pn (r) ≈ - 2n + 1/2n  n/22 a    q   e^-b k^2/b^3(-1)^n/2/2^n ×
    ∫_-1^1 e^b r  sinθ( r  sinθ - 3/b)
    [ 1 - n   (n+1)/2 (cosθ)^2 ] dcosθ


To evaluate this integral, we resort to the saddle point method, which works well for strongly peaked functions such as the exponential in the expression above.


    V_pn (r) ≈ - 2n + 1/2n  n/22 a    q   e^-b k^2/b^3(-1)^n/2/2^n ( r - 3/b) √(2 π/b r) e^b r×
    [1  - 1/4br( 1/r-3/b + n(n+1)+3/8) + O  (1/b^2r^2) ]


Here, the first term in the parenthesis is the O(1) term which does not depend on b or r. The second term is the order O(1/b) expansion in the saddle point approximation. We calculate the expression for Δ J_n for both these orders of expansion below. Evaluating V_pn and its derivative at r = R = 1 (fractional radius in our calculations) gives us our final expression for Δ J_n:


    Δ J_n = - 2n + 1/2n  n/22 a   Ω_o^2  π^2/b^3(-1)^n/2/2^n√(2 π/b) [ -5/2 + b + 3/2b] j_n(π) - π[ 1 - 3/b] j'_n(π)/(n+1)j_n(π) + π j'_n(π)



    Δ J_n = - 2n + 1/2n  n/22 a    q  π^2/b^3(-1)^n/2/2^n√(2 π/b) ×
    [ ( -5/2 + b + 3/2b) + n (n+1)/4b( 7/2 -b - 9/2b) ] j_n(π) - π[ (1 - 3/b) ( 1 - n(n+1)/4b) ] j'_n(π)/(n+1)j_n(π) + π j'_n(π)




§ EFFECT OF PLANET'S OBLATENESS

We discussed the consequence of relaxing our assumption about the spherical shape of the planet on the gravitational moments of a celestial body. Here, we show the calculation in detail instead of simply quoting the results. The result we present is essentially derived from the continuity of the gravitational potential at the surface of the planet, which is no longer spherical but is as described by:


    r → r (1 + ε P_2 ( cosθ) )


Here, ε is small (typically << 1)and P_2 is the second order Legendre polynomial. The boundary condition for the potential at the surface is:


    ∑_Even   n V_pn (r) P_n (cos θ) + ∑_Even   n A_n j_n (kr) P_n (cos θ) = - ∑_Even   nΔ J_n ( R/r)^n+1 P_n (cos θ)


where we have equated the sum of the particular and the homogeneous solutions for the potential interior to the body to the potential exterior to the body. For a spherical planet, we equate terms corresponding to the n^th Legendre polynomial and set on r = R, which gives us:


    V_pn (R) + A_n j_n(kR) = - Δ J_n


However, for an oblate surface, r = R no longer represents the boundary and is replaced by r = R ( 1 + ε P_2 ( cos θ) ), which introduces a θ dependence in the expression. Since ε is small, we can do a Taylor series expansion of V_pn and j_kr around r = R to obtain an approximate expression for these quantities when the planet is oblate.


    V_pn≈ V_pn (R) + d V_pn/d rε P_2 (cos θ)  R 
     
    j_n ≈ j_n (kR) + d j_n (kr)/d rε P_2 (cos θ) R


Substituting these expressions in the boundary condition gives terms with products of two Legendre polynomials P_2 and P_n. The product of two Legendre polynomials can be written as a sum of Legendre polynomials weighted by the Clebsch-Gordon coefficients <cit.>.


    P_p (cos θ)   P_q (cos θ) = ∑_i = 0^q B_i B_p-i B_q-i/ B_p+ q -i 2p + 2q - 4i + 1/2p + 2q -2i +1  P_p+q-2i (cos θ)  for  p ≥ q


where the B coefficients are given by:


    B_m = ( 1/2 ) ( 1/2 +1 ) ... ( 1/2 + m - 1 ) /m ! with  B_0 = 1


Given that we equate coefficients of a particular n^th order Legendre polynomial at the boundary, including the effect of oblateness via P_2 couples the n^th gravitational moment with the order n-2 and the n+2 moments by including terms that emerge from the multiplication of P_2 with P_n. The equality of interior and exterior potential for an oblate planet therefore gives us:


    V_pn (R)  + j_n (kR)
    + ε[          1/4( R V_pn'(R) + π A_n j_n'(kR)  ) +; 3/8( R V_p(n-2)'(R) + π A_n-2  j_n-2'(kR) )  +;     3/8( R V_p(n+2)'(R) + π A_n+2 j_n+2'(kR) ) ]
    =
    - Δ J_n,oblate


Here, the Clebsch-Gordon coefficients are exact in the limit of large n but serve as a useful approximation for the range of n we are interested in. The expression in the parenthesis is related to gravitational moments for a spherical planet by the continuity of the gradient of the potential and is simply given by:


    R V_pn'(R) + π  A_n  j_n'(kR) = (n+1) Δ J_n


Substituting these relations into our expression for Δ J_n,oblate, we obtain the final form of our result:


    Δ J_n,oblate = Δ J_n,sph  - ε[ 1/4 (n+1) Δ J_n,sph + 3/8 (n-1) Δ J_n-2,sph  +  3/8 (n+3) Δ J_n+2,sph]




§ NUMERICAL CALCULATIONS

All the numerical calculations are carried out both in Python and Mathematica. Here, we describe the details of the procedure that we adopt for numerical calculations of Δ J_n. Instead of calculating Δ J_n directly for sinc-like differential flows profiles, we write them as a sum of 15 exponential profiles (Equation <ref>).
Numerical solution for exponentially decaying flow profiles is rapid. In addition, since we use the same 15 exponentials for all profiles and its only their coefficients that vary, we just need to perform this calculation once for each of these exponential functions. Thereafter, Δ J_n for any profile can simply be written as a sum of Δ J_n,i obtained from the i^th exponential profile and weighted by the respective coefficient due to the theory's linear nature:


    Δ J_n = ∑_i=1^15α_i Δ J_n,i


For the exponential profiles, we then calculate Δ J_n as follows:

    
  * We integrate Equation <ref> to find δ V_rot using Sympy integrate (Python) or simply `Integrate' in Mathematica ©.
    
  * Given δ V_rot, we can solve Equation <ref> to find V_p using any of the differential equation solving routines (Sympy `dsolve' or Mathematica©`DSolve'). The calculation is analytical thus far and has already been described in detail in Appendix B.
    
  * This is step where we deviate from our analytical solutions. Instead of using an approximation for the Legendre polynomials, we numerically integrate Equation <ref> to find the projection of V_p onto spherical coordinates. The numerical integration is performed using `NIntegrate' in Mathematica©or Numpy `trapz' in Python (the trapezium method only introduces small errors but is much faster and allows us to perform MCMC). In addition, Equation <ref> tells us that calculating dV_pn / dr simply involves differentiating V_p with respect to r inside the integral (which is an analytical step), followed by numerical integration of the expression. We can then use V_pn and V'_pn to calculate Δ J_n from Equation <ref>.
    


