Modeling Jupiter with a Multi-layer Spheroidal Liquid Mass Rotating Differentially

With the aim of improving the Jupiter equilibrium liquid model consisting of two distorted spheroids (“spheroidals”) of our last paper, we generalize it here to any number l of layers, demanding that the calculated gravitational moments, J2n (n=1,.., 4), agree with those surveyed by the Juno mission, which is fulfilled with a much higher accuracy than for l=2. The layers are of constant density and concentric (but otherwise free from any specific constriction between their semiaxes), each rotating with its own distribution of differential angular velocity, in accordance with our law in a past work. We point out that the angular velocity profiles are a consequence of the equilibrium itself, rather than being imposed ad initio. Although planetary structure aspects are not contemplated in our models, we arrange matters so that they can be compared with Gudkova’s and Guillot’s, paying attention on the distributions of mass and pressure. Our procedure is exact, in contrast with the self-consistent CMS (Concentric Maclaurin Spheroids) method developed by Hubbard, whose inexactitude resides in maintaining a single constant angular velocity while the spheroids are deformed. Our model predicts a differential rotation for Jupiter with periods for pole and equator of 9h38m and 10h14m corresponding to a mean period of 9h55m.


Introduction
In a past paper [1], a Jupiter model consisting of two liquid concentric distorted spheroids ("spheroidals"), core and envelope, rotating differentially was proposed, demanding agreement between the calculated gravity moments J2n and those surveyed by the Juno mission, and that the model was in hydrostatic equilibrium. For the envelope, the mean rotation period found was 9.92 h, the average being taken over the values at the equator and the pole. The core resulted highly distorted, too small, and rotating very fast. The calculated J 2n fell slightly off center of their uncertainty regimes. This model resulted to be stable. In an attempt to improve it, in the current paper the body is generalized to l layers. The minimization procedure used for the calculation of the J 2n is achieved with a mean accuracy of ∼ 10 −10 , falling much nearer to the uncertainty centers than for l=2, where the accuracy was ∼ 10 −1 . Next, the corresponding equilibrium figures are established.
The theory of stellar structure has evolved considerably. Although the interior of a star cannot be seen, there is little doubt of what is going on there. In the case of giant planets, the tendency is to imitate the stellar theory as much as possible. In recent times, there have been some important advances, but the uncertainty in the picture of the interior of planets still persists. Important problems occur principally, among others, in the not yet fully understood over-presence of heavy elements, the lack of exact knowledge of equations of state (EOS) for H/He or heavy elements at high pressures, as well as phase changes.
The Juno mission allows the determination with good precision of Jupiter gravity, opening the opportunity to check the plausibility of the several interior structure models proposed thus far. Efforts to explain the observed gravity are made, for instance, by [2], on the abundance of heavy elements, or by [3], which compare various EOS. An alternative to EOS is related to a technique employing simulations known as Density Functional Molecular Dynamics [4] which, at a quantum level, avoids the approximations and uncertainties of EOS.
The Jupiter rotation state is also an open field. For example, [5] suggest the possibility that its atmosphere rotates differentially, but that internally rotates as a solidbody. Some of the cited references, and others like [6][7][8][9][10][11][12][13] use the CMS (Concentric Maclaurin Spheroids) method, developed by [14,15] to determine the J 2n . This method is not exact because the expression used for equilibrium is not valid (see sections 4 and 7). The interior structure models are not affected by the imprecision if they are not forced to reproduce the J 2n . However, if they are, as in the case of simulation procedures like those of [5], in which models are picked up from all the simulated ones that reproduce the observed J 2n , as inferred with the CMS procedure, this structure models might actually not have the correct mass distributions for reproducing the observed J 2n .
Our method to establish the J 2n is based on an exact equilibrium equation, described in detail in [16], and applied in [1]. This time, an application of the l-layer generalization based on existing internal structure models is made. For this purpose, we have chosen the models of [17] and [18], although more recent ones are available. Our main interest is to describe how to use the model, without intending to establish the suitableness of a particular one. Moreover, their models present quantitative results that are adequate for our method. One of our results concerns with a differential rotation for the envelope, with periods of pole and equator of T p =9h38m and T e =10h14m. We point out that the angular velocity is a consequence of the theory, and hence it does not require of any constraint. Furthermore, the layers are deformed spheroids ('spheroidals'), with adjustable size and deformation by the numerical procedure.

