A THREE-DIMENSIONAL NUMERICAL SOLUTION FOR THE SHAPE OF A ROTATIONALLY DISTORTED POLYTROPE OF INDEX UNITY

Consider an isolated mass of gas with total mass M that is rotating rapidly about the z-axis with angular velocity $\Omega _0\hat{\textbf {\em z}}$. When the density ρ = ρ₀ is uniform, i.e., ρ₀ is constant, the rotating fluid body is in the shape of an axisymmetric oblate spheroid with eccentricity $\mathcal {E}_0$ given by the well-known equation (Lamb 1932):

$$\begin{equation} \frac{\Omega _0^2}{2\pi G\rho _0}= \frac{\sqrt{1-\mathcal {E}_0^2}}{\mathcal {E}_0^3}\big(3-2\mathcal {E}_0^2\big)\sin ^{-1}\mathcal {E}_0-\frac{3\big(1-\mathcal {E}_0^2\big)}{\mathcal {E}_0^2}. \end{equation} \tag{ 4 }$$

When the density ρ within a rapidly rotating body is non-uniform, it was shown analytically that the leading-order solution is also in the shape of an axisymmetric oblate spheroid (Kong et al. 2010). In this study, we assume that the effect of rotation is sufficiently strong that the rotating body is in the shape of an oblate spheroid with moderate or large eccentricity that is still within the stable limit.

With a polytrope of index unity, the pressure–density relation of the gas is given by

$$\begin{equation} p=K\rho ^2, \end{equation} \tag{ 5 }$$

where K is a constant. The mechanical equilibrium equation is

$$\begin{eqnarray} &&-\frac{\mbox{\boldmath $\nabla $}p}{\rho }-\mbox{\boldmath $\nabla $}V_g-\mbox{\boldmath $\nabla $}V_c=0\; \mbox{in}\; \mathcal {D}, \end{eqnarray} \tag{ 6 }$$

where V_g is the gravitational potential, V_c is the centrifugal potential, and $\mathcal {D}$ is the domain of the rotating body, subject to the boundary condition

$$\begin{eqnarray} &&p=0\; \mbox{at the bounding surface of } \;\; \mathcal {D}. \end{eqnarray} \tag{ 7 }$$

The shape of a rotationally distorted body is characterized by its eccentricity $\mathcal {E}$ defined as

$$\mathcal {E}=\frac{\sqrt{R_e^2-R_p^2}}{R_e},$$

where $0 < \mathcal {E}< 1$ and, R_e and R_p are the equatorial and polar radii of an oblate spheroid, respectively. In Equation (6), the gravitational potential V_g is given by

$$\begin{equation} V_g=-G\int\!\! \int\!\! \int _{\mathcal {D}}{\frac{\rho (\textbf {\em r}^\prime ){d}^{3}\textbf {\em r}^\prime }{\left|\textbf {\em r}-\textbf {\em r}^\prime \right|}} \end{equation} \tag{ 8 }$$

where $\textbf {\em r}$ is the position vector. The centrifugal potential V_c is

$$\begin{equation} V_c=- \frac{\Omega _0^2 c^2}{2}(1+\xi ^2)(1-\eta ^2). \end{equation} \tag{ 9 }$$

Here we have used oblate spheroidal coordinates, (ξ, η, ϕ), defined by the coordinate transformation with Cartesian coordinates

$$\begin{eqnarray*} x&=&c\sqrt{(1+\xi ^2)(1-\eta ^2)}\cos {\phi }, \nonumber \\ y&=&c\sqrt{(1+\xi ^2)(1-\eta ^2)}\sin {\phi }, \nonumber \\ z&=&c\xi \eta , \nonumber \end{eqnarray*}$$

where c is the focal length of an oblate spheroid with its bounding surface ${\cal S}$ described by

$$\xi =\xi _o=\sqrt{\frac{1}{\mathcal {E}^2}-1}.$$

The shape parameter, the size of eccentricity $\mathcal {E}$, is a priori unknown.

On the bounding surface ${\cal S}$ of a rotating body, the free-surface condition must be imposed at the equilibrium

