1. Introduction

In ^{Cisneros et al. I (2015}) equilibrium figures were
obtained from a self-gravitating homogeneous liquid mass endowed with an internal
motion of differential vorticity, whose surface equation was that of either a
distorted ellipsoid or a distorted spheroid; however, as we noticed *a
posteriori*, the figures, called ellipsoidal and spheroidal after Jeans
(1920), were erroneously deduced, although their qualitative features remained more
or less unaﬀected (see § 5.1). In the current paper, we investigate the equilibrium
of a heterogeneous model consisting of two concentric homogeneous spheroidal masses
of diﬀerent density, the ‘nucleus’ and the ‘atmosphere’ (whose semi-axes are not
correlated) through a simple and precise rotation law. Aﬃxing sub indexes
*n* and *a* to the pertinent quantities, the
body’s relative density, (*ρ*_{n}*−
ρ*_{a})/*ρ*_{a}
, is denoted by *ε*, with *ρ*_{n
}*> ρ*_{a}, so that the distribution
of density is somewhat closer to that prevailing in real stars. The quantities in
our development are normalized as explained in ^{Cisneros et al. I (2015)}.

2.Bernoulli’s theorem, the classical homogeneous figures and ours

Before going on with our heterogeneous model, we wish to discuss Bernoulli’s theorem as is commonly employed for obtaining the classical homogeneous figures (^{Lyttleton 1953}; ^{Dryden 1956}), particularly the Dedekind ellipsoids which, along with our models, are static. The steady-state equations of motion for a self-gravitating fluid are

where **v** is the velocity field, *V* is the potential and *p* is the pressure. These equations can also be written as

For a streamline (tangent to **v**) we get, after taking the dot product of **v** in the first equation (2),

that is

*k* being constant on a streamline. Commonly, *k* changes from one streamline to another, being an overall constant only when rot **v** = 0. The changes in *k* are dictated by [cf. equations (2) and (4)]

The Maclaurin and Jacobi figures are characterized by the rotational velocity field

where ω is the constant angular velocity. We can see that relation (6) satisfies the continuity equation div **v** = 0. According to equation (6), the *k* in Bernoulli’s equation is given by

that is

From which we deduce

*c* being any constant. Therefore, Bernoulli’s equation (4) becomes

which agrees with the well-known equation.^{1}

While the Jacobi and the Maclaurin figures rotate as a rigid body, a Dedekind ellipsoid (1, e_{2}, e_{3} are the major, middle, and minor semi-axes) remains fixed in space, its equilibrium being due to an internal motion of uniform vorticity ζ, its velocity field being

We can convince ourselves that V is tangent to ellipses x^{2} + y^{2}/e^{2}_{2} = const.. Hence, streamlines are ellipses perpendicular to the z-axis. For ζ constant, the velocity field satisfies the continuity equation, and *k* obeys the relation

so that

The solution of the partial diﬀerential equations (13) is

where c is an arbitrary constant. In this case, Bernoulli’s equation

reduces to

This is essentially equation (10) for obtaining the Jacobi ellipsoids if we put

a known result that establishes the equivalence between Jacobi’s angular velocity and Dedekind’s vorticity (^{Chandrasekhar 1969}). Of course, if we had assumed from the beginning that *k* = const., it would have been impossible to reach equation (15).

2.1. General Case with Cylindrical Symmetry

Each fluid point rotates, now, with non-constant angular velocity ω. The velocity field is given again by

but the continuity equation leads, presently, to

or

which means that, generally, the angular velocity must be a function of the kind

Assuming, additionally, that the velocity field is symmetric about the z-axis, we can express it as

Here, it is more convenient to use cylindrical coordinates (R, φ, z), and so the problem is independent of one coordinate: φ. In this system, the velocity field (tangent to circles) has a φ-component alone:

Equation (5) will have two terms only:

Making the variable change

equations (22) become

From the last of equations (24), we deduce

where *f* is an arbitrary function. Therefore, the first equation (24) implies that

or

i. e., the angular velocity can be at most a function of R alone, and the same for k:

In other words, the angular velocity distribution has cylindrical symmetry. Since *k* does not depend exclusively on *z*, disk-like distributions are impossible, contrary to what we stated in Paper II. Substituting k of equation (27) into Bernoulli’s equation (4), we obtain

Since the function *f* depends on *R* but not on *z* (i. e., it is constant on cylinders), we can determine it using only the surface equation, on which *p* = 0:

where *z* is a function of *R* for a figure with cylindrical symmetry. Equation (29) allows to determine the function *f* (= −*V*), and with equation (26) ω is established:

which is the general angular velocity distribution law for any axial-symmetric, incompressible, self-gravitating fluid with *p* = 0 on its surface. Equation (30) is, essentially, Newton’s second law for a unit mass particle:

Using the familiar variable r (= x^{2} + y^{2} = R^{2}), equation (30) can be written as

In the special case of a Maclaurin spheroid, the potential can be expressed as

At any interior point; on the surface

Hence, according to equation (31), the angular velocity is given by

which agrees with the well-known result. Taking into account that

and that υ_{1} = 2 π − υ_{3}, we can also write

2.2. the homogeneus spheroidal mass

The surface equation of the homogeneous spheroidal mass (the ellipsoidal mass will not be considered here) is

where *d* is a parameter larger than −1/4; the rotation semi-axis is not
*e _{3}*, but

For establishing the angular velocity at equilibrium (31), we must take the derivative of the potential with respect to *r* at the surface. Since *V* is only known numerically, this process has to be carried out numerically. Nonetheless, we can use another approach. We approximate the potential by a polynomial in *r* of the form

so that the angular velocity becomes

For *d* not too large, the mean absolute error in *V* is about 10^{−7}, which grows with increasing *d* (for *d* = 2, the mean error is 10^{−4}).

3. The heterogeneous spheroidal mass

The composite model consists of a core, or nucleus, of density ρ_{n}, whose shape
is

surrounded by an envelope, or atmosphere, of density ρ_{a}, of the form

Here *e _{1}* is the ratio of the nucleus and atmosphere major axes;

*e*and

_{n}*e*are proportional to the rotation semi-axes, which are

_{a}

The net potentials at any point of the nucleus and the atmosphere are (^{Montalvo et al. 1983}),

respectively, where *V _{nn}* is the self-potential of the nucleus;

*V*is the potential on the nucleus due to the atmosphere;

_{na}*V*is the potential on the atmosphere due to the nucleus;

_{an}*V*is the self-potential of the atmosphere; ε is the body’s relative density difference:

_{aa}

Let us now suppose that core and envelope rotate with (variable) angular velocities
*ω _{n}* and

*ω*, so that their velocity fields are

_{a }

Both fields must obey the continuity equation div **v** = 0. Thus, according to
equation (20), we have

In the equilibrium state the angular velocities (42) must fulfill Bernoulli’s equation:

where *f _{n}* and

*f*are functions to be determined.

_{a }*f*and

_{n}*f*are established from the surface conditions

_{a}*p*=

_{n}*p*(core surface) and

_{a}*p*=0 (envelope surface), where

_{a}*p*and

_{n}*p*are pressures on points of the nucleus and the atmosphere:

_{a}

Using equation (31), these relations can be written as

where *r = x ^{2}*

*+ y*and

^{2}*V*,

_{N}*V*are potentials on core and envelope surfaces, taken as functions of

_{A}*r*only. Thus, for establishing

*ω*and

_{a}*ω*, the potential derivatives must be available at each point of the surface. Since the potential is given numerically, its derivative must be numerically computed. For a homogeneous body (only interior points are needed), we find that a reasonable good approximated expression for the potential at points on the surface is:

_{n}

where *α _{i}* are parameters to be fixed by fit procedures, and thus
the angular velocity is given by equation (35)

In the case of our heterogeneous mass, equation (46) is not acceptable, especially at points where external potentials are needed (they must approach 0, as r →∞). To obtain equilibrium figures for the heterogeneous case, the potential derivatives will be established by numerical means. For this purpose, we use the approximation

where *h* is a small quantity. The precision of expression (47) is of order
*h ^{2}*. A not too small value for

*h*will be used, since potentials are calculated with an accuracy of about 10

^{−7}; a reasonable

*h*could be ≈ 5 × 10

^{−3}.

4. Model geometry

The heterogeneous mass can easily be built by making the two surfaces similar: coordinates of corresponding surface points are proportional, so that

for all (*x, y, z*), from what it follows that *c* =
1/*e _{1}*,

*e*=

_{n}*e*

_{1}*e*and

_{a}*d*=

_{n}*d*. Thus, the whole mass is characterized by only two parameters, say,

_{a}*e*and

_{a}*d*(and, of course, the relative size

_{a}*e*of major semi-axes, and the relative density diﬀerence ε). Hence, the series, should they exist, would be relatively simple to handle. Yet we prefer to explore a more general case, and let

_{1 }*e*,

_{n}*e*,

_{a}*d*,

_{n}*d*vary freely, thus leading to more involved series (actually, six-parametric, if we take into account the parameters

_{a}*e*and ε). The model eﬀectively yields series which are subject to certain limits, both physical and geometrical; for instance, the major semi-axis of the atmosphere must be greater than that of the nucleus:

_{1}*e*< 1; besides, the rotation semi-axes

_{1}*z*,

_{Mn}*z*must be smaller than the corresponding major semi-axes:

_{Ma}

or [cf. equation (38)]

Finally, the rotation semi-axis of the atmosphere must not be smaller than that of the nucleus, otherwise the atmosphere would intrude in the nucleus, and the present equilibrium conditions would not apply:

The particular configuration for which z_{M n} = *z*_{Ma}
(*e _{1}* ≠ 1), i. e., when the poles coincide, is
termed a ‘contact figure’;

*e*and

_{a}*e*are related by

_{n}

5. Numerical results

5.1. The homogeneous spheroidal mass

Clearly, our calculations of (^{Cisneros et al. 2015}, Paper I) must be aﬀected if the right *k* in Bernoulli’s equation (4) is used, but the qualitative features of the new figures remain more or less alike. In the case of homogeneous figures for which the surface equation is

we come as before to continuous *e _{3}* -series for each

*d*value and, furthermore, we again find limits. For

*d*positive, and low

*e*-values up to a definite limit,

_{3}*ω*

^{2}increases from pole to equator; thereafter, the tendency is inverted (Figure 1). For

*d*negative (but greater than −1/4), another limit appears: as d becomes more and more negative, the

*e*-series becomes shorter, because

_{3}*ω*

^{2}takes negative values, and the figures come into a forbidden region (Figure 1), a behavior that was also noticed in previous work.

5.2. The Heterogeneous Spheroidal Mass

The present series are more involved than those for the homogeneous case, so we must proceed
care-fully and not try to get an overwhelming bulk of models that would render
it diﬃcult to give a clear panorama of their regularities. For this purpose, a
basic series is constructed by fixing five parameters: ε= 1,
*e _{1}* = 0.5,

*e*= 0.1 and

_{n}*d*=

_{n}*d*= −1/8, while the remaining one e

_{a}_{a}is let to vary in its allowed range.

5.2.1. The Case d_{n}= d_{a}= −1/8

First, we take *d _{n}*

*= d*= −1/8 and allow

_{a}*e*to vary from 0.1 up to its maximum value, finding that equilibrium is possible for the angular velocity distribution given by equation (31) for the envelope and core. In Figure 2 we plot the angular velocity distribution in the nucleus (R < 0.5), and the atmosphere (R < 1); Figure 3 is a sketch of the model geometry.

_{a}

In Figure 2 (see also Figure 3), the series begins with a contact figure in which the
poles of core and envelope touch each other, but not the equators, since the
equator of the atmosphere is twice as large as that of the nucleus
(*e _{1}* = 0.5). We see that the core rotates
slower than its envelope, even at the pole. In both, the angular velocity
decreases from pole to equator, i. e., the central parts rotate faster than
the outer ones. As

*e*increases (the atmosphere expands),

_{a}_{a}-values,

_{a}= 0.78 (z

_{M}= 0.84). At the same time, less dramatic changes in the core occur with increasing e

_{a}: the pole

^{2}becomes zero at a point, and the other, less drastic, where the angular velocity reverses its tendency from monotonically decreasing to increasing. Certainly, the change takes place gradually, and we refer to a limit when ω

^{2}increases monotonically for the first time. Still another limit might be disclosed by building more series for larger e

_{n}values. Indeed, we find similar series starting with a contact figure and ending up with a top figure (e

_{a}≈ 0.78) if e

_{n}≤ 0.3975 (z

_{Mn}= 0.43). For the small interval 0.4004 > e

_{n}> 0.3975 there are series that do not be-gin with a contact figure and end with a one-member series at e

_{n}= 0.4004, e

_{a}= 0.969.

5.2.2. The Case d_{n} = d_{a} = 1/8

Here we take d_{n} = d_{a} = 1/8 and let e_{a} vary from 0.1
(contact figure) up. Once more, we get a set of series somewhat diﬀerent
from the former one (Figure 2), but
also sharing several properties:

The series begins with a contact figure with the atmosphere rotating faster than the nucleus, as formerly.

Both diﬀerential angular velocities always monotonically decrease from pole to equator. In other words, there are not transition figures.

The limiting figure occurs at e

_{a}= 0.967 (z_{Ma}= 0.917).The contact figure has a limit for e

_{n}= 0.398 (cf. § 2.1.1).The series do not end at e

_{n}= 0.398; they continue for larger e_{n}and are more and more narrow, ending with an isolated figure with e_{n}= 0.4827, e_{a}= 0.969 (z_{Mn}= 0.457) (see Figure 4).

Hence, the plus sign of d_{a}, d_{n} hinders the monotonous property of the
angular velocity, and raises the upper limit of e_{a}. Additionally,
the series does not disappear after a last contact figure, but at a higher
e_{n} limit.

5.2.3. The Case d_{n}= 1/8, d_{a}= −1/8

This time, only d_{n} sign is modified, so that our set of parameters is

Once more, the series starts with the contact figure e_{a} = 0.0876; e_{a}
is established by z_{Ma} = z_{Mn} condition (equation (51)).
The series resembles more the d_{a} = d_{n} = −1/8 case than
the d_{a} = d_{n} = 1/8 one, in the sense that it has
angular velocity transition limits and a limiting contact figure at
e_{n} ≈ 0.4 (z_{Mn} = 0.38). For greater e_{n}
values, there are, however, more series without an initial contact figure,
which become narrower as e_{n} → 0.5 (z_{Mn} = 0.47). For
e_{n} = 0.5 there only exists a series with the member
e_{a} = 0.97 (z_{Ma} = 0.92).

5.2.4. The Case d_{n} = −1/8, d_{a} = 1/8

Here, d_{a} alone is modified relative to first case, so that the parameters values are now

With these parameters, we also establish series starting with a contact figure, when
e_{n} < 0.4 (z_{Mn} = 0.43). Series without a contact
figure are possible when 0.4 < e_{n} < 0.421; for
e_{n} = 0.421 (z_{Mn} = 0.46) the series has a unique
member at e_{a} = 0.97. The atmosphere does not present angular
velocity transition limits, and the nucleus has only one, at e_{n} =
0.2, e_{a} = 0.6.

5.3. Consequences of the ε-value

To study the ε impact on our series, we assumed values for the parameters according to

and constructed the series changing e_{a} for a given e_{n} of a set 0.1 . .
. 0.5, as formerly done. As an illustration, only cases d_{n} = ∓1/8,
d_{a} = ±1/8 were considered, and cases with d_{n} =
d_{a} = −1/8, 1/8 were disregarded. We refer briefly to these
d-values as −+, +−, −−, and ++ (first sign for d_{n}, second for
d_{a}). Generally speaking, we did not find any new property, (other
than those recognized above) when ε was modified, i. e., for each ε, the set of
series can or cannot show angular velocity distribution ‘turning points’
(transition to monotonically increasing distribution), there is a final series
having a contact figure, and there is a last (one-member) series. Clearly, as ε
changes, the values characterizing a last contact figure, a last series, and a
transition figure, must suﬀer modifications. The contact limiting figure has an
e_{n} -value that first decreases and then increases, as ε increases
from 1 to 24 (see Figure 5). On the other
hand, e_{n} in the last series continually decreases as ε increases. The
e_{n} span between last contact figure series and last series
initially increases and thereafter continuously decreases (Figure 5).

**5.4. Consequences of the e**
_{1}
**-value**

To study changes of our basic series (Figure 2)
regarding the nucleus and atmosphere relative size e_{1}, we varied it
in steps of 0.1 from 0.9 down. We fixed d_{a} = ±1/8, d_{n} =
±1/8, ε = 1 and constructed e_{n} series varying e_{a} for each
e_{1} value of the set. For example, when e_{1} = 0.9 we
built the series for e_{n} = 0.1, 0.2, . . . . Qualitatively again, we
found series with or without transition limits, e_{a} upper limits, last
contact figure for e_{n} series, and final one-member series.
Quantitatively, there are diﬀerences regarding the above results. The
ω^{2} magnitudes and the variation range change; the e_{n}
gap between last contact figure and last one-member series grows as
e_{1} decreases, starting at e_{1} = 0.5; likewise, beyond
e_{1} = 0.5 the gap approaches 0 (see Figures 6 and 7). This behavior is somewhat similar to the ε-eﬀect
(Figures 4 and 5): as e_{1} or
ε increases the gap be-comes narrower, and finally it disappears.