The Figures' Shapez
As in previous works [16,19], the layers' surfaces are analytically proposed from the beginning and given in the normalized form by the equations

The Gravitational Moments
The gravitational potential of a mass distribution of density ( ) r ρ (cylindrical symmetry assumed) is established as a series of gravitational moments given by a is the body's equatorial radius, τ the volume, P 2n the Legendre polynomial of order 2n and ϑ the azimuthal angle measured from the pole to the equator. Let τ i be the total volume limited by the th i layer's outer surface, and l the number of layers; notice that τ 1 =τ. Since ρ (r)=ρ i for points r of th i layer only, equation (3) can be written as and cylindrical coordinates (R, ϕ, z) have been used for convenience. The integrals in equation (5) are known analytically and given in the Appendix of [1] ( [1], n ≤ 4).

The Equilibrium State
In papers [16,1], the equilibrium equation for a rotating fluid was stated, specializing it for a two-layer model in [1]. The generalization to l layers is immediate. In the steady state, the generalization of Bernoulli equation is valid for each point of the medium: p is the pressure, V the gravitational potential and f an arbitrary function that is determined by the boundary conditions in the layers' interfaces: The total potential V i (R, z) in each interface point is a function of R only since z can be expressed in terms of R by equation (1). Functions f i are related to the angular velocity Ω (=ω 2 /Gρ 1 ) through ( [16]).
Hence, according to equation (8) Ωi for layer i is given by Here, r is an abbreviation for R2. Finally, from equation (7) the pressure in each point of a layer is For numerical procedures, the gravitational potential is normalized to Gρ1a2. In the center, the pressure is

The Numerical Procedure
The numerical method employed in our approach was already described in [1], so it will succinctly be reproduced here for a rapid reference. The method is carried out in two steps: the construction of the mass distribution that reproduces the observed J 2n , and the establishing of the angular velocity profile that sustains equilibrium. It must be emphasized that the last is not a restriction that needs to be imposed to the model; rather it obeys a law that we deduced and explained in detail in [16]. Hence, in our model the angular velocity is of a predictive character.
For the first step, the inputs are the following facts for Jupiter where a, b, M J are the equatorial and polar radii and the mass of Jupiter. Higher moments were not considered because the error bars increase considerably ( [20,21]), besides they are for the procedure irrelevant. With this input, the method produces as output a mass distribution having the observed parameters. It is characterized by the following variables: equatorial e i and polar z Mi radii, distortion parameters d ni and densities ρ i , of all layers. This comes from an optimization procedure as is explained next.
Since the theoretical moments that we call , are known analytically [1], it is demanded that variables e i , z Mi , d i , ρ i on which depend be such that This problem is solved numerically by minimizing the sum in equation (14) from which the model parameters result.
Knowing the model parameters, we come to the second step. The unknown here is the angular velocity that must satisfy a law depending on the gravitational potential (equations (8 to 10) coming from the equilibrium conditions), which can be determined for the already known mass distribution. For a specific application of the method, see sections 5 and 6.

Bizyaev et al. Model
Although this topic is not related specifically with Jupiter, it deals on rotating equilibrium figures of ideal fluids' density stratifications worked out in an exact fashion [22], without generally supposing that the angular velocity is constant. In particular, they consider an application of their method for a stratification consisting of confocal spheroids (all shells have common focal points). The results are already known ( [23], [24]) and are presented and applied here because their method is more general and can also be used for continuous stratifications.
Constraining the l-layer model to be composed of confocal spheroids, the focal length c is fixed already by the Jupiter polar and equatorial radii: where c n is the normalized focal length normalized Jupiter polar radius). It is important to remark that since all layers are limited by confocal surfaces, the gravitational potential at each surface point is a linear function of r (=R 2 =x 2 + y 2 ), and the angular velocity must be constant within each layer (equation (10), dV/dr=const). It is the exceptional case of a fluid in which the layers can rotate as solid bodies. For other instances, constant layer's angular velocity is not possible or it must be differential. Given the stratification, and hence the potentials, Bizyaev et al. come to the rotating state that sustains equilibrium: meaning z Mni the normalized polar radius, ∆ i =ρ i − ρ i− e and ∆ 0 =ρ 0 =atmosphere density. Knowing the layer's equatorial radii e i , by the confocality condition, the polar ones are fixed: and with them all angular velocities (16) (in terms of the densities ρ k ). Certainly, not all layers' geometrical configurations interest here, but only those reproducing the ascertained gravitational moments in the Juno mission. Hence, the attempt is made to find a configuration that reproduces the observed gravitational moments with the procedure described in section 3. The result should be to get the correct layer's radii e i (or z Mni ) and densities ρ i , and with them the necessary angular velocities (16) that maintain equilibrium. Many efforts to find an expected stratification was made, without success. We believe that this is due to the condition that a confocal stratification must exist that is able to reproduce the observed moments. However, experimenting with J 2n -values (especially J 2 ) that differ too much from the observed ones, a configuration can be found.