$$\begin{equation} V_g+V_c= \mbox{constant } \mbox{at the outer surface} \; {\cal S}. \end{equation} \tag{ 10 }$$

Eliminating p from Equations (5) and (6) yields the equation

$$\begin{equation} 2K\mbox{\boldmath $\nabla $}\rho +\mbox{\boldmath $\nabla $}V_g+\mbox{\boldmath $\nabla $}V_c=0, \end{equation} \tag{ 11 }$$

which, after applying $\mbox{\boldmath $\nabla $}\mbox{\boldmath $\cdot $}$ to the both sides, becomes

$$\begin{equation} K\mbox{\boldmath $\nabla $}^2\rho +2\pi G\rho =\Omega _0^2. \end{equation} \tag{ 12 }$$

After scaling Equation (12) with the length scale R_e and the mass scale M, we obtain the governing equation in the dimensionless form

$$\begin{equation} \mbox{\boldmath $\nabla $}^2\rho +\alpha \rho =\beta , \end{equation} \tag{ 13 }$$

where α and β are defined as

$$\begin{equation} \alpha = \frac{2\pi G R_e^2}{K} \;\; \mbox{and} \;\; \beta =\frac{\Omega _0^2R_e^5}{MK}. \end{equation} \tag{ 14 }$$

The boundary condition (7), after using the EOS, becomes

$$\begin{equation} \rho = 0\; \mbox{at the outer surface} \; {\cal S}. \end{equation} \tag{ 15 }$$

Since α is always positive, Equation (13) represents an inhomogeneous Helmholtz equation (Chandrasekhar 1933). The condition of the free surface in Equation (10) can also be written in the dimensionless form

$$\begin{eqnarray} &&\int\!\! \int\!\! \int _{\mathcal {D}}{\frac{\rho (\textbf {\em r}^\prime ) {d}^{3}\textbf {\em r}^\prime }{\left|\textbf {\em r}-\textbf {\em r}^\prime \right|}} + \frac{\beta c^2}{4 \alpha }(1+\xi ^2)(1-\eta ^2)\nonumber\\ &&\quad = \mathrm{constant} \;\; \mbox{at} \;\; \xi =\sqrt{\frac{1}{\mathcal {E}^2}-1}, \end{eqnarray} \tag{ 16 }$$

which must be satisfied at the equilibrium. Our primary task is to solve the inhomogeneous Helmholtz equation (13) for both the shape parameter $\mathcal {E}$ and the density distribution ρ(ξ, η) that satisfy both Equations (15) and (16).

When the shape parameter $\mathcal {E}$ and the density distribution ρ(ξ, η) become available for a highly flattened planet, we can readily compute the total mass of the rotating body

$$\begin{equation} M= \int \!\!\int\!\! \int _{\mathcal {D}} \rho (\textbf {\em r}^\prime ) {d}^{3}\textbf {\em r}^\prime , \end{equation} \tag{ 17 }$$

and the exterior gravitational potential V_g which is expanded in terms of spherical harmonics P_2n ,

$$\begin{equation} V_g(r,\theta )=-\frac{GM}{r} \left[\!1-\sum _{n=1}^{\infty }J_{2n}\left(\frac{R_e}{r}\right)^{2n}P_{2n}(\cos \theta )\!\right], \;\; r \ge R_e, \end{equation} \tag{ 18 }$$

where (r, θ, ϕ) are spherical polar coordinates with θ = 0 being at the axis of rotation, and J₂, J₄, J₆, ⋅⋅⋅, are the zonal gravitational coefficients. It is important to point out that Equation (18) is only valid for r ⩾ R_e where there is no mass and the potential satisfies Laplace's equation. This is because the region R_p < r < R_e has empty space in part and mass near low latitudes and, consequently, the potential satisfying Poisson's equation with the density is a highly complicated function of latitude and radius.

