Proper Frequencies of a Rotating Two-layer Spheroidal Liquid Mass: A Jupiter Model

We establish the first six proper frequencies for the series of figures of our heterogeneous models. For this purpose, we use the second-order virial equations as developed by Chandrasekhar. The main aim is to know if the figures are stable. We consider all types of figures: those with an oblate or prolate nucleus, with ρn > ρa, as well as the cases with ρn < ρa. We attempt here also to fit our two-layer model to the planet Jupiter, using its known gravitation obtained by the Juno mission, along with other accepted facts. It is found that the model reproduces reasonably well the observed gravitational moments. Further, it comes out that Jupiter rotates differentially, with a mean period of 9.92 hr, agreeing with the recognized one.


Introduction
In our paper Cisneros et al. (2015), we obtained the first six proper oscillation frequencies of equilibrium figures from a self-gravitating homogeneous spheroidal (i.e., a distorted spheroid) liquid mass rotating differentially. In our first attempt, we proposed empirical expressions for the angular velocity profile; however, the frequencies were not altered by much once they were recalculated with the true rotation law, as we deduced (see the first of Equation (8)).
The figures were organized in series, that is, as a set of figures with a common distortion parameter value d and variable polar radius (flattening), in much the same way as the MacLaurin series (d=0). As in the former, our series have regions of instability (complex frequencies from a certain flattening), whereas some others presented neutral frequencies (= 0), thus leaving open the possibility of new triaxial figures branching off from the series.
In the current paper, we examine the stability of a mass composed of two concentric spheroidals ("nucleus" and "atmosphere"), the equilibrium figures of which were established with the correct rotation law. The series obtained are more involved than the homogeneous ones, because they depend on six parameters.
As an application of our previous studies and those that follow, we built a model for Jupiter based on recent observations made by the Juno mission. Specifically, we take the precisely known gravitational moments and look for the best two-layer geometry that reproduces them. It turns out that the model has good performance in this sense, and a differential angular velocity that sustains the model's equilibrium can be found. The core is relatively small, flattened at the poles, and deviates greatly from the spheroid's shape. Moreover, it rotates rapidly. The envelope also rotates differentially, but much slower, and differs little from the spheroid. Additionally, it results in the model being stable. Our model belongs to the class worked out by Hubbard (Hubbard 2012(Hubbard , 2013) and Debras (Debras & Chabrier 2017) in the sense that they are multilayer ideal fluids in equilibrium state sustained by rotation; the main goal with these models is to reproduce (planets') gravitational fields, coming more or less close to them. The fluid has no other properties such as density. They are not constructed from planet evolution considerations. Other models having an evolutionary basis like those of Wahl et al. (2017) and Guillot et al. (2018) have a determinate structure resulting from taking into account the state equations, the chemical composition of different regions (layers), and thermodynamical considerations; thus, the models are intended to give insight into the actual structure of the planet. As an application, some of their resulting figures reproduce acceptably well the gravitation deduced from the Juno mission. Our model (and that of Hubbard's) is not yet proposed to be applicable to a specific case, but only to serve as a basis for dealing with realistic cases.
The basic mathematical formalism needed is presented briefly in Sections 2 and 3 with minor adaptations, because it was already done extensively in previous works (Cisneros et al. 2015(Cisneros et al. , 2016. The main purpose is to have a rapid reference for the numerical and analytical procedures required here.

The Equilibrium Equations
As in previous works (Cisneros et al. 2016(Cisneros et al. , 2017, the nucleus and atmosphere surfaces are known analytically from the beginning and given in normalized form by the equations Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
The net potential at any core or envelope point is, respectively, where V nn is the nucleus self-potential, V an is the atmosphere potential on the nucleus, V na is the core potential on the envelope, V aa is the envelope self-potential, and ε is the relative density difference of the nucleus and atmosphere: In steady state (equilibrium state), the Bernoulli equations are satisfied: where f n and f a are functions to be determined and are related to the angular velocities through w = ¢ f 2 a a 2 and w = ¢ f 2 n n 2 . They are established from the surface conditions p n =p a (core surface) and p a =0 (envelope surface), p n and p a being pressures on the nucleus and atmosphere points: We have put r=x 2 +y 2 , and V N and V A are potentials on the core and envelope surfaces, taken as functions of r only. Our composite models depend on the parameters ε, e 1 , e n , e a , d n , and d a , which change freely, resulting in a wide variety of series, some of which have as special members those with a prolate nucleus. There are also particular series with a contact model. This was defined for z Mn =z Ma , i.e., both poles coincide, and hence e a and e n are related via

Virial Technique to Second Harmonic
The interest here is on the proper oscillations of the figure and for this reason, we propose for the small perturbation x at a point x and time t . 1 0 t If the eigenvalue λ is positive, the equilibrium state will be unstable; if it is negative, the motion will be damped; and if λ is imaginary, we have a stable oscillatory motion. To make the analysis easier, we take the approximation (Chandrasekhar & Lebovitz 1962b;Tassoul & Ostriker 1968 where L ij is a set of nine constants to be determined. The nine second-order virial equations are thus (Cisneros et al. 2015) ò 12 =−1, ò 21 =1, and ò ij =0 for the remaining cases. The supermatrix W lk;ij is defined as (Chandrasekhar & Lebovitz 1962a) ò ò r r Since ρ(x)=ρ n if xòτ n and ρ(x)=ρ a if xòτ a , the double volume integrals (13) and (14) can be expressed in terms of τ n and τ (= τ n +τ a ) as ò ò ò ò ò ò ò ò where ρ a =1 after normalizing. According to our figure's symmetry, only 7 of the 81 supermatrix components W pq;ij are unequal and different from zero: Moreover, all these components can be computed in terms of W 33;33 (for example), W 11 , and W 33 (Chandrasekhar & Lebovitz 1962a): For ease of writing, we abbreviate This is a homogeneous equation system in L ij , the determinant of which must be zero in order to have a nontrivial solution. Instead of working with the whole determinant, we solve subsystems, as is commonly done. This is possible by first noticing that the unknowns L 13 , L 23 , L 31 , and L 32 -and only these-are present in four of the equations, which bring us to the initial oscillation type.

Transverse-shear Modes
These modes are related to the displacements according to  the determinant of which must be zero: 1 J 1 and J 3 are the inertia moments relative to the x 1 and x 3 axes (see Equation (34)). Equation (20) is the eigenfrequency dispersion relation.

Toroidal Modes
The remaining five equations are l l    and these can be used to obtain the toroidal modes. Subtracting Equation (22) from (23) and adding Equations (25) and (26), we obtain a system in the two variables (L 11 −L 22 ) and (L 12 +L 21 ), the solution of which is subject to the vanishing of the determinant (λ=iσ and 2W 12;12 =W 11;11 −W 11;22 ):

Pulsatory Mode
We come to the last mode following a typical procedure for incompressible fluids (Chandrasekhar 1969, p. 84). Adding Equations (22) and (23), and subtracting twice Equation (24), one gets l l l Subtracting now Equation (25) from (26) Finally, we eliminate L 12 −L 21 and L 11 +L 22 from Equation (29) with the help of Equation (30) and Chandrasekhar's Theorem 17 (Chandrasekhar 1969;Cisneros et al. 2015) based on the figure's shapes and the continuity equation

Goldreich's Criterion for Local Stability
Goldreich's criterion for the local stability of gaseous stars is (Goldreich & Schubert 1967) where s 2 =x 2 +y 2 ; ω(~W) is the angular velocity, with k s and k z being the wave numbers relative to the radial and z dimensions, respectively, and where it is generally assumed that ∂ω/∂z>0 (in the present work, ∂ω/∂z=0), so that as usual ω increases from the equator to the pole; this might imply that for sufficiently large k s /k z , the value of the second term in the above inequality can override the first, making its left-hand side negative. Hence, to obtain local stability, the angular velocity distribution should be such that ∂ω/∂z=0, and Goldreich's criterion reduces to and we will test if our cylindrical distribution Ω(r) can satisfy it.

Concentric Spheroids
Although we are not concerned here with the spheroids, they will be touched on briefly in order to compare with earlier results about the stability of compound confocal spheroids (Cisneros et al. 2000). It is the only case known to us that allows solid body rotation. Certainly, the total potential in the atmosphere's surface depends on the elliptic coordinate λ being constant because of confocality, and thus the potential is a linear function of r (see Montalvo et al. 1983). Hence, by Equation (8), the angular velocity must be constant.
Here, it is not supposed that there is any relation between the shape of the nucleus and that of the atmosphere, except that they are concentric. Therefore, in general, λ varies from point to point in the envelope's surface, and the angular velocity will not be constant, i.e., the body rotates differentially.
We have established the frequencies for some particular series. For example, a series of models with an oblate nucleus for ε=1.0, e 1 =0.5, and e n =0.1 begins at e a =0.1 (contact model) and ends at e a =0.9, and has the frequencies  One sees that there are complex frequencies σ 4 and σ 5 whenever e a 0.3, that is, the model is unstable for a very flattened nucleus and atmosphere. A series of models with an oblate nucleus for ε=1.0, e 1 =0.5, and e n =0.3, i.e., a larger core, begins with e a =0.3 and ends with e a =0.9, and has the frequencies  The whole series is, this time, stable.

Figures with e > 0
We have computed the transverse-shear modes (σ 1 , σ 2 , σ 3 ), the toroidal modes (σ 4 , σ 5 ), and the pulsatory mode σ 6 for several series. ε was varied from small (∼1) to large values (∼100). d n and d a are allowed to take positive and negative values. e 1 , e n , and e a were chosen to obtain series with a prolate or oblate nucleus. From the comprehensive series set, we select just a few. They are not representative of all cases and are presented as examples only. Since s p n r = G 2 a is a normalized quantity, it is preferable to display results for specific cases, when we know the concrete value of the density ρ a .
A series of models with a prolate nucleus for ε=0.2, e 1 =0.5, and e n =0.48 (z Mn =0.52), d n =d a =−1/8, begins with e a =0.48 (contact figure) and ends with e a =0.63, and has the frequencies  All models in the series are stable, and near e a =0.54 is a neutral mode for σ 4 . A series of models with an oblate nucleus for ε=0.3, e 1 =0.3, and e n =0.1 begins with e a =0.1 (contact model) and ends with e a =0.7, and has the frequencies  Unstable models result when e a 0.2.

Figures when e < 0
It is believed that such models are physically impossible, because the core is lighter than the envelope (ρ n <ρ a ):the Rayleigh-Taylor instability prevents them from existing. 3 Nonetheless, they exist. We found in our previous work (Cisneros et al. 2017) that these models are in the equilibrium state, because the fundamental equations of motion are satisfied. (Self-)gravity is surely present and is balanced by pressure and (differential) rotation. Rotation plays an important role in maintaining equilibrium; pressure alone cannot counterbalance gravitation. Moreover, there is evidence that in the case of Rayleigh-Taylor instability, rotation tends to reduce it (Baldwin et al. 2015;Scase et al. 2017). By the way, the term instability is applied here to the equilibrium condition.
We have established the frequencies for some particular series. For example, a series of models with a prolate nucleus for ε=−0.3, e 1 =0.5, and e n =0.54 (Z Mn =0.58), is a normalized frequency. In this series, all models are stable (locally too). The lowest frequency σ 4 =0.02 is near the neutral mode. For other ε values, the frequencies behave in a similar way. However, they decrease a little as ε becomes more negative.

A Jupiter Model
We will construct a model for Jupiter based on our heterogeneous two-layer figure. This is characterized by the following parameters:the outer layer's density ρ a , the core's density ρ n (or better, the relative density difference ε), the normalized envelope's semiaxis z Ma , the normalized core's semiaxes e 1 and z Mn , and the distortion parameters d a and d n of the atmosphere and nucleus. These seven quantities are established using known observed facts of Jupiter. One part is ==´- Equatorial radius 7.1492 10 m, Polar radius 6.6854 10 7 m, Period 9.92597 hr, Angular velocity 1.7584 10 rad s . The other important part of Jupiter's observed data is the known gravitational field obtained from the Juno mission. From this field, which is obtained in a very precise way, the gravitational moments are deduced. One can find them in Table  S3 of Bolton et al. (2017) and are given here as In a more detailed study (Folkner et al. 2017), the gravitational moments are estimated for the two Juno orbits individually and combined, as given above. For the first orbit (PJ1), Doppler data were taken in an 8.3 hr interval around the perijove, and in a 5.9 hr interval for the second (PJ2). The perijoves were near the pole, but at different places for the two orbits. The uncertainties in the estimated J 2n for PJ2 were lower than those for PJ1. Here, we used only the combined results. We will now use another quantity, Jupiter's mean density ρ m , which is not directly taken from observation, because it depends on the model used, commonly a spheroid, to establish the volume. We think that it is not too far from reality, as Jupiter's shape does not differ greatly from the spheroid. We use the value r = -( ) 1326.5 kg m . 46 m 3 In our last paper (Cisneros et al. 2017), we built analytical expressions for the quadrupole and octupole in terms of the model's parameters. We also showed the relation with J 2 and J 4 . However, we proceed here directly and use for J 2n the definition where r and ϑ are spherical coordinates, and P 2n are 2n-order Legendre polynomials. We evaluate the integral more conveniently in cylindrical coordinates, by noticing that ρ(r)=ρ a if r is in the atmosphere, and ρ(r)=ρ n if r is in the nucleus: ò ò r e j j = -  In the Appendix, we give the expressions for these integrals in terms of the parameters of our model.  ) and density ρ a and the relative density difference ε r e = = ( ) 529.8, 285.04. 56 a As we see, the observed J 2n are acceptably well reproduced by our model, because they are within the range of uncertainty (actually, J 6 is outside the uncertainty interval; however, it almost touches it).
Our two-layer model explains reasonably well the observed gravitation of Jupiter and gives us insight into its internal structure. Nonetheless, it is possible that the interior might be more complicated and better described by a multilayer model. For this reason, we will speak of a mean outer region, instead of an atmosphere, and a mean central region, instead of a nucleus. According to it, Jupiter has an extended outer part of mean density ρ a =529.8, and a small, highly condensed, central portion with mean density ρ n ≈286ρ a . The model gives one additional information. Using Equations (54) and (56) as parameters, we ask whether the equilibrium state is possible, i.e., if there exist an angular velocity distribution for the core and envelope necessary to sustain equilibrium. We find numerically that indeed this is the case. The numerical results for Ω a and Ω n can be reproduced approximately by the expressions where ω a and ω n are properly the angular velocities of the atmosphere and nucleus. In the envelope's pole and equator, they are w w = = 0.000181245, 0.000170567.

ap ae
Hence, Jupiter's outer and central parts rotate differentially. According to our model, the mean angular velocity in Jupiter's surface is ω am =(ω ap +ω ae )/2=0.000175906 and corresponds to the period T m =2π/ω am =9.92195 hr. The periods for the pole and equator are T p =9.62967 hr and T e = 10.2325 hr. As we see, the calculated mean period is pretty nearly the observed one. Regarding Jupiter's central part, which cannot be seen, it rotates differentially, substantially faster than the outer part, with the angular velocities in the pole and equator ω np =0.00402584, ω ne =0.00313825, corresponding to periods of T p =26.0119 minutes and T e =33.3688 minutes. Moreover, the inner region deviates significantly from a spheroidal (d n =1.36706), whereas the exterior part is nearly a spheroid (d a =0.00186; Equation (49)).
The core's angular velocity is really high, and one wonders if it can be sustained by self-gravitation and pressure, i.e., if the equilibrium is stable. Hence, the next obvious step in building our model concerns its stability. Possessing the equilibrium model, we proceeded to calculate the lower six proper frequencies using the technique described above. Thus, all modes are real numbers, and our model is stable.
Since ω n is high anyway, one can ask how near it is to instability. We have artificially raised it by 4% and recalculated the proper modes:σ 4 and σ 5 became complex, meaning instability. This mode changes axial symmetry into a triaxial one (see Cisneros et al. 2015), in the approximation used (Equation (11)). Admitting our model, this means that, if for any reason the angular velocity of Jupiter's central part would increase, it might become unstable.

Discussion
Homogeneous Jupiter models cannot explain the observed properties. For example, the spheroid (e 1 =1, d=0) has the gravitational moments, according to the limit In both cases, the gravitational moments deviate significantly from the ascertained ones (see Equation (45)), with those of the spheroid being worse than those of the spheroidal figure. In the last, J 2n approximates the observed one better than the former for higher n.
Compound models like those of Hubbard (2012Hubbard ( , 2013 and Debras & Chabrier (2017) come much closer to Jupiter's gravitational field, besides the fact that their refined analytic procedures have merit in themselves. We feel that with such models, one can have deeper insight into the mass and angular velocity distribution in Jupiter's internal structure.
Our model predicts Jupiter's differential rotation, ω a , being nearly constant, with periods at the pole and equator of 9.63 and 10.23 hr and a mean period of 9.92 hr, the last being practically the same as that usually reported. Although the difference of 0.60 hr is small, it can, in principle, be measured.In this way, we could know if Jupiter does rotate differentially.
With this simple two-layer model, we do not claim to explain the real internal structure of Jupiter. We have worked with a plain ideal fluid mass, ignoring the chemical composition, state equations, temperature, magnetic field, and so on. The model was worked out in a complete fashion that fulfills the equilibrium requirement and stability for obtaining a geometric configuration that explains Jupiter's gravity. Moreover, regarding the number of layers, more realistic examples would take three (Guillot 1999), four (Guillot et al. 2018), or even five shells (Gudkova & Zharkov 1999). In a future paper, a generalization to n layers of our model will be presented, anticipating that from three layers on, we get exactly the moments J 2n and having each layer a determinate size.
Our technique is useful only to establish J 2n (like Hubbard 2013), imposing the restriction of equilibrium (and stability). To build models to study Jupiter's interior, taking as a basis models like those here presented, we must bear in mind that magnetic fields are present, and that fluids have in general a certain viscosity. This is important in our case, because of the existence of differential rotation; moreover, the rotation state is not the same for all shells. New results show that magnetic fields can smooth out the differential rotation to convert it to solid body rotation (Guillot et al. 2018).
In terms of the normalized variables, J 2n can be constructed by means of the functions (see Equation (48)