Gudkova Models
For the Gudkova models [17], for example model MJ8, the relevant information here is presented in Table 1, which we call model M1. R is the normalized distance from Jupiter center, ρ and P are density and pressure at the layers' interfaces. In this model, the density increases rather moderately in the first two layers (1.0-0.9 and 0.9-0.81), however, becomes considerably larger in the central layers, reaching the value 13.020 g cm −3 at the center that is about 250 times the value in the atmosphere. Regarding the pressure, its behavior is similar to the density's. The low initial growth comes to the two Mbar level at a depth of ∆R=0.2 and reaches so high values as 50 Mbar at the center. We now construct a 5-layer spheroidal model taking the above radii as the equatorial ones since our model is not spherical. We demand, firstly that the mass distribution be such that the gravitational moments J 2 , J 4 , J 6 and J 8 agree with those provided by the Juno mission. Putting aside the restriction of a layer fixed radius, we find several models with the exact J 2n , that is, the solution is not unique. Certainly, with the knowledge of only some gravitational moments, it is not possible, in general, to find a single mass distribution alone that reproduces J 2 ,..., J 8 and gives the potential (infinite series) solely. However, we do not accept any of these models, besides those that are in equilibrium under the action of gravity, pressure and rotation. In most cases, equilibrium is not satisfied because the centripetal force necessary for rotation cannot exist within a region. Although the error bars are more or less narrow, the observed J 2n are not exactly known. For this reason, we select randomly the moments within the error bars and seek the best mass distribution that reproduces them. We find several models that are somewhat alike. Next, their equilibrium is tested, i. e., we attempt to find an angular velocity profile sustaining equilibrium. A typical equilibrium model is M2 given in table 2. For this model one gets as gravity moments: 6 6 14696. 5   The model with parameters of M2 has low distortion (small d) in the three outer layers and, hence, they differ little from the spheroid (d=0); additionally, their flattening is low. The two innermost layers, on the contrary, are considerably more distorted and flattened. The density grows toward the center up to about ten times: from 0.209 to 2.030 g cm −3 . The pressure has a remarkable behavior: from the atmosphere inwards, it increases rapidly reaching the two Mbar level in the bottom of the first layer and 5 Mbar in the second one; it continues rising at a lower rate, reaches a maximum of 19 Mbar and then decreases to 17 Mbar at the center. The pressure here has no restriction, like in a barotropic relation or an equation of state, it only must satisfy expression (7). The angular velocity of a layer is influenced by that of the layer above it (equation (10)). As R decreases, the angular velocity increases, causing a greater increase near the center, resulting in a greater centripetal force, and hence in a reduction of the pressure gradient. The angular velocities are almost constant, not so the atmosphere (outermost layer) which clearly rotates differentially, with the pole running faster than the equator. The mean period of the atmosphere is 9.93 h, pretty near the accepted one [25]. Therefore, our model reproduces poorly Gudkova's general parameters (or conversely).
One source for the discrepancies probably resides in the layer concept. Our model's layer is a bounded region of constant density delimited by non-spherical surfaces. For Gudkova's (and some other researchers in the planetary structure field) a layer is a finite region also characterized by certain chemical and physical properties, such as composition, equation of state, entropy and so on. Their layers are commonly limited by spherical surfaces, and the density changes from a point to point. In order that our model comes closer to Gudkova's, and from a gravitational viewpoint, the structure model will accordingly be Liquid Mass Rotating Differentially considered composed by three (or four) shells: an external shell, including the thin atmosphere, a middle layer and a central core. In the outer part the density increases almost linearly to the center; in the one or two shells core, the density is approximately constant. Gudkova's models [17] have densities that increase nearly linearly with the (normalized) radius R , from the surface to the core ( 1 . 0 ≈ R for a one-layer core or 15 . 0 ≈ R for a two-layer one); the core has practically constant layer density. Consequently, a spheroidal multilayer body is built here, say of 15 shells, divided into three regions: an outer part consisting of a constant density thin atmosphere; a twelve or thirteen layered region with a nearly linear density profile, and a one or two constant density layer core.
A set of four gravitational moments , , , is randomly generated within the confidence intervals' limits of the observed J 2n . As explained in [1], a best mass distribution is looked for that reproduced the observed moments, that is, for which (see also section 3) where is the theoretical moment calculated with the functions given in the Appendix of [1]. It depends of polar and equatorial radii, distortion parameter and relative density difference of neighboring layers. The procedure's succeeding does not guarantee an equilibrium state for the particular matter distribution, i. e. to find a differential angular velocity distribution that sustains gravity and pressure. Regularly, the configurations with proper J 2n are not equilibrium figures. To obtain an equilibrium model with good J 2n , a huge number of gravitational moments' sets is randomly generated and from each one, the correct mass distribution is obtained. From these procedures we found a 15-layer equilibrium model characterized by model M3 of In this model, the layers' surfaces are not severely distorted, although the innermost deviates more from the spheroid's shape than the peripheral one. The density attains higher values than in the model M2, beginning at 0.344 and ending at 5.781 g cm −3 , with a similar behavior as the pressure that reaches about 30 Mbar in the center. Regarding rotation that supports equilibrium of the distorted model, the central part rotates about four times faster than the outer one; moreover, the angular velocity profile is markedly more differential in the first. The model's mean rotation period is 9.63 min that practically agrees with the observed one. Certainly, models with higher inner densities were found; however, they did not correspond to Gudkova's layers or were non-equilibrium models. Figures with a nearly linear density profile in the outer part and same size as Gudkova's layers, seem to have a central density limit of about 6.000 g cm −3 , a value close to a half of model M1. We presume that models with higher central densities cannot be found satisfying the particular conditions.
In our model, density increases inward more rapidly first than in Gudkova's, since, for instance, the layer e 1 =0.89 to 0.72, with a midpoint ≈ 0.81, has density 1.311, whereas for R=0.81 Gudkova gets 0.920. Near the core, the relative difference between the two models decreases visibly. It is in the nucleus where Gudkova's model density rises sensibly, deviating much more from ours. The not too severely different behavior of the two density profiles in the most of the models, could be attributed to a lack of distortion (and rotation) in Gudkova's model. Though, the excessive density jump in the core that may be influenced by the remaining layers, is due perhaps to some weakening in the structure considerations; model MJ6 comes closer to model M3 in the sense that the central density is slightly lower and core mass is small, about 3 M⊕ (in our model is about 1.5 M⊕). Our rotating distorted model reproduces well the observed gravitational field, but does not take yet into account structure equations; in particular, temperature (that surely affects density and pressure, low in our model) is absent.