A three-dimensional finite-element method is employed to solve the shape and internal structure of a non-spherical body through the following three different stages. The first stage is to construct a three-dimensional finite-element mesh by making a tetrahedralization of an oblate spheroid with a guessed eccentricity $\mathcal {E}_G$. This is because we do not have a priori knowledge of the shape of a rapidly rotating liquid planet. A sketch of the finite-element mesh for an oblate spheroid with $\mathcal {E}=0.5$ is illustrated in Figure 1(a). In comparison to a spectral or finite difference method, the finite-element method is free of the pole and central numerical singularities.

Zoom In Zoom Out Reset image size

In the second stage, we solve Equation (13) for ρ together with the boundary condition (15) with a given α and β for a guessed eccentricity $\mathcal {E}=\mathcal {E}_G$. A Galerkin-weighted residual approach is adopted in the finite-element formulation of Equation (13). Multiplying Equation (13) by the corresponding weight functions w_ρ, then integrating the resulting equation over the spheroidal domain $\mathcal {D}$, after making use of integration by parts, we derive the weak formulation of Equation (13)

$$\begin{equation} \int _{\mathcal {D}}{-\mbox{\boldmath $\nabla $}\rho\, \mbox{\boldmath $\cdot $}\,\mbox{\boldmath $\nabla $}w_\rho \,{d}V + \alpha \int _{\mathcal {D}}\rho w_\rho } {d}^{}V = \beta \int _{\mathcal {D}}{ w_\rho } {d}^{}V, \end{equation} \tag{ 19 }$$

where the boundary integral vanishes as the weight function w_ρ is zero on the bounding spheroidal surface $\mathcal {S}$. Within each tetrahedron, as shown in Figure 1(b), 10 nodes are used in representing the density ρ

$$\begin{equation} \rho =\sum _{j=1}^{10}\rho _j \Phi _j, \end{equation} \tag{ 20 }$$

where ρ_j is the value of the density ρ on the jth node in the element and Φ_j is the quadratic function defined as

$$\begin{eqnarray} &&\Phi _1=L_1(2L_1-1), \Phi _2=L_2(2L_2-1),\nonumber \\ &&\Phi _3=L_3(2L_3-1), \Phi _4=L_4(2L_4-1),\nonumber \\ &&\Phi _5=4L_1L_2, \Phi _6=4L_1L_3, \Phi _7=4L_1L_4,\nonumber \\ &&\Phi _8=4L_2L_3, \Phi _9=4L_2L_4, \Phi _{10}=4L_3L_4, \end{eqnarray} \tag{ 21 }$$

where L_j , j = 1, ⋅⋅⋅, 4 are the volume coordinates of the finite element and Φ_j has the properties

$$\begin{equation} \Phi _j(\textbf {\em r}_i)=\delta _{ij}, \;\; \sum _{j=1}^{10}\Phi _j=1, \end{equation} \tag{ 22 }$$

where $\textbf {\em r}_i$ is the position vector for the ith node in the element. The weight functions are selected to be the same as the corresponding shape functions. Applying a standard procedure of the finite-element method, we obtain a system of linear equations for the coefficients ρ_j on the whole three-dimensional mesh, which is then solved by a Krylov subspace iterative method. The values of ρ_j on the mesh provide a numerical solution ρ(r) that satisfies Equations (13) and (15). For the results reported in this paper, a spheroidal domain is typically divided into about 10⁶ ∼ 2 × 10⁶ tetrahedral elements.

It is critically important to note that a numerical solution found in the second stage does not, in general, satisfy the free-surface condition (16). Although condition (16) is mathematically simple, it is numerically problematic. In order to find a numerical solution satisfying Equation (16), we introduce an auxiliary function, || dV_t / dη||₂, defined as

$$\begin{eqnarray} \left|\left|\frac{\,{d}V_t}{\,{d}\eta } \right|\right|_{2} &=& \frac{1}{2 \pi } \int\!\! \int _{\mathcal {S}} \Bigg| \frac{\partial }{\partial \eta } \left[\int \!\!\int\!\! \int _{\mathcal {D}}{\frac{\rho (\textbf {\em r}^\prime ){d}^{3}\textbf {\em r}^\prime }{\left|\textbf {\em r}-\textbf {\em r}^\prime \right|}}\right. \nonumber\\ &&\left. + \frac{\beta c^2}{4 \alpha }(1+\xi ^2)(1-\eta ^2) \right]_{\xi =\xi _o} \Bigg|^2 \,{d}\mathcal {S}, \end{eqnarray} \tag{ 23 }$$

where $\int \int _\mathcal {S} \,{d}\mathcal {S}$ represents the surface integral over the bounding surface $\mathcal {S}$ of the oblate spheroid. Theoretically, we would expect that

$$\left|\left|\frac{\,{d}V_t}{\,{d}\eta } \right|\right|_{2} >0$$

in the neighborhood of an equilibrium while

$$\left|\left|\frac{\,{d}V_t}{\,{d}\eta } \right|\right|_{2} =0$$

at the equilibrium. Far away from the equilibrium, we would have

$$\left|\left|\frac{\,{d}V_t}{\,{d}\eta } \right|\right|_{2} =\mbox{O}(1).$$

Numerically, however, the mathematical condition || dV_t /dη||₂ = 0 cannot be exactly satisfied at the equilibrium. Instead, || dV_t / dη||₂ = 0 will be replaced by a numerical condition that || dV_t / dη||₂ at the equilibrium reaches a minimum that is generally non-zero.

The third stage is an iterative procedure that repeats the first and second stages, such that condition (16) is approximately satisfied. In other words, by computing many numerical solutions for different values $\mathcal {E}_G$ at fixed K, α, and β, we are able to determine a particular set of ρ and $\mathcal {E}$ that satisfies not only Equations (13) and (15) but also Equation (16). All the numerical computations presented in this paper are fully three-dimensional, but all the numerical solutions turn out to be axisymmetric (∂/∂ϕ = 0). The proposed three-dimensional method is potentially suited for describing general properties of three-dimensional deformations of both rotationally and tidally distorted bodies.

In modeling a rapidly rotating gaseous planet or star, some parameters will be regarded as being well determined by observations while other parameters have to be treated as unknown. In the case of Jupiter, we regard the equatorial and polar radii at the one-bar surface R_e and R_p —which are R_e = 71, 492 km and R_p = 66, 854 km (Seidelmann et al. 2007)—as the known parameters. They yield an eccentricity $\mathcal {E}_J=0.3543$ for the shape of Jupiter. We also regard the rotational period T_J of Jupiter, T_J = 9.925 hr (Seidelmann et al. 2007), as a known parameter. This corresponds to the angular velocity Ω_J = 1.7585 × 10⁻⁴ s⁻¹ for the rotation of Jupiter. The two-dimensional distributions of the density ρ(ξ, η) and the pressure p(ξ, η) within Jupiter are regarded as unknown functions to be determined via our three-dimensional numerical calculation. In addition, the size of K in the EOS p = Kρ² is also unknown.

The objective of our numerical calculation is, via a hybrid inverse method, to find the density ρ(ξ, η) and the pressure p(ξ, η) by matching K in the EOS p = Kρ² with the observed shape $\mathcal {E}_J$ and the observed angular velocity Ω_J . For a trial value K = K_G and a fixed Ω_J , we calculate, through the three-stage iterative scheme, an equilibrium solution ρ_G satisfying Equation (13) as well as Equations (15) and (16) at a particular value $\mathcal {E}_{G}$ at the equilibrium. However, $\mathcal {E}_{G}$ is generally inconsistent with the observed value $\mathcal {E}_J$ for Jupiter. By repeating this process at many different values of K_G , i.e., through an iterative scheme with K, we are able to find a particular value K_J = K_G such that the shape parameter $\mathcal {E}_{G}$ is consistent with the observed value $\mathcal {E}_J$. In other words, a numerically expensive two-parameter ($\mathcal {E}-K$) iterative process is required in determining the rotationally distorted internal structure of Jupiter.

A relationship between K and $\mathcal {E}$ at the equilibrium, resulting from the two-parameter ($K-\mathcal {E}$) iterative process, is shown in Figure 3(a) where the star symbol denotes a particular numerical solution that matches the numerical model with the observed shape of Jupiter $\mathcal {E}=\mathcal {E}_J$ when K = K_J = 200, 102 Pa m⁶ kg⁻². In other words, for the fixed rate of rotation Ω_J and a polytrope of unit index n = 1, the EOS for Jupiter must be of the form