Guillot Models
Regarding Guillot models [18,5], the relevant data are not given in the form of model M1, but as a relation between density ρ and pressure P that we write in Table 4: These results are to be compared with M3, expressed as a relation ρ -P, where ρ is a mean density ((ρ i−s +ρ i )/2) taken in each interface of our model and P is the pressure at the pole of the interface. It leads to model M5 of Table 5. After a closer inspection, it is found that pressure considered as a function of density is always below for model M4 as for M5. From the surface to ρ ≈ 4gcm −3 it rises much more slowly than in M5, afterwards P increases sensibly in the last small interval, reaching the values of M5. In other words, there is a gap between the pressures of both models that widens first up to ρ ≈ 4gcm −3 and afterwards diminishes more and more. If we use a mean polytrope P=Kρ n to describe roughly the models, then (M4) has a higher mean index n=3.2 than M5, n=1.3, but constant K of the former is smaller than for the last. However, we must emphasize that M5 was built supposing a nearly linear density profile in the outer part of the model, and we have insufficient data to know if M4 fulfill this condition. Nonetheless, recent similar models [26] shows a linear tendency of ρ against R. Although the ρ-P gap between the two models is not narrow, the order of magnitude of the variables is similar, and practically the same near the surface and the core. Discrepancies in the rest of the models could be attributed to a lack of physical and chemical conditions in our model; or to a need of more refinement in the chemical restrictions and (differential) rotation and deformation in the model M4.
Perhaps a more accurate model will be between M4 and M5.