$$p= 200,102 \,\mathrm{Pa} \,\mathrm{m}^6 \,\mathrm{kg}^{-2} \, \rho ^2$$

in order that its shape matches the observed value $\mathcal {E}=\mathcal {E}_{J}$. The detail of the three-stage iterative procedure at K = K_J is depicted in Figure 3(b). It shows that the function || dV_t / dη||₂ at fixed K = K_J decreases from || dV_t / dη||₂ = O(1) far away from the equilibrium to || dV_t / dη||₂ = O(10⁻³) at the equilibrium when reaching the minimum at $\mathcal {E}=\mathcal {E}_J$. Each triangle in Figure 3(b) represents a three-dimensional solution at the fixed K_J using different three-dimensional finite-element meshes. The particular solution at $\mathcal {E}=\mathcal {E}_J$ and K = K_J represents a numerical model for Jupiter with a polytrope of unit index n = 1. Figure 4 shows the density and pressure distribution of this numerical solution in the equatorial plane as well as the eccentricity of constant density surfaces while Figure 5 depicts the two-dimensional density distribution, ρ(η, ξ), in a meridional plane.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The total mass M_J of Jupiter is M_J = 1.8986 × 10²⁷ kg (Williams 2010), giving rise to the mean density $\bar{\rho }_J=1.326\times 10^{3}\,\mathrm{kg\,m}^{-3}$. It needs to be emphasized, however, that the Jupiter model with a polytrope of unit index n = 1 does not require an a priori value for the planet's mass. Rather, total mass is determined as part of the three-dimensional numerical solution. This provides a useful constraint on the model whose predicted mass can be compared with the observed mass of Jupiter M_J . For the two-dimensional density distribution ρ(ξ, η) shown in Figure 5, we find the total mass of our Jupiter model is 100.4% of M_J . This difference could be attributable to the inadequacy of the n = 1 polytrope to represent the real EOS of the Jovian interior. The zonal gravitational coefficients J₂, J₄, J₆ predicted by the model using the expression (18) are listed in Table 1. Our model predictions of the zonal gravitational coefficients are in reasonably good agreement with observed values (Jacobson 2003) especially when considering that there is only one parameter in the model.

Our model, which assumes only the size and shape of Jupiter's surface and the planet's rotation rate together with an EOS of the form p = Kρ², does a reasonable job of predicting the low-order gravitational coefficient J₂ of Jupiter with a better than 2% accuracy and the planet's mass with about 4% accuracy. This opens the possibility of inferring the mass and gravitational coefficients of extra-solar bodies whose rotation rates, sizes, and shapes can be measured. Table 1 in the third column also gives values of the zonal gravitational coefficients J₂, J₄, J₆ computed from the theory of figures. They were obtained by an application of the Zharkov–Trubitsyn third-order theory of figures to Jupiter (Zharkov & Trubitsyn 1978), with a polytropic density of index one, an equatorial radius of 71,492 km, and the IAU System III rotation rate of ${\rm 9^h 55^m 29\buildrel{\mathrm{s}}\over{.}71}$.

The Appendix presents a simplified theory for estimating the value of the polytropic constant K. The theoretical estimate for K from Equation (A4) gives K = 200, 544 Pa m⁶ kg⁻². This estimate is quite close to K = 200, 102 Pa m⁶ kg⁻² obtained by our hybrid inverse method via three-dimensional numerical simulation.

Our non-spherical model is particularly suitable for studying the fast rotating B-type star, α Eridani (Harmanec 1988; Perryman et al. 1997), which is strongly distorted by rotational effects. Since the internal pressure in α Eridani is likely radiation-field-dominated, its mean polytropic index n would be larger than unity. We still adopt the polytropic index n = 1 because our main objective is to illustrate a new method of determining the rotationally distorted shape and internal structure from the known or deduced physical parameters of the rapidly rotating star.

Harmanec (1988) estimated, from the luminosity–mass relation of main sequence stars, that the mass of α Eridani is M_A = 6.07 M_☉(1.21 × 10³¹ kg). From spectral line broadening it is deduced that Achernar's rotational speed projected on its equator is V_eqsin i = 225 km s⁻¹ (Slettebak 1982), where i is the inclination angle between the rotational axis and the line of sight from the Earth which cannot be observed directly. According to model 4 of Carcifi et al. (2008), which best fits the interferometry observation, the equatorial radius of α Eridani is 10.2 R_☉(R_e = 7.10 × 10⁶ km). Using this radius and taking i = 65° (Carcifi et al. 2008), we obtain an angular velocity Ω_A = 3.493 × 10⁻⁵ s⁻¹. In our calculation for α Eridani, we regard its total mass M_A , its equatorial radius R_e , and its angular velocity Ω_A as the known parameters. The two-dimensional distributions of the density ρ(ξ, η) and the pressure p(ξ, η) within α Eridani are regarded as unknown functions to be determined via our three-dimensional numerical calculation. In addition, the shape parameter $\mathcal {E}$ and the size of K in the EOS p = Kρ² are also unknown.

The objective of our numerical calculation for α Eridani is, via a hybrid inverse method, to find the shape parameter $\mathcal {E}$, the density ρ(ξ, η), and the pressure p(ξ, η) by matching K in the EOS p = Kρ² with the deduced mass M_A . For a trial value K = K_G and a fixed Ω_A , we calculate, through the three-stage iterative scheme, an equilibrium solution ρ_G satisfying Equation (13) as well as Equations (15) and (16) at the equilibrium, which gives rise to a particular total mass M_G . However, this size of the mass M_G is generally inconsistent with the value M_A for α Eridani. By repeating this process at many different values of K_G , i.e., through an iterative scheme with K, we are able to find a particular value K_A = 1.8241 × 10⁹ Pa m⁶ kg⁻² such that its total mass M_G is consistent with the value M_A . A relationship between K and M at the equilibrium, resulting from the two-parameter (K−M) iterative process, is shown in Figure 6(a) where the star symbol denotes a particular numerical solution that matches the numerical model with the observed mass of α Eridani M_A at K = K_A = 1.8241 × 10⁹ Pa m⁶ kg⁻². In other words, for the fixed rate of rotation Ω_A = 3.493 × 10⁻⁵ s⁻¹ and a polytrope of unit index n = 1, the EOS for α Eridani must be of the form

$$p= 1.8241 \times 10^9 \,\mathrm{Pa} \,\mathrm{m}^6 \,\mathrm{kg}^{-2} \; \rho ^2$$

in order that its mass matches the value M_A = 6.07 M_☉(1.21 × 10³¹ kg). The detail of the three-stage iterative procedure at K = K_A = 1.8241 × 10⁹ Pa m⁶ kg⁻² is depicted in Figure 6(b). It shows that the function || dV_t / dη||₂ at fixed K = K_A decreases from || dV_t / dη||₂ = O(10) far away from the equilibrium to || dV_t / dη||₂ = O(10⁻²) at the equilibrium when reaching the minimum at $\mathcal {E}=\mathcal {E}_A=0.7222$. The particular solution obtained at $\mathcal {E}_A=0.7222$ and K_A = 1.8241 × 10⁹ Pa m⁶ kg⁻² represents a numerical model for α Eridani with a polytrope of unit index n = 1. Figures 7(a) and (b) show the density and pressure distribution of this numerical solution in the equatorial plane while the eccentricity of constant density surfaces within α Eridani is shown in Figure 7(c). Figure 8 depicts the two-dimensional density distribution, ρ(η, ξ), in a meridional plane.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

It is found, via our model calculation, that α Eridani must be highly flattened with eccentricity $\mathcal {E}=0.7222$ or R_e /R_p = 1.45, which is slightly larger than the observational estimate R_e /R_p = 1.4 (Carcifi et al. 2008). Evidently, perturbation theories based on an expansion using a small rotation parameter around spherical geometry based on $\mathcal {E}\ll 1$ are inapplicable to the present problem with $\mathcal {E}= \mbox{O}(1)$.