Discussion
The models used here for a comparison date from 1999, although more recent ones are available [5]. They were chosen because they present numerical (graphical) results which are necessary for our model. From the papers, indistinctly if they are recent or not, one cannot get exact quantities, especially if they are given in graphical form. However, the principal interest here was to show the way to the application of our model, without the purpose yet to favor one particular internal structure model. Later, we will integrate into our model equations of internal structure.
With the data extracted from the models of Gudkova and Guillot, we have seen that our model predicts too lower pressures and densities near the center of Jupiter than that Gudkova's, but nearly of the same magnitude as that of Guillot's. Conversely, in the remaining of the model, the coincidences are better with the first than with second. To get a global picture of this, let us represent the models roughly by average polytropes (logarithmic scale): 6  it is seen that Gudkova's and our model have practically the same mean tendency and pressures, P Ga and P Ou , are slightly separated from one another; nonetheless, we know that there is a high discrepancy between the two in the central region.
On the other hand, Guillot's model pressure grows too much more rapidly with density than in the other cases, although it begins at a lot lower values; however, the order of magnitude of P Gt is alike P Ou . Why is the increasing rate of P Gt stronger than of P Ga ? Surely, it resides in the assumptions made in both models; but here cannot be decided what differences in the models induce, for example, a rapid pressure growth, or a very high central density or pressure. Last behavior is not probably due to a simple or multilayered core since Gudkova [17] gets similar ρ in both cases. On the other hand, with the data of the two models and that observed for Jupiter, our model predicts a Jupiter differential rotation of the atmosphere, that is given by the points where ω is expressed in rad h −s . Our models lack presently of physical and chemical properties of the fluid; the only supposition being that the layers are of constant density ρ i . Moreover, the stratification is deformed and rotates differentially to sustain equilibrium.
For an arbitrary homogeneous fluid with cylindrical symmetry, the angular velocity in the surface is necessarily given by 2 2 dV dR ω = − (19) where V is the gravitational potential at a surface point at distance R from the rotation axis. In Maclaurin spheroid case, the slightest deformation [14] causes that the potential will differ from the simple quadratic form in (x, y, z) (see Section 4) and the angular velocity will no longer be constant. Nonetheless, Hubbard does not allow to change ω. Certainly, for a multilayered object, ω in each interface point is expressed by an equation like (10). The model M2 or M3 based in Gudkova's and Guillot's models foresees a Jupiter surface rotating differentially with periods in pole and equator of T p =9h38m and T e =10h14m, with a mean period T m =9h55m. The method used here to predict the angular velocities in Jupiter surface is similar to that of [25,27] in the sense that they use the Jupiter gravitational field only to deduce a rotation period T=9h55m; they assume from the beginning solid body rotation. They, like us, do not put any constraints on the internal structure, hence the method is not structure model specific. Here was demanded only that density changes approximately linearly with the equator radius (inspired by Gudkova's and Guillot's models).
In the recent time, it is common to calculate the J 2n moments with the CMS procedure developed by [15]. It is not an exact method because his basic expression does not correspond to the equilibrium state: U=V + Q=const., where 2 2 (1 (cos ) / 3 2 Q r P ω θ = − . According to [16], a rotating axial-symmetric fluid is in equilibrium generally if it rotates differentially, independently of its shape. For a rotating stratification, each layer rotates differentially with its own angular velocity. As a special case, when the layers are confocal spheroids, their angular velocities are constant (Section 4). On the other hand, differential rotation cannot be common to all shells. For instance, let us suppose only two spheroid layers rotating with the same angular velocity. Equation (10) reduces to (Ω 2 =Ω 1 , Ω 1 is the common angular velocity) 2 2 1 dV dr Ω = − (20) Indices 1 and 2 refer to external and internal surfaces. By means of equations (8) and (9), we conclude that 2 1 dV dV dr dr = (21) It is to be noticed that both sides of equation (21) are evaluated at the same distance R (R 2 =r) from the rotation axis, but at different z (z 1 > z 2 ): point 1 is farther than point 2 from the center, hence the derivative at point 1 is smaller than at point 2, and the equality (21) is not valid.
According to our procedure, whose theoretical basis is presented in [16], there is no freedom in the choice of ω: it is determined by the gravitational potential at each point of the fluid (equations (9), (10) and (20)), and in general it is of differential type.