Challenges Constructing Ellipsoidal-Like Configurations

First, let's review the three different approaches that we have described for constructing Riemann S-type ellipsoids. Then let's see how these relate to the technique that has been used to construct infinitesimally thin, nonaxisymmetric disks.

Riemann S-type Ellipsoids

Usually, the density, <math>~\rho</math>, and the pair of axis ratios, <math>~b/a</math> and <math>~c/a</math>, are specified. Then, the Poisson equation is solved to obtain <math>~\Phi_\mathrm{grav}</math> in terms of <math>~A_1</math>, <math>~A_2</math>, and <math>~A_3</math>. The aim, then, is to determine the value of the central enthalpy, <math>~H_0</math> — alternatively, the thermal energy density, <math>~\Pi</math> — and the two parameters, <math>~\Omega_f</math> and <math>~\lambda</math>, that determine the magnitude of the velocity flow-field. Keep in mind that, as viewed from a frame of reference that is spinning with the ellipsoid (at angular frequency, <math>~\Omega_f</math>), the adopted (rotating-frame) velocity field is,

<math>~\lambda \biggl[ \boldsymbol{\hat\imath} \biggl( \frac{a}{b}\biggr) y - \boldsymbol{\hat\jmath} \biggl( \frac{b}{a} \biggr) x \biggr] \, .</math>

Hence, the inertial-frame velocity is given by the expression,

<math>~\bold{u} + \bold{\hat{e}}_\varphi \Omega_f \varpi \, .</math>

While we will fundamentally rely on the <math>~(\Omega_f, \lambda)</math> parameter pair to define the velocity flow-field, in discussions of Riemann S-type ellipsoids it is customary to also refer to the following two additional parameters: The (rotating-frame) vorticity,

<math>~\boldsymbol{\zeta} \equiv \boldsymbol{\nabla \times}\bold{u}</math>	<math>~=</math>	<math>~ \boldsymbol{\hat\imath} \biggl[ \frac{\partial u_z}{\partial y} - \frac{\partial u_y}{\partial z} \biggr] + \boldsymbol{\hat\jmath} \biggl[ \frac{\partial u_x}{\partial z} - \frac{\partial u_z}{\partial x} \biggr] + \bold{\hat{k}} \biggl[ \frac{\partial u_y}{\partial x} - \frac{\partial u_x}{\partial y} \biggr] </math>
	<math>~=</math>	<math>~\bold{\hat{k}} \biggl[ - \lambda \biggl(\frac{b}{a} + \frac{a}{b}\biggr) \biggr] \, ;</math>

and the dimensionless frequency ratio,

<math>~\equiv</math>

<math>~\frac{ \zeta}{\Omega_f} \, .</math>

2^nd-Order TVE Expressions

As we have discussed in detail in an accompanying chapter, the three diagonal elements of the <math>~(3 \times 3)</math> 2^nd-order tensor virial equation are sufficient to determine the equilibrium values of <math>~\Pi</math>, <math>~\Omega_3</math>, and <math>~\zeta_3</math>.

Indices

2^nd-Order TVE Expressions that are Relevant to Riemann S-Type Ellipsoids

<math>~ \biggl[ \frac{3\cdot 5}{2^2\pi a b c\rho} \biggr] \Pi +\biggl\{ \Omega_3^2 + 2 \biggl[ \frac{b^2}{b^2+a^2}\biggr] \Omega_3 \zeta_3 ~-~(2\pi G\rho) A_1 \biggr\} a^2 + \biggl[ \frac{a^2}{a^2 + b^2}\biggr]^2 \zeta_3^2 b^2 </math>

<math>~ \biggl[ \frac{3\cdot 5}{2^2\pi a b c \rho} \biggr]\Pi + \biggl[ \frac{b^2}{b^2+a^2}\biggr]^2 \zeta_3^2 a^2 + \biggl\{ \Omega_3^2 + 2 \biggl[ \frac{a^2}{a^2 + b^2}\biggr] \Omega_3 \zeta_3 ~-~( 2\pi G \rho) A_2 \biggr\}b^2 </math>

<math>~ \biggl[ \frac{3\cdot 5}{2^2\pi abc\rho} \biggr]\Pi - (2\pi G \rho)A_3 c^2 </math>

The <math>~(i, j) = (3, 3)</math> element gives <math>~\Pi</math> directly in terms of known parameters. The <math>~(1, 1)</math> and <math>~(2, 2)</math> elements can then be combined in a couple of different ways to obtain a coupled set of expressions that define <math>~\Omega_3</math> and <math>~f \equiv \zeta_3/\Omega_3</math>.

<math>~\biggl[ \frac{b^2 a^2}{b^2+a^2}\biggr] f \Omega_3^2 </math>	<math>~=</math>	<math>~ \pi G\rho \biggl[ \frac{(A_1 - A_2)a^2b^2}{ b^2 - a^2} - A_3 c^2\biggr] \, ; </math>
[ EFE, Chapter 7, §48, Eq. (34) ]

and,

<math>~ \Omega_3^2 \biggl\{1 + \biggl[ \frac{a^2b^2}{(a^2 + b^2)^2}\biggr] f^2 \biggr\} </math>	<math>~=</math>	<math>~ \frac{2\pi G\rho}{ (a^2-b^2) } \biggl[ A_1 a^2 - A_2 b^2 \biggr] \, . </math>
[ EFE, Chapter 7, §48, Eq. (33) ]

Ou's (2006) Detailed Force Balance

In a separate accompanying chapter, we have described in detail how Ou(2006) used, essentially, the HSCF technique to solve the detailed force-balance equations. Beginning with the,

Eulerian Representation
of the Euler Equation
as viewed from a Rotating Reference Frame

<math>\biggl[\frac{\partial\vec{v}}{\partial t}\biggr]_{rot} + ({\vec{v}}_{rot}\cdot \nabla) {\vec{v}}_{rot}= - \frac{1}{\rho} \nabla P - \nabla \Phi_\mathrm{grav}

- {\vec{\Omega}}_f \times ({\vec{\Omega}}_f \times \vec{x}) - 2{\vec{\Omega}}_f \times {\vec{v}}_{rot} \, ,</math>

it can be shown that, for the velocity fields associated with all Riemann S-type ellipsoids,

<math>~({\vec{v}}_{rot}\cdot \nabla) {\vec{v}}_{rot}</math>	<math>~=</math>	<math>~ -\nabla \biggl[ \frac{1}{2} \lambda^2(x^2 + y^2) \biggr] \, ; </math>
<math>~- {\vec{\Omega}}_f \times ({\vec{\Omega}}_f \times \vec{x})</math>	<math>~=</math>	<math>~ +\nabla\biggl[\frac{1}{2} \Omega_f^2 (x^2 + y^2) \biggr] \, ; </math>
<math>~- 2{\vec{\Omega}}_f \times {\vec{v}}_{rot} </math>	<math>~=</math>	<math>~ - \nabla\biggl[ \Omega_f \lambda\biggl( \frac{b}{a} x^2 + \frac{a}{b}y^2 \biggr) \biggr] \, . </math>

Hence, within each steady-state configuration the following Bernoulli's function must be uniform in space:

<math>~ H + \Phi_\mathrm{grav} - \frac{1}{2} \Omega_f^2(x^2 + y^2) - \frac{1}{2} \lambda^2(x^2 + y^2) + \Omega_f \lambda \biggl(\frac{b}{a}x^2 + \frac{a}{b}y^2 \biggr) </math>

Ou(2006), p. 550, §2, Eq. (6)

where <math>~C_B</math> is a constant. So, at the surface of the ellipsoid (where the enthalpy H = 0) on each of its three principal axes, the equilibrium conditions demanded by the expression for detailed force balance become, respectively:

On the x-axis, where (x, y, z) = (a, 0, 0):

<math>~2\biggl[ \frac{C_B}{a^2} + (\pi G\rho)I_\mathrm{BT} \biggr]</math>

<math>~ (2\pi G \rho) A_1 - \Omega_f^2 - \lambda^2 + 2\Omega_f \lambda \biggl(\frac{b}{a} \biggr) </math>

On the y-axis, where (x, y, z) = (0, b, 0):

<math>~2\biggl[ \frac{C_B}{a^2} + (\pi G\rho)I_\mathrm{BT} \biggr]</math>

<math>~ (2\pi G \rho) A_2 \biggl( \frac{b^2}{a^2}\biggr) - \Omega_f^2 \biggl( \frac{b^2}{a^2} \biggr) - \lambda^2\biggl( \frac{b^2}{a^2} \biggr) + 2\Omega_f \lambda \biggl(\frac{b}{a}\biggr) </math>

On the z-axis, where (x, y, z) = (0, 0, c):

<math>~\Rightarrow ~~~ 2 \biggl[ \frac{C_B}{a^2} + (\pi G\rho)I_\mathrm{BT}\biggr]</math>

<math>~=</math>

<math>~ (2\pi G \rho) A_3 \biggl( \frac{c^2}{a^2}\biggr) </math>

This third expression can be used to replace the left-hand-side of the first and second expressions. Then via some additional algebraic manipulation, the first and second expressions can be combined to provide the desired solutions for the parameter pair, <math>~(\Omega_f, \lambda)</math>, namely,

<math>~\frac{\Omega_f^2}{(\pi G \rho)}</math>

<math>~\frac{1}{2} \biggl[M + \sqrt{ M^2 - 4N^2} \biggr] \, ,</math>

and

<math>~\frac{\lambda^2}{(\pi G \rho)}</math>

<math>~\frac{1}{2} \biggl[M - \sqrt{ M^2 - 4N^2} \biggr] \, ,</math>

Ou(2006), p. 551, §2, Eqs. (15) & (16)

where,

<math>~M</math>	<math>~\equiv</math>	<math>~ 2\biggl[ A_1 - A_2 \biggl( \frac{b^2}{a^2}\biggr) \biggr]\biggl[ \frac{a^2}{a^2 - b^2} \biggr] \, ,</math> and,
<math>~N</math>	<math>~\equiv</math>	<math>~ \frac{1}{a b ( a^2 - b^2 )}\biggl[ A_3 ( a^2 - b^2 )c^2 - (A_2 - A_1) a^2 b^2 \biggr] \, . </math>

Hybrid Scheme

In a separate chapter we have detailed the hybrid scheme. For steady-state configurations, the three components of the combined Euler + Continuity equations give,

Hybrid Scheme Summary for Steady-State Configurations

<math>~\boldsymbol{\hat{k}:}</math>	<math>~ \bold\nabla \cdot (\rho v_z \bold{u}) </math>	<math>~=</math>	<math>~\bold{\hat{k}} \cdot (\rho \bold{a}) \, ;</math>
<math>~\bold{\hat{e}_\varpi:}</math>	<math>~ \bold\nabla \cdot (\rho v_\varpi \bold{u}) </math>	<math>~=</math>	<math>~\bold{\hat{e}}_\varpi \cdot (\rho \bold{a}) + \frac{v_\varphi^2}{\varpi} \, ;</math>
<math>~\bold{\hat{e}_\varphi:}</math>	<math>~ \bold\nabla \cdot (\rho \varpi v_\varphi \bold{u}) </math>	<math>~=</math>	<math>~\bold{\hat{e}}_\varphi \cdot (\rho \varpi \bold{a}) \, .</math>

In this context, the vector acceleration that drives the fluid flow is, simply,

<math>~-\nabla(H + \Phi_\mathrm{grav} ) \, .</math>

Then, for the velocity flow-patterns in Riemann S-type ellipsoids, we have,

<math>~\nabla \cdot (\rho v_z \bold{u})</math>	<math>~=</math>	<math>~0</math> (because <math>~v_z = 0</math>);
<math>~\nabla \cdot (\rho v_\varpi \bold{u})</math>	<math>~=</math>	<math>~\frac{\lambda^2}{\varpi^3} \biggl[\frac{a}{b} - \frac{b}{a} \biggr] \biggl\{ y^4 \biggl(\frac{a}{b}\biggr) - x^4 \biggl(\frac{b}{a}\biggr) \biggr\}\rho \, ; </math>
<math>~\nabla \cdot (\rho \varpi v_\varphi \bold{u})</math>	<math>~=</math>	<math>~ 2 \lambda xy \Omega_f \biggl[\frac{a}{b} - \frac{b}{a} \biggr]\rho \, ; </math>
<math>~\varpi v_\varphi</math>	<math>~=</math>	<math>~ - \biggl[ \lambda \biggl(\frac{b}{a}\biggr) - \Omega_f\biggr]x^2 - \biggl[ \lambda \biggl(\frac{a}{b}\biggr) - \Omega_f\biggr]y^2 \, . </math>

Vertical Component: Given that <math>~\bold{\hat{k}}\cdot (\rho \bold{a}) = 0</math>, we deduce that,

Azimuthal Component: Irrespective of the <math>~(x, y, z)</math> location of each fluid element, this component requires,

<math>~ - a b \lambda \Omega_f </math>

<math>~ \pi G \rho \biggl[ \frac{( A_1 - A_2 )a^2b^2}{b^2 - a^2} - c^2 A_3 \biggr] \, . </math>

Radial Component: After inserting the "azimuthal component" relation and marching through a fair amount of algebraic manipulation, we find that Irrespective of the <math>~(x, y, z)</math> location of each fluid element, this component requires,

<math>~ \frac{2\pi G \rho }{(a^2 - b^2) } \biggl[ A_1 a^2 - A_2 b^2 \biggr] </math>

<math>~ \biggl[ \lambda^2 + \Omega_f^2\biggr] \, . </math>

Compressible Structures

Here we draw heavily on the published work of Korycansky & Papaloizou (1996, ApJS, 105, 181; hereafter KP96) that we have reviewed in a separate chapter, and on the doctoral dissertation of Saied W. Andalib (1998).

Part I

Returning to the above-mentioned,

Eulerian Representation
of the Euler Equation
as viewed from a Rotating Reference Frame

<math>\frac{\partial \bold{u}}{\partial t} + (\bold{u}\cdot \nabla) \bold{u} = - \frac{1}{\rho} \nabla P - \nabla \Phi_\mathrm{grav}

- {\vec{\Omega}}_f \times ({\vec{\Omega}}_f \times \vec{x}) - 2{\vec{\Omega}}_f \times \bold{u} \, ,</math>

we next note — as we have done in our broader discussion of the Euler equation — that,

<math> (\bold{u} \cdot\nabla)\bold{u} = \frac{1}{2}\nabla(\bold{u} \cdot \bold{u}) - \bold{u} \times(\nabla\times\bold{u}) = \frac{1}{2}\nabla(u^2) + \boldsymbol\zeta \times \bold{u} , </math>

where, as above, <math>\boldsymbol\zeta \equiv \nabla\times \bold{u}</math> is the vorticity. Making this substitution, we obtain what is essentially equation (7) of KP96, that is, the,

Euler Equation
written in terms of the Vorticity and
as viewed from a Rotating Reference Frame

<math>\frac{\partial \bold{u}}{\partial t} + (\boldsymbol\zeta+2{\vec\Omega}_f) \times {\bold{u}}= - \frac{1}{\rho} \nabla P - \nabla \biggl[\Phi + \frac{1}{2}u^2 - \frac{1}{2}|{\vec{\Omega}}_f \times \vec{x}|^2 \biggr]</math> .

Hence, in steady-state, the Euler equation becomes,

<math> \nabla F_B + \vec{A} = 0 , </math>

where, the scalar "Bernoulli" function,

<math> F_B \equiv \frac{1}{2}u^2 + H + \Phi - \frac{1}{2}|\Omega\hat{k} \times \vec{x}|^2 ; </math>

and,

<math> \vec{A} \equiv ({\boldsymbol\zeta}+2{\vec\Omega}_f) \times {\bold{u}} . </math>

For later use …

Curl of steady-state Euler equation:

<math>~0</math>	<math>~=</math>	<math>~\nabla F_B + \bold{A}</math>
<math>~\Rightarrow~~~0</math>	<math>~=</math>	<math>~\nabla \times \biggl[ \nabla F_B + \bold{A} \biggr]</math>
<math>~\Rightarrow~~~0</math>	<math>~=</math>	<math>~\nabla \times \bold{A} \, .</math>

This last step is justified because the curl of any gradient is zero.

KP96 only deal with two-dimensional motion in the equatorial plane and, hence, there is no vertical motion:
Hence, <math>~\bold{u}</math> lies in the equatorial plane; both <math>~\vec\zeta</math> and <math>~\vec\Omega_f</math> only have z-components; and, <math>~\bold{A}</math> is perpendicular to both <math>~\vec\Omega_f</math> and <math>~\bold{u}</math>. Also, given that <math>~\bold{A}</math> necessarily lies in the equatorial plane, its curl will only have a z-component, that is,

<math>~\nabla \times \bold{A} = 0</math>

<math>~\Leftrightarrow</math>

<math>~[\nabla \times \bold{A}]_z = 0 \, .</math>

"Dot" <math>~\bold{u}</math> into the steady-state Euler equation:

<math>~0</math>	<math>~=</math>	<math>~\nabla F_B + \bold{A}</math>
<math>~\Rightarrow~~~0</math>	<math>~=</math>	<math>~\bold{u} \cdot \biggl[ \nabla F_B + \bold{A} \biggr]</math>
<math>~\Rightarrow~~~0</math>	<math>~=</math>	<math>~\bold{u} \cdot \nabla F_B \, .</math>

This last step is justified because <math>~\bold{A}</math> is necessarily always perpendicular to <math>~\bold{u}</math>.

We will leave discussion of the Euler equation, for the moment, and instead look at the continuity equation. As viewed from the rotating frame of reference,

<math>~\frac{\partial (\rho \bold{u})}{\partial t} + \nabla\cdot (\rho\bold{u})</math>

If we are able to write the momentum density (in the rotating frame) in terms of a stream-function, <math>~\Psi</math>, such that,

<math>~ \nabla \times (\boldsymbol{\hat{k}} \Psi) = \boldsymbol{\hat\imath} \biggl[ \frac{\partial \Psi}{\partial y} \biggr] - \boldsymbol{\hat\jmath} \biggl[ \frac{\partial \Psi}{\partial x}\biggr] \, ,</math>

then satisfying the steady-state continuity equation is guaranteed because the divergence of a curl is always zero. Note, as well, that when written in terms of this stream-function, the z-component of the vorticity will be,

<math>~\zeta_z </math>	<math>~=</math>	<math>~ \frac{\partial u_y}{\partial x} - \frac{\partial u_x}{\partial y} </math>
	<math>~=</math>	<math>~ \frac{\partial }{\partial x}\biggl[- \frac{1}{\rho} \frac{\partial \Psi}{\partial x} \biggr] - \frac{\partial }{\partial y} \biggl[\frac{1}{\rho} \frac{\partial \Psi}{\partial y} \biggr] \, . </math>

Note that the steady-state continuity equation may be rewritten in the form,

<math>~\nabla\cdot \bold{u}</math>

<math>~ - \bold{u} \cdot \nabla [ \ln \rho] \, . </math>

It can also be shown that the condition, <math>~[\nabla \times \bold{A}]_z = 0</math> can be rewritten as,

<math>~\nabla\cdot \bold{u}</math>

<math>~ - \bold{u} \cdot \nabla [ \ln(2\Omega_f + \zeta_z] \, . </math>

By combining these last two expressions, we appreciate that,

<math>~ \bold{u} \cdot \nabla \ln \biggl[ \frac{(2\Omega_f + \zeta_z)}{\rho} \biggr] </math>

This means that, in the steady-state flow whose spatial structure we are seeking, the velocity vector, <math>~\bold{u}</math> (and also the momentum density vector, <math>~\rho \bold{u}</math>) must everywhere be tangent to contours of constant vortensity, <math>~[(2\Omega_f + \zeta_z)/\rho]</math>.

We need a function <math>~g(\Psi) </math> such that,

<math>~g(\Psi) </math>	<math>~=</math>	<math>~ \frac{1}{\rho} \biggl\{ \zeta_z + 2\Omega_f \biggr\} </math>
	<math>~=</math>	<math>~ \frac{1}{\rho} \biggl\{ 2\Omega_f - \frac{\partial }{\partial x}\biggl[\frac{1}{\rho} \frac{\partial \Psi}{\partial x} \biggr] - \frac{\partial }{\partial y} \biggl[\frac{1}{\rho} \frac{\partial \Psi}{\partial y} \biggr] \biggr\} \, . </math>

Let's try, <math>~\Psi = \rho^2</math>, and

<math>~\rho</math>	<math>~=</math>	<math>~\rho_c \biggl\{1 - \biggl[ \frac{y^2}{b^2} + \frac{x^2}{a^2}\biggr]\biggr\} </math>
<math>~\Rightarrow~~~\frac{\partial^2 \rho}{\partial x^2} = -\frac{\partial}{\partial x}\biggl\{ \frac{2\rho_c x}{a^2}\biggr\} = - \frac{2\rho_c}{a^2}</math>	and,	<math>~\frac{\partial^2 \rho}{\partial y^2} = - \frac{\partial}{\partial y}\biggl\{ \frac{2\rho_c y}{b^2} \biggr\} = - \frac{2\rho_c}{b^2} \, .</math>

Then,

<math>~g(\Psi) </math>	<math>~=</math>	<math>~ \frac{1}{\rho} \biggl\{ 2\Omega_f - \frac{\partial }{\partial x}\biggl[\frac{1}{\rho} \frac{\partial \rho^2}{\partial x} \biggr] - \frac{\partial }{\partial y} \biggl[\frac{1}{\rho} \frac{\partial \rho^2}{\partial y} \biggr] \biggr\} </math>
	<math>~=</math>	<math>~ \frac{1}{\rho} \biggl\{ 2\Omega_f - 2 \biggl[ \frac{\partial^2 \rho}{\partial x^2} \biggr] - 2 \biggl[\frac{\partial^2 \rho}{\partial y^2} \biggr] \biggr\} </math>
	<math>~=</math>	<math>~ \frac{1}{\rho} \biggl\{ 2\Omega_f + 4\rho_c \biggl[ \frac{1}{a^2} + \frac{1}{b^2} \biggr] \biggr\} \, . </math>

Hence,

<math>~ \biggl\{ 2\Omega_f + 4\rho_c \biggl[ \frac{1}{a^2} + \frac{1}{b^2} \biggr] \biggr\} \Psi^{-1 / 2} \, . </math>

Next, given that,

we conclude that, to within an additive constant,

<math>~ F_B(\Psi)</math>	<math>~=</math>	<math>~ - \int g(\Psi) d\Psi = - g_0 \int \Psi^{-1 / 2} d\Psi </math>
	<math>~=</math>	<math>~ - 2g_0 \Psi^{1 / 2} </math>
	<math>~=</math>	<math>~ - 2g_0 \rho \, , </math>

where,

<math>~\equiv</math>

<math>~ \biggl\{ 2\Omega_f + 4\rho_c \biggl[ \frac{1}{a^2} + \frac{1}{b^2} \biggr] \biggr\} \, . </math>

Here's what to do:

Given <math>~g(\Psi)</math>, write out the functional forms of <math>~\bold{A}</math> and <math>~F_B(\Psi)</math>. Then see if <math>~\nabla F_B = - \bold{A}</math>.

Part II

Consider a steady-state configuration that is the compressible analog of a Riemann S-type ellipsoid; even better, give the configuration a "peanut" shape rather than the shape of an ellipsoid. As viewed from a frame that is spinning with the configuration's overall angular velocity, <math>~\vec\Omega_f = \boldsymbol{\hat{k}} \Omega_f</math>, generally we expect the configuration's internal (and surface) flow to be represented by a set of nested stream-lines and at every <math>~(x, y)</math> location the fluid's velocity (and its momentum-density vector) will be tangent to the stream-line that runs through that point. It is customary to represent the stream-function by a scalar quantity, <math>~\Psi(x, y)</math>, appreciating that each stream-line will be defined by a curve, <math>~\Psi = \mathrm{constant}</math>; also, the local spatial gradient of <math>~\Psi(x,y)</math> will be perpendicular to the local stream-line, hence it will be perpendicular to the local velocity vector as well. If we specifically introduce the stream-function via the relation,

<math>~ \nabla \times (\boldsymbol{\hat{k}} \Psi) = \boldsymbol{\hat\imath} \biggl[ \frac{\partial \Psi}{\partial y} \biggr] - \boldsymbol{\hat\jmath} \biggl[ \frac{\partial \Psi}{\partial x}\biggr] \, ,</math>

it will display all of the just-described attributes and we are also guaranteed that the steady-state continuity equation will be satisfied everywhere, because the divergence of a curl is always zero.

We also have demonstrated that the vector, <math>~\bold{A}</math>, has the right properties if,

<math>~ \bold{u} \cdot \nabla \ln \biggl[ \frac{(2\Omega_f + \zeta_z)}{\rho} \biggr] </math>

This means that, at every location in the plane of the fluid flow, the gradient of the vortensity also must be perpendicular to the velocity vector. This constraint can be immediately satisfied if we simply demand that the vortensity be a function of the stream-function, <math>~\Psi</math>, that is, we need,

<math>~\frac{(2\Omega_f + \zeta_z)}{\rho}</math>

In other words, the scalar vortensity is constant along each stream-line. And, once we have determined the mathematical expression for this function, we will know that,

<math>~[ \boldsymbol{\hat{k}} g(\Psi) ] \times \rho\bold{u} \, ;</math>

Furthermore, we should be able to determine the mathematical expression for <math>~F_B(x,y)</math> because,

As a check, we should find that,

<math>~\nabla F_B + \bold{A}</math>

Part III

Here we closely follow Chapter 4 of Saied W. Andalib (1998).

From §4.1 (p. 80): "Euler's equation, the equation of continuity, the Poisson equation and the equation of state … govern the dynamics and evolution of these equilibrium configurations."

Equation of Continuity

In steady state, <math>~\partial (\rho\bold{u})/\partial t = 0</math>. Hence the rotating-frame-based continuity equation becomes,

<math>~\nabla \cdot (\rho\bold{u})</math>

If we insist that the momentum-density vector be expressible in terms of the curl of a vector — for example,

<math>~\rho\bold{u}</math>	<math>~=</math>	<math>~\nabla \times \boldsymbol{\mathfrak{J}} \, ,</math>
Saied W. Andalib (1998), §4.1, p. 80, Eq. (4.1)

then satisfying this steady-state continuity equation is guaranteed because the divergence of a curl is always zero. "The task of satisfying the steady-state equation of continuity then shifts to identifying an appropriate expression for the vector potential, <math>~\boldsymbol{\mathfrak{J}} \, .</math>" In the most general case, in terms of this vector potential the three Cartesian components of the momentum-density vector are,

<math>~\rho u_x</math>	<math>~=</math>	<math>~ \frac{\partial \mathfrak{J}_z}{\partial y} - \frac{\partial \mathfrak{J}_y}{\partial z} \, ; </math>
<math>~\rho u_y</math>	<math>~=</math>	<math>~ \frac{\partial \mathfrak{J}_x}{\partial z} - \frac{\partial \mathfrak{J}_z}{\partial x} \, ; </math>
<math>~\rho u_z</math>	<math>~=</math>	<math>~ \frac{\partial \mathfrak{J}_y}{\partial x} - \frac{\partial \mathfrak{J}_x}{\partial y} \, . </math>

Here, we will follow Andalib's lead and only look for fluid flows with no vertical motions. That is to say, we will set <math>~\rho u_z = 0</math>, in which case this last expression establishes the constraint,

<math>~\frac{\partial \mathfrak{J}_y}{\partial x}</math>

<math>~ \frac{\partial \mathfrak{J}_x}{\partial y} \, . </math>

"A general solution to this equation can be found only if there exists a scalar function <math>~\Gamma(x, y, z)</math> such that …"

<math>~\mathfrak{J}_y = \frac{\partial \Gamma}{\partial y}</math>

and,

<math>~\mathfrak{J}_x = \frac{\partial \Gamma}{\partial x} \, ;</math>

note that this adopted functional behavior works because the constraint becomes,

<math>~\frac{\partial^2 \Gamma}{\partial x \partial y}</math>

<math>~\frac{\partial^2 \Gamma}{\partial y \partial x} \, .</math>

Hence, the expressions for the x- and y-components of the momentum-density vector may be rewritten, respectively, as,

<math>~\rho u_x</math>	<math>~=</math>	<math>~ \frac{\partial \mathfrak{J}_z}{\partial y} - \frac{\partial}{\partial z} \biggl[ \frac{\partial \Gamma}{\partial y} \biggr] = + \frac{\partial}{\partial y}\biggl[ \mathfrak{J}_z - \frac{\partial\Gamma}{\partial z} \biggr] \, ; </math>
<math>~\rho u_y</math>	<math>~=</math>	<math>~ \frac{\partial }{\partial z}\biggl[ \frac{\partial \Gamma}{\partial x} \biggr] - \frac{\partial \mathfrak{J}_z}{\partial x} = - \frac{\partial}{\partial x}\biggl[ \mathfrak{J}_z - \frac{\partial\Gamma}{\partial z} \biggr] \, . </math>

If we again follow Andalib's lead and only look for models in which the x-y-plane flow is independent of the vertical coordinate, z, then, <math>~\mathfrak{J}_z</math> and <math>~\partial\Gamma/\partial z</math> must be functions of x and y only. Therefore, <math>~\mathfrak{J}_z</math> is independent of z and <math>~\Gamma</math> is at most linear in z. Now, rather than focusing on the determination of <math>~\Gamma(x,y)</math>, we can just as well define the scalar function,

<math>~\equiv</math>

<math>~ \mathfrak{J}_z - \frac{\partial\Gamma}{\partial z} \, , </math>

in which case "… the components of the momentum density may be written as:"

<math>~ \boldsymbol{\hat\imath} \frac{\partial \Psi}{\partial y} - \boldsymbol{\hat\jmath} \frac{\partial \Psi}{\partial x} \, . </math>

It is straightforward to demonstrate that this expression for the momentum-density vector does satisfy the steady-state continuity equation. "The function <math>~\Psi(x, y)</math> will serve a similar role as the velocity potential for incompressible fluids."

Related Useful Expressions

Given that, by our design, the fluid motion will be confined to the x-y-plane, the fluid vorticity will have only a z-component; that is,

<math>~\vec\zeta = \nabla \times \bold{u} = \boldsymbol{\hat{k}} \zeta_z</math>,

where,

<math>~ \frac{\partial u_y}{\partial x} - \frac{\partial u_x}{\partial y} \, . </math>

And when it is written in terms of <math>~\Psi(x,y)</math>, this z-component of the vorticity will be obtained from the expression,

<math>~ \frac{\partial }{\partial x}\biggl[- \frac{1}{\rho} \frac{\partial \Psi}{\partial x} \biggr] - \frac{\partial }{\partial y} \biggl[\frac{1}{\rho} \frac{\partial \Psi}{\partial y} \biggr] \, . </math>

This is useful to know because, in the Euler equation (see immediately below) we will encounter a term that involves the cross product of the vector, <math>~(\vec\zeta + 2\vec\Omega) </math>, with the rotating-frame-based velocity vector. Appreciating as well that the vector, <math>~\vec\Omega = \boldsymbol{\hat{k}} \Omega_f</math>, only has a nonzero z-component, we recognize that this term may be written as,

<math>~(\vec\zeta + 2\vec\Omega) \times \bold{u}</math>	<math>~=</math>	<math>~\boldsymbol{\hat{k}} \biggl[ \frac{ (\zeta_z + 2\Omega_f)}{\rho} \biggr] \times \biggl[ \boldsymbol{\hat\imath} \frac{\partial \Psi}{\partial y} - \boldsymbol{\hat\jmath} \frac{\partial \Psi}{\partial x} \biggr]</math>
	<math>~=</math>	<math>~\biggl[ \frac{ (\zeta_z + 2\Omega_f)}{\rho} \biggr] \biggl[ \boldsymbol{\hat\jmath} \frac{\partial \Psi}{\partial y} + \boldsymbol{\hat\imath} \frac{\partial \Psi}{\partial x} \biggr]</math>
	<math>~=</math>	<math>~\biggl[ \frac{ (\zeta_z + 2\Omega_f)}{\rho} \biggr] \nabla\Psi \, .</math>
Saied W. Andalib (1998), §4.2, p. 83, Eq. (4.13)

Euler Equation

We begin with the,

Euler Equation
written in terms of the Vorticity and
as viewed from a Rotating Reference Frame

<math>\frac{\partial \bold{u}}{\partial t} + (\boldsymbol\zeta+2{\vec\Omega}_f) \times {\bold{u}}= - \frac{1}{\rho} \nabla P - \nabla \biggl[\Phi + \frac{1}{2}u^2 - \frac{1}{2}|{\vec{\Omega}}_f \times \vec{x}|^2 \biggr]</math> .

Next, we rewrite this expression to incorporate the following three realizations:

For a barotropic fluid, the term involving the pressure gradient can be replaced with a term involving the enthalpy via the relation, <math>~\nabla H = \nabla P/\rho</math>.
The expression for the centrifugal potential can be rewritten as, <math>~\tfrac{1}{2}|\vec\Omega_f \times \vec{x}|^2 = \tfrac{1}{2}\Omega_f^2 (x^2 + y^2)</math>.
In steady state, <math>~\partial \bold{u}/\partial t = 0</math>.

This means that,

<math>~ (\boldsymbol\zeta+2{\vec\Omega}_f) \times {\bold{u}}</math>

<math>- \nabla \biggl[H + \Phi_\mathrm{grav} + \frac{1}{2}u^2 - \frac{1}{2}\Omega_f^2 (x^2 + y^2) \biggr] \, .</math>

If the term on the left-hand-side of this equation can be expressed in terms of the gradient of a scalar function, then it can be readily grouped with all the other terms on the right-hand-side, which already are in the gradient form.

Striving for Gradient Form

As we have already demonstrated above, the term on the left-hand-side of the Euler equation can be rewritten as,

<math>~(\vec\zeta + 2\vec\Omega) \times \bold{u}</math>

<math>~\biggl[ \frac{ (\zeta_z + 2\Omega_f)}{\rho} \biggr] \nabla\Psi \, .</math>

If the term inside the square brackets on the right-hand-side were a constant — that is, independent of position — then it could immediately be moved inside the gradient operator and we will have accomplished our objective. But, while <math>~\Omega_f</math> is constant "… generally the vorticity and density are both functions of <math>~x</math> and <math>~y</math>." As Andalib has explained, "The expression … can be cast in the form of a gradient only if

<math>~\biggl[ \frac{ (\zeta_z + 2\Omega_f)}{\rho} \biggr]</math>

where <math>~g</math> is an arbitrary function." Specifically in this case, the term on the left-hand-side of the Euler equation may be written as,

<math>~(\vec\zeta + 2\vec\Omega) \times \bold{u}</math>	<math>~=</math>	<math>~\nabla F_B(\Psi) </math>
	<math>~=</math>	<math>~\biggl[ \frac{dF_B(\Psi)}{d\Psi} \biggr] \nabla \Psi \, .</math>

That is, we accomplish our objective by recognizing that the sought-after function, <math>~F_B(\Psi)</math>, is obtained from <math>~g(\Psi)</math> via the relation,

For example, try …

<math>~F_B(\Psi)</math>	<math>~=</math>	<math>~C_0 \Psi + \frac{1}{2} C_1 \Psi^2 </math>
Saied W. Andalib (1998), §4.3, p. 85, Eq. (4.22)
<math>~\Rightarrow ~~~ g(\Psi)</math>	<math>~=</math>	<math>~C_0 + C_1 \Psi \, .</math>

This means that,

<math>~\biggl[ \frac{ (\zeta_z + 2\Omega_f)}{\rho} \biggr]</math>	<math>~=</math>	<math>~~C_0 + C_1 \Psi</math>
<math>~\Rightarrow ~~~ 2\Omega_f - C_0 \rho </math>	<math>~=</math>	<math>~~C_1 \rho \Psi - \zeta_z </math>
	<math>~=</math>	<math>~~C_1 \rho \Psi + \frac{\partial }{\partial x}\biggl[\frac{1}{\rho} \frac{\partial \Psi}{\partial x} \biggr] + \frac{\partial }{\partial y} \biggl[\frac{1}{\rho} \frac{\partial \Psi}{\partial y} \biggr] \, .</math>
Saied W. Andalib (1998), §4.3, p. 85, Eq. (4.24)

Having transformed the left-hand-side term into the gradient of the scalar function, <math>~F_B(\Psi)</math>, the Euler equation can now be written as,

<math>\nabla \biggl[H + \Phi_\mathrm{grav} + F_B(\Psi) + \frac{1}{2}u^2 - \frac{1}{2}\Omega_f^2 (x^2 + y^2) \biggr] </math>	<math>~=</math>	<math>~ 0</math>
<math>~\Rightarrow~~~ H + \Phi_\mathrm{grav} + F_B(\Psi) + \frac{1}{2}u^2 - \frac{1}{2}\Omega_f^2 (x^2 + y^2) </math>	<math>~=</math>	<math>C_B \, , </math>

where we will refer to <math>~C_B</math> as the Bernoulli constant.

Strategy

STEP 0: Choose the pair of model-sequence parameters, <math>~(C_0, C_1)</math>, that are associated with the function, <math>~F_B(\Psi)</math>. Hold these fixed during iterations.

STEP 1: Guess a density distribution, <math>~\rho(x,y)</math>. For example, pick the equatorial-plane (uniform) density distribution of a Riemann S-type ellipsoid with an equatorial-plane axis-ratio, <math>~b/a</math> and meridional-plane axis-ratio, <math>~c/a</math>; use the same <math>~b/a</math> ratio to define two points on the configuration's surface throughout the iteration cycle.

STEP 2: Given <math>~\rho(x,y)</math>, solve the Poisson equation to obtain, <math>~\Phi_\mathrm{grav}(x,y)</math>. In the first iteration, this should exactly match the <math>~A_1, A_2, A_3</math> values associated with the chosen Riemann S-type ellipsoid.

STEP 3: Guess a value of <math>~\Omega_f</math> — perhaps the spin-frequency associated with your "initial guess" Riemann ellipsoid — then solve the following two-dimensional, elliptic PDE to obtain <math>~\Psi(x,y)</math> …

<math>~ 2\Omega_f - C_0 \rho </math>

<math>~~C_1 \rho \Psi + \frac{\partial }{\partial x}\biggl[\frac{1}{\rho} \frac{\partial \Psi}{\partial x} \biggr] + \frac{\partial }{\partial y} \biggl[\frac{1}{\rho} \frac{\partial \Psi}{\partial y} \biggr] \, .</math>

Boundary Condition

Moving along various rays from the center of the configuration, outward, the surface is determined by the location along each ray where <math>~H(x,y)</math> goes to zero for the first time. We set <math>~\Psi = 0</math> at these various surface locations. At each of these locations, the velocity vector must be tangent to the surface. This requirement also, then, sets the value of <math>~\partial \Psi/\partial y</math> and <math>~\partial \Psi/\partial x</math> at each location.

STEP 4: Determine (rotating-frame) velocity from knowledge of <math>~\Psi(x,y)</math> and <math>~\rho(x,y)</math>.

<math>~\bold{u}</math>	<math>~=</math>	<math>~ \frac{1}{\rho}\biggl\{ \boldsymbol{\hat\imath} \biggl[ \frac{\partial \Psi}{\partial y} \biggr] - \boldsymbol{\hat\jmath} \biggl[ \frac{\partial \Psi}{\partial x}\biggr] \biggr\} </math>
<math>~\Rightarrow~~~ u^2 = \bold{u} \cdot \bold{u}</math>	<math>~=</math>	<math>~ \frac{1}{\rho^2}\biggl\{ \biggl[ \frac{\partial \Psi}{\partial y} \biggr]^2 + \biggl[ \frac{\partial \Psi}{\partial x}\biggr]^2 \biggr\} \, .</math>

STEP 5: Using the "scalar Euler equation,"

<math>~ H + \Phi_\mathrm{grav} + C_0 \Psi + \frac{1}{2} C_1 \Psi^2 + \frac{1}{2}u^2 - \frac{1}{2}\Omega_f^2 (x^2 + y^2) </math>	<math>~=</math>	<math>C_B \, , </math>
Saied W. Andalib (1998), §4.3, p. 85, Eq. (4.23)

Set <math>~H = 0</math> at two different points on the surface of the configuration — usually at <math>~(x,y) = (a,0)</math> and <math>~(x,y) = (0,b)</math> — to determine values of the two constants, <math>~\Omega_f^2</math> and <math>~C_B</math>.
At all points inside the configuration, determine <math>~H(x,y)</math>.

STEP 6: Use the barotropic equation of state to determine the "new" mass-density distribution from the knowledge of the enthalpy, <math>~H(x,y)</math>.

Compare

Incompressible

<math>~\bold{u}</math>	<math>~=</math>	<math>~\lambda \biggl[ \boldsymbol{\hat\imath} \biggl( \frac{a}{b}\biggr) y - \boldsymbol{\hat\jmath} \biggl( \frac{b}{a} \biggr) x \biggr] \, .</math>
<math>~\Rightarrow ~~~ u^2 \equiv \bold{u}\cdot \bold{u}</math>	<math>~=</math>	<math>~ \lambda^2 \biggl[ b^2 \biggl(\frac{x^2}{a^2}\biggr) + a^2 \biggl(\frac{y^2}{b^2}\biggr) \biggr] \, . </math>
<math>~\Rightarrow~~~ \boldsymbol{\zeta} \equiv \boldsymbol{\nabla \times}\bold{u}</math>	<math>~=</math>	<math>~ \boldsymbol{\hat\imath} \biggl[ \frac{\partial u_z}{\partial y} - \frac{\partial u_y}{\partial z} \biggr] + \boldsymbol{\hat\jmath} \biggl[ \frac{\partial u_x}{\partial z} - \frac{\partial u_z}{\partial x} \biggr] + \bold{\hat{k}} \biggl[ \frac{\partial u_y}{\partial x} - \frac{\partial u_x}{\partial y} \biggr] </math>
	<math>~=</math>	<math>~\bold{\hat{k}} \biggl[ - \lambda \biggl(\frac{b}{a} + \frac{a}{b}\biggr) \biggr] \, .</math>
<math>~\Rightarrow~~~ \boldsymbol{\zeta} \times \bold{u} </math>	<math>~=</math>	<math>~ - \lambda^2 \biggl(\frac{b}{a} + \frac{a}{b}\biggr) \biggl[ \boldsymbol{\hat\imath} \biggl(\frac{b}{a}\biggr)x + \boldsymbol{\hat\jmath} \biggl(\frac{a}{b}\biggr)y \biggr] </math>
	<math>~=</math>	<math>~- \nabla \biggl\{ \frac{\lambda^2(a^2 + b^2)}{2} \biggl[ \frac{x^2}{a^2} + \frac{y^2}{b^2} \biggr] \biggr\} \, . </math>

Returning to the above-mentioned,

Eulerian Representation
of the Euler Equation
as viewed from a Rotating Reference Frame

<math>\frac{\partial \bold{u}}{\partial t} + (\bold{u}\cdot \nabla) \bold{u} = - \frac{1}{\rho} \nabla P - \nabla \Phi_\mathrm{grav}

- {\vec{\Omega}}_f \times ({\vec{\Omega}}_f \times \vec{x}) - 2{\vec{\Omega}}_f \times \bold{u} \, ,</math>

we next note — as we have done in our broader discussion of the Euler equation — that,

<math> (\bold{u} \cdot\nabla)\bold{u} = \frac{1}{2}\nabla(\bold{u} \cdot \bold{u}) - \bold{u} \times(\nabla\times\bold{u}) = \frac{1}{2}\nabla(u^2) + \boldsymbol\zeta \times \bold{u} \, . </math>

Making this substitution, we obtain

<math>~ \frac{\partial \bold{u}}{\partial t} + \biggl[ \frac{1}{2}\nabla(u^2) + \boldsymbol\zeta \times \bold{u} \biggr] </math>

<math>~ - \frac{1}{\rho} \nabla P - \nabla \Phi_\mathrm{grav} - {\vec{\Omega}}_f \times ({\vec{\Omega}}_f \times \vec{x}) - 2{\vec{\Omega}}_f \times \bold{u} </math>

it can be shown that, for the velocity fields associated with all Riemann S-type ellipsoids,

<math>~({\vec{v}}_{rot}\cdot \nabla) {\vec{v}}_{rot}</math>	<math>~=</math>	<math>~ -\nabla \biggl[ \frac{1}{2} \lambda^2(x^2 + y^2) \biggr] \, ; </math>
<math>~- {\vec{\Omega}}_f \times ({\vec{\Omega}}_f \times \vec{x})</math>	<math>~=</math>	<math>~ +\nabla\biggl[\frac{1}{2} \Omega_f^2 (x^2 + y^2) \biggr] \, ; </math>
<math>~- 2{\vec{\Omega}}_f \times {\vec{v}}_{rot} </math>	<math>~=</math>	<math>~ - \nabla\biggl[ \Omega_f \lambda\biggl( \frac{b}{a} x^2 + \frac{a}{b}y^2 \biggr) \biggr] \, . </math>

As a check,

<math>~(\bold{u} \cdot\nabla)\bold{u}</math>	<math>~=</math>	<math>~\frac{1}{2}\nabla(u^2) + \boldsymbol\zeta \times \bold{u} </math>	<math>~=</math>	<math>~\nabla \biggl\{\frac{\lambda^2}{2} \biggl[ b^2 \biggl(\frac{x^2}{a^2}\biggr) + a^2 \biggl(\frac{y^2}{b^2}\biggr) \biggr] - \frac{\lambda^2(a^2 + b^2)}{2} \biggl[ \frac{x^2}{a^2} + \frac{y^2}{b^2} \biggr] \biggr\} </math>
			<math>~=</math>	<math>~\nabla \frac{\lambda^2}{2} \biggl\{ -(x^2 + y^2) \biggr\} \, . </math>

Yes! This expression matches the one that appears just a few lines earlier in this discussion.

Now, let's switch the order of terms in the steady-state Euler equation to permit an easier comparison with our attempt to develop the compressible models.

<math>~ \cancel{\frac{\partial \bold{u}}{\partial t}} + \biggl[ \frac{1}{2}\nabla(u^2) + \boldsymbol\zeta \times \bold{u} \biggr] </math>	<math>~=</math>	<math>~ - \frac{1}{\rho} \nabla P - \nabla \Phi_\mathrm{grav} - {\vec{\Omega}}_f \times ({\vec{\Omega}}_f \times \vec{x}) - 2{\vec{\Omega}}_f \times \bold{u} </math>
<math>~ \Rightarrow~~~ 0 </math>	<math>~=</math>	<math>~ - \nabla \biggl[ H + \Phi_\mathrm{grav} + \frac{1}{2} u^2 - \frac{1}{2} \Omega_f^2 (x^2 + y^2) \biggr] - (\boldsymbol\zeta + 2{\vec{\Omega}}_f )\times \bold{u} \, . </math>

Compressible

Now in our above discussion of Andalib's work, the steady-state form of the Euler equation was formulated as,

<math>\nabla \biggl[H + \Phi_\mathrm{grav} + F_B(\Psi) + \frac{1}{2}u^2 - \frac{1}{2}\Omega_f^2 (x^2 + y^2) \biggr] </math>

It is easy to appreciate that,

<math>~\nabla F_B(\Psi) </math>

<math>~(\boldsymbol\zeta + 2{\vec{\Omega}}_f )\times \bold{u} \, .</math>

As we have shown, in the case of the incompressible (Riemann S-type ellipsoid) models,

<math>~(\boldsymbol\zeta + 2{\vec{\Omega}}_f )\times \bold{u}</math>

<math>~\nabla \biggl\{ \frac{\lambda^2(a^2 + b^2)}{2} \biggl[ \frac{x^2}{a^2} + \frac{y^2}{b^2} \biggr] + \biggl[ \Omega_f \lambda\biggl( \frac{b}{a} x^2 + \frac{a}{b}y^2 \biggr) \biggr] \biggr\} \, . </math>

If we attempt to directly relate these two expressions, we must acknowledge that,

<math>~F_B(\Psi)</math>	<math>~=</math>	<math>~ \frac{\lambda^2(a^2 + b^2)}{2} \biggl[ \frac{x^2}{a^2} + \frac{y^2}{b^2} \biggr] + \biggl[ \Omega_f \lambda\biggl( \frac{b}{a} x^2 + \frac{a}{b}y^2 \biggr) \biggr] </math>
	<math>~=</math>	<math>~ \frac{x^2}{a^2} \biggl[\frac{1}{2}\lambda^2(a^2 + b^2) + \Omega_f \lambda a b \biggr] + \frac{y^2}{b^2} \biggl[ \frac{1}{2}\lambda^2(a^2 + b^2) + \Omega_f \lambda a b \biggr] \, . </math>

As we have discussed above, Andalib (1998) found that some interesting model sequences could be constructed if he adopted the functional form,

Evidently, the incompressible (uniform-density) Riemann S-type ellisoids can be retrieved from our derived compressible-model formalism if we set, <math>~C_1 = 0</math>, and,

with,

<math>~\biggl[ \frac{1}{2}\lambda^2(a^2 + b^2) + \Omega_f \lambda a b \biggr] </math>

<math>~- \frac{\zeta_z}{2}\biggl( \frac{a^2 b^2}{a^2 + b^2}\biggr)\biggl[ 2\Omega_f - \zeta_z \biggr] \, .</math>

Afterthought: Because we want <math>~\Psi(x,y)</math> to go to zero at the surface, it likely will be better to set,

<math>~\Psi \equiv 1 - \biggl[ \frac{x^2}{a^2} + \frac{y^2}{b^2} \biggr] \, ,</math>

then adjust the sign of <math>~C_0</math> and add a constant (zeroth-order term) to the definition of <math>~F_B(\Psi)</math>.

Trial #1

Restricting our discussion to nonaxisymmetric, thin disks, let's assume <math>~\rho</math> is uniform throughout the configuration and that,

<math>~ \Psi_0 \biggl[1 - \biggl( \frac{x^2}{a^2} + \frac{y^2}{b^2} \biggr) \biggr] \, . </math>

This means that,

<math>~\frac{ \partial \Psi}{\partial x} = - \Psi_0 \biggl( \frac{2x}{a^2}\biggr)</math>

and,

<math>~\frac{ \partial \Psi}{\partial y} = - \Psi_0 \biggl( \frac{2y}{b^2}\biggr) \, .</math>

The momentum density vector is governed by the relation,

<math>~\rho\bold{u}</math>	<math>~=</math>	<math>~ \boldsymbol{\hat\imath} \biggl[ \frac{\partial \Psi}{\partial y} \biggr] - \boldsymbol{\hat\jmath} \biggl[ \frac{\partial \Psi}{\partial x}\biggr] </math>
	<math>~=</math>	<math>~ \Psi_0 \biggl\{- \boldsymbol{\hat\imath} \biggl[ \frac{2y}{b^2} \biggr] + \boldsymbol{\hat\jmath} \biggl[ \frac{2x}{a^2}\biggr] \biggr\} \, . </math>
[<math>~\Psi</math> has units of "density × length² per time"]

First of all, let's see if the steady-state continuity equation is satisfied:

<math>~\nabla \cdot (\rho \bold{u})</math>	<math>~=</math>	<math>~ \frac{\partial (\rho u_x)}{\partial x} + \frac{\partial (\rho u_y)}{\partial y} </math>
	<math>~=</math>	<math>~\rho \biggl\{ \frac{\partial u_x}{\partial x} + \frac{\partial u_y}{\partial y} \biggr\} </math>
	<math>~=</math>	<math>~\Psi_0 \rho \biggl\{ \frac{\partial }{\partial x}\biggl[ - \frac{2y}{b^2} \biggr] + \frac{\partial }{\partial y}\biggl[ \frac{2x}{a^2} \biggr] \biggr\} </math>
	<math>~=</math>	<math>~0 \, .</math> Q.E.D.

Next, let's determine the z-component of the vorticity and the vortensity:

<math>~\zeta_z</math>	<math>~=</math>	<math>~ \frac{\partial u_y}{\partial x} - \frac{\partial u_x}{\partial y} </math>
	<math>~=</math>	<math>~\Psi_0 \biggl\{ \frac{\partial }{\partial x}\biggl[ \frac{2x}{\rho a^2} \biggr] + \frac{\partial }{\partial y}\biggl[ \frac{2y}{\rho b^2} \biggr] \biggr\} </math>
	<math>~=</math>	<math>~ \frac{2\Psi_0}{\rho} \biggl[ \frac{1}{a^2} + \frac{1}{b^2} \biggr] \, ; </math>
<math>~\Rightarrow~~~ \frac{(\zeta_z + 2\Omega_f)}{\rho}</math>	<math>~=</math>	<math>~\frac{1}{\rho} \biggl\{ \frac{2\Psi_0}{\rho} \biggl[ \frac{1}{a^2} + \frac{1}{b^2} \biggr] + 2\Omega_f \biggr\}\, . </math>

This means that,

<math>~ (\vec\zeta + 2\vec\Omega) \times \bold{u} </math>	<math>~=</math>	<math>~\biggl\{ \frac{(\zeta_z + 2\Omega_f)}{\rho} \biggr\} \boldsymbol{\hat{k}} \times (\rho \bold{u} )</math>
	<math>~=</math>	<math>~\Psi_0 \biggl\{ \frac{(\zeta_z + 2\Omega_f)}{\rho} \biggr\} \boldsymbol{\hat{k}} \times \biggl\{- \boldsymbol{\hat\imath} \biggl[ \frac{2y}{b^2} \biggr] + \boldsymbol{\hat\jmath} \biggl[ \frac{2x}{a^2}\biggr] \biggr\}</math>
	<math>~=</math>	<math>~-2 \Psi_0 \biggl[ \frac{(\zeta_z + 2\Omega_f)}{\rho} \biggr] \biggl[ \boldsymbol{\hat\imath} \biggl( \frac{x}{a^2}\biggr) + \boldsymbol{\hat\jmath} \biggl( \frac{y}{b^2} \biggr) \biggr] \, .</math>

But this entire process was designed to ensure that,

<math>~ (\vec\zeta + 2\vec\Omega) \times \bold{u} </math>

<math>~\nabla F_B(\Psi) \, ,</math>

where, <math>~F_B(\Psi) \equiv C_0 \Psi</math>. Let's see if we get the same expression …

<math>~\nabla F_B(\Psi)</math>	<math>~=</math>	<math>~\biggl[ \boldsymbol{\hat\imath} \frac{\partial}{\partial x} + \boldsymbol{\hat\jmath} \frac{\partial}{\partial y} \biggr]C_0 \Psi</math>
	<math>~=</math>	<math>~ C_0 \Psi_0\biggl[ \boldsymbol{\hat\imath} \frac{\partial}{\partial x} + \boldsymbol{\hat\jmath} \frac{\partial}{\partial y} \biggr] \biggl[1 - \biggl( \frac{x^2}{a^2} + \frac{y^2}{b^2} \biggr) \biggr] </math>
	<math>~=</math>	<math>~ - 2C_0 \Psi_0\biggl[ \boldsymbol{\hat\imath} \biggl( \frac{x}{a^2} \biggr) + \boldsymbol{\hat\jmath} \biggl( \frac{y}{b^2} \biggr) \biggr] \, . </math>

This is indeed the same expression as above if we set,

<math>~ \frac{(\zeta_z + 2\Omega_f)}{\rho} \, .</math>

Hooray!!

Finally, let's make sure that the elliptic PDE "constraint" equation is satisfied.

<math>~ 2\Omega_f </math>	<math>~=</math>	<math>~C_0 \rho+~C_1 \rho \Psi + \frac{\partial }{\partial x}\biggl[\frac{1}{\rho} \frac{\partial \Psi}{\partial x} \biggr] + \frac{\partial }{\partial y} \biggl[\frac{1}{\rho} \frac{\partial \Psi}{\partial y} \biggr] </math>
	<math>~=</math>	<math>~(\zeta_z + 2\Omega_f)+~\cancel{C_1 \rho \Psi} - \frac{\partial }{\partial x}\biggl[\frac{\Psi_0}{\rho} \frac{2x}{a^2} \biggr] - \frac{\partial }{\partial y} \biggl[\frac{\Psi_0}{\rho} \frac{2y}{b^2} \biggr] </math>
	<math>~=</math>	<math>~\zeta_z + 2\Omega_f - \biggl[\frac{\Psi_0}{\rho} \frac{2}{a^2} \biggr] - \biggl[\frac{\Psi_0}{\rho} \frac{2}{b^2} \biggr] </math>
	<math>~=</math>	<math>~\zeta_z + 2\Omega_f - \frac{2\Psi_0}{\rho}\biggl[\frac{1}{a^2} + \frac{1}{b^2} \biggr] </math>
	<math>~=</math>	<math>~2\Omega_f \, .</math> Yes!

Trial #2

Still restricting our discussion to nonaxisymmetric, thin disks, let's try, <math>~\Psi = \Psi_0 (\rho/\rho_c)^2</math>, and

<math>~\rho</math>	<math>~=</math>	<math>~\rho_c \biggl\{1 - \biggl[ \frac{y^2}{b^2} + \frac{x^2}{a^2}\biggr]\biggr\} </math>
<math>~\Rightarrow~~~\frac{\partial^2 \rho}{\partial x^2} = -\frac{\partial}{\partial x}\biggl\{ \frac{2\rho_c x}{a^2}\biggr\} = - \frac{2\rho_c}{a^2}</math>	and,	<math>~\frac{\partial^2 \rho}{\partial y^2} = - \frac{\partial}{\partial y}\biggl\{ \frac{2\rho_c y}{b^2} \biggr\} = - \frac{2\rho_c}{b^2} \, .</math>

IMPORTANT NOTE (by Tohline on 22 September 2020): As I have come to appreciate today — after studying the relevant sections of both EFE and BT87 — this is an example of a heterogeneous density distribution whose gravitational potential has an analytic prescription. As is discussed in a separate chapter, the potential that it generates is sometimes referred to as a Ferrers potential, for the exponent, n = 1.

In our accompanying discussion we find that,

<math>~\frac{ \Phi_\mathrm{grav}(\bold{x})}{(-\pi G\rho_c)} </math>

<math>~ \frac{1}{2} I_\mathrm{BT} a_1^2 - \biggl(A_1 x^2 + A_2 y^2 +A_3 z^2 \biggr) ~+ \biggl( A_{12} x^2y^2 + A_{13} x^2z^2 + A_{23} y^2z^2\biggr) ~+ \frac{1}{6} \biggl(3A_{11}x^4 + 3A_{22}y^4 + 3A_{33}z^4 \biggr) \, , </math>

where,

for <math>~i \ne j</math>
<math>~A_{ij}</math>	<math>~\equiv</math>	<math>~-\frac{A_i-A_j}{(a_i^2 - a_j^2)} </math>
[ EFE, §21, Eq. (107) ]

for <math>~i = j</math>
<math>~2A_{ii} + \sum_{\ell = 1}^3 A_{i\ell}</math>	<math>~=</math>	<math>~\frac{2}{a_i} </math>
[ EFE, §21, Eq. (109) ]

More specifically, in the three cases where the indices, <math>~i=j</math>,

<math>~3A_{11}</math>	<math>~=</math>	<math>~ \frac{2}{a_1^2} - (A_{12} + A_{13}) \, , </math>
<math>~3A_{22}</math>	<math>~=</math>	<math>~ \frac{2}{a_2^2} - (A_{21} + A_{23}) \, , </math>
<math>~3A_{33}</math>	<math>~=</math>	<math>~ \frac{2}{a_3^2} - (A_{31} + A_{32}) \, . </math>

This means that,

<math>~\frac{ \partial \Psi}{\partial x} = \biggl( \frac{\Psi_0}{\rho_c^2}\biggr) 2\rho \frac{ \partial \rho}{\partial x} = -2\rho \biggl( \frac{2\rho_c x}{a^2} \biggr)\biggl( \frac{\Psi_0}{\rho_c^2}\biggr)</math>

and,

<math>~\frac{ \partial \Psi}{\partial y} = \biggl( \frac{\Psi_0}{\rho_c^2}\biggr)2\rho \frac{ \partial \rho}{\partial y} = -2\rho \biggl( \frac{2\rho_c y}{b^2} \biggr)\biggl( \frac{\Psi_0}{\rho_c^2}\biggr) \, .</math>

The Elliptic PDE Constraint Equation

<math>~ 2\Omega_f </math>	<math>~=</math>	<math>~C_0 \rho+~\cancel{C_1 \rho \Psi} + \frac{\partial }{\partial x}\biggl[\frac{1}{\rho} \frac{\partial \Psi}{\partial x} \biggr] + \frac{\partial }{\partial y} \biggl[\frac{1}{\rho} \frac{\partial \Psi}{\partial y} \biggr] </math>
	<math>~=</math>	<math>~C_0 \rho - \frac{\partial }{\partial x}\biggl[ \biggl( \frac{4 x}{a^2} \biggr)\biggl( \frac{\Psi_0}{\rho_c}\biggr) \biggr] - \frac{\partial }{\partial y} \biggl[ \biggl( \frac{4 y}{b^2} \biggr)\biggl( \frac{\Psi_0}{\rho_c}\biggr) \biggr] </math>
	<math>~=</math>	<math>~C_0 \rho - \frac{4\Psi_0}{\rho_c}\biggl[ \frac{1}{a^2} + \frac{1}{b^2} \biggr] </math>

The momentum density vector is governed by the relation,

<math>~\rho\bold{u}</math>	<math>~=</math>	<math>~ \boldsymbol{\hat\imath} \biggl[ \frac{\partial \Psi}{\partial y} \biggr] - \boldsymbol{\hat\jmath} \biggl[ \frac{\partial \Psi}{\partial x}\biggr] </math>
	<math>~=</math>	<math>~ - \boldsymbol{\hat\imath} \biggl[ 2\rho \biggl( \frac{2\rho_c y}{b^2} \biggr)\biggl( \frac{\Psi_0}{\rho_c^2}\biggr) \biggr] + \boldsymbol{\hat\jmath} \biggl[ 2\rho \biggl( \frac{2\rho_c x}{a^2} \biggr)\biggl( \frac{\Psi_0}{\rho_c^2}\biggr)\biggr] </math>
	<math>~=</math>	<math>~ \rho \biggl\{ - \boldsymbol{\hat\imath} \biggl[ \frac{4\Psi_0 y}{\rho_c b^2} \biggr] + \boldsymbol{\hat\jmath} \biggl[ \frac{4\Psi_0 x}{\rho_c a^2} \biggr] \biggr\} \, . </math>
[<math>~\Psi</math> has units of "density × length² per time"]

As above, let's see if the steady-state continuity equation is satisfied:

<math>~\nabla \cdot (\rho \bold{u})</math>	<math>~=</math>	<math>~ \frac{\partial (\rho u_x)}{\partial x} + \frac{\partial (\rho u_y)}{\partial y} </math>
	<math>~=</math>	<math>~ -\frac{\partial }{\partial x}\biggl[ \frac{4\Psi_0 y \rho}{\rho_c b^2} \biggr] + \frac{\partial }{\partial y}\biggl[ \frac{4\Psi_0 x \rho }{\rho_c a^2} \biggr] </math>
	<math>~=</math>	<math>~\frac{4\Psi_0}{\rho_c} \biggl\{ - \frac{y}{b^2} \cdot \frac{\partial \rho}{\partial x} + \frac{ x }{ a^2} \cdot \frac{\partial \rho}{\partial y} \biggr\} </math>
	<math>~=</math>	<math>~\frac{4\Psi_0}{\rho_c} \biggl\{ \frac{y}{b^2} \biggl[ \frac{2\rho_c x}{a^2} \biggr] - \frac{ x }{ a^2} \biggl[ \frac{2 \rho_c y}{b^2} \biggr] \biggr\} </math>
	<math>~=</math>	<math>~0 \, .</math> Yes!

Next, as above, let's determine the z-component of the vorticity and the vortensity:

<math>~\zeta_z</math>	<math>~=</math>	<math>~ \frac{\partial u_y}{\partial x} - \frac{\partial u_x}{\partial y} </math>
	<math>~=</math>	<math>~\frac{ 4\Psi_0}{\rho_c} \biggl\{ \frac{\partial }{\partial x}\biggl[ \frac{x}{a^2} \biggr] + \frac{\partial }{\partial y}\biggl[ \frac{y}{b^2 } \biggr] \biggr\} </math>
	<math>~=</math>	<math>~\frac{ 4\Psi_0}{\rho_c} \biggl[ \frac{1}{a^2} + \frac{1}{b^2} \biggr] \, . </math>
<math>~\Rightarrow~~~ \frac{(\zeta_z + 2\Omega_f)}{\rho}</math>	<math>~=</math>	<math>~\frac{1}{\rho} \biggl\{ \frac{ 4\Psi_0}{\rho_c} \biggl[ \frac{1}{a^2} + \frac{1}{b^2} \biggr] + 2\Omega_f \biggr\}\, . </math>

This means that,

<math>~ (\vec\zeta + 2\vec\Omega) \times \bold{u} </math>	<math>~=</math>	<math>~\biggl\{ \frac{(\zeta_z + 2\Omega_f)}{\rho} \biggr\} \boldsymbol{\hat{k}} \times (\rho \bold{u} )</math>
	<math>~=</math>	<math>~\frac{ 4\Psi_0 (\zeta_z + 2\Omega_f)}{\rho_c} \boldsymbol{\hat{k}} \times \biggl\{ - \boldsymbol{\hat\imath} \biggl[ \frac{y}{b^2} \biggr] + \boldsymbol{\hat\jmath} \biggl[ \frac{x}{a^2} \biggr] \biggr\} </math>
	<math>~=</math>	<math>~- \frac{ 4\Psi_0 (\zeta_z + 2\Omega_f)}{\rho_c} \biggl\{ \boldsymbol{\hat\imath} \biggl[ \frac{x}{a^2} \biggr] + \boldsymbol{\hat\jmath} \biggl[ \frac{y}{b^2} \biggr] \biggr\} </math>

Now, let's examine the gradient of <math>~F_B(\Psi)</math>.

<math>~\nabla F_B(\Psi)</math>	<math>~=</math>	<math>~\biggl[ \boldsymbol{\hat\imath} \frac{\partial}{\partial x} + \boldsymbol{\hat\jmath} \frac{\partial}{\partial y} \biggr]C_0 \Psi</math>
	<math>~=</math>	<math>~ \frac{ C_0 \Psi_0}{\rho_c^2} \biggl[ \boldsymbol{\hat\imath} \frac{\partial \rho^2}{\partial x} + \boldsymbol{\hat\jmath} \frac{\partial \rho^2}{\partial y} \biggr] </math>
	<math>~=</math>	<math>~ \frac{ 2C_0 \rho \Psi_0}{\rho_c^2} \biggl\{ -\boldsymbol{\hat\imath} \biggl[ \frac{2 \rho_c x}{a^2} \biggr] - \boldsymbol{\hat\jmath} \biggl[ \frac{2 \rho_c y}{b^2} \biggr] \biggr\} </math>
	<math>~=</math>	<math>~ - \frac{ 4 C_0 \rho \Psi_0}{\rho_c} \biggl\{ \boldsymbol{\hat\imath} \biggl[ \frac{x}{a^2} \biggr] + \boldsymbol{\hat\jmath} \biggl[ \frac{y}{b^2} \biggr] \biggr\} \, , </math>

which is identical to the immediately preceding expression if we set,

<math>~ \frac{(\zeta_z + 2\Omega_f)}{\rho} \, . </math>

Continuing with the above examination of the elliptic PDE "constraint" equation, we find that,

<math>~2\Omega_f</math>	<math>~=</math>	<math>~C_0 \rho - \zeta_z </math>
	<math>~=</math>	<math>~(2\Omega_f + \zeta_z) - \zeta_z </math>
	<math>~=</math>	<math>~2\Omega_f \, .</math> Hooray!

Trial #3

If the density distribution has been specified, then <math>~\Psi</math> is the "stream-function" from which all rotating-frame velocities are determined. Specifically,

<math>~ \frac{1}{\rho} \biggl[ \boldsymbol{\hat\imath} \frac{\partial \Psi}{\partial y} - \boldsymbol{\hat\jmath} \frac{\partial \Psi}{\partial x} \biggr] \, . </math>

Most importantly, as has been detailed above, the term on the left-hand-side of the steady-state Euler equation becomes,

<math>~(\vec\zeta + 2\vec\Omega) \times \bold{u}</math>	<math>~=</math>	<math>~\biggl[ \frac{ (\zeta_z + 2\Omega_f)}{\rho} \biggr] \nabla\Psi \, ,</math>
Saied W. Andalib (1998), §4.2, p. 83, Eq. (4.13)

where,

<math>~ -~\frac{\partial }{\partial x}\biggl[\frac{1}{\rho} \frac{\partial \Psi}{\partial x} \biggr] - \frac{\partial }{\partial y} \biggl[\frac{1}{\rho} \frac{\partial \Psi}{\partial y} \biggr] \, . </math>

Still restricting our discussion to infinitesimally thin, nonaxisymmetric disks, let's assume that,

<math>~\rho_c \biggl[1 - \biggl( \frac{x^2}{a^2} + \frac{y^2}{b^2} \biggr) \biggr] </math>

<math>~\Rightarrow~~~\frac{\partial \rho}{\partial x} = -\biggl( \frac{2\rho_c x}{a^2} \biggr)</math>	and,	<math>~\frac{\partial \rho}{\partial y} = -\biggl( \frac{2\rho_c y}{b^2} \biggr) \, ;</math>
<math>~\Rightarrow~~~\frac{\partial^2 \rho}{\partial x^2} = -\frac{\partial}{\partial x}\biggl\{ \frac{2\rho_c x}{a^2}\biggr\} = - \frac{2\rho_c}{a^2}</math>	and,	<math>~\frac{\partial^2 \rho}{\partial y^2} = - \frac{\partial}{\partial y}\biggl\{ \frac{2\rho_c y}{b^2} \biggr\} = - \frac{2\rho_c}{b^2} \, .</math>

And, let's assume that,

<math>~\Psi_0 \biggl( \frac{\rho}{\rho_c} \biggr)^q \, ,</math>

<math>~\Rightarrow~~~ \frac{\partial\Psi}{\partial x} </math>	<math>~=</math>	<math>~ \rho^{q-1}\biggl(\frac{q \Psi_0}{\rho_c^q}\biggr) \frac{\partial\rho}{\partial x} = - \rho^{q-1}\biggl(\frac{q \Psi_0}{\rho_c^q}\biggr) \biggl(\frac{2\rho_c x}{a^2} \biggr) = - \biggl( \frac{\rho}{\rho_c}\biggr)^{q-1} \biggl(\frac{2q \Psi_0x}{a^2} \biggr) \, ; </math>
and, similarly, <math>~\frac{\partial\Psi}{\partial y} </math>	<math>~=</math>	<math>~ - \biggl( \frac{\rho}{\rho_c}\biggr)^{q-1} \biggl(\frac{2q \Psi_0 y}{b^2} \biggr) \, . </math>

This means that,

<math>~\bold{u}</math>	<math>~=</math>	<math>~-~ \frac{1}{\rho} \biggl\{ \boldsymbol{\hat\imath} \biggl[ \biggl( \frac{\rho}{\rho_c}\biggr)^{q-1} \biggl(\frac{2q \Psi_0 y}{b^2} \biggr) \biggr] - \boldsymbol{\hat\jmath} \biggl[ \biggl( \frac{\rho}{\rho_c}\biggr)^{q-1} \biggl(\frac{2q \Psi_0x}{a^2} \biggr) \biggr] \biggr\} </math>
<math>~\Rightarrow~~~\bold{u}\cdot \bold{u} </math>	<math>~=</math>	<math>~ \frac{1}{\rho^2} \biggl\{ \biggl[ \biggl( \frac{\rho}{\rho_c}\biggr)^{q-1} \biggl(\frac{2q \Psi_0 y}{b^2} \biggr) \biggr]^2 + \biggl[ \biggl( \frac{\rho}{\rho_c}\biggr)^{q-1} \biggl(\frac{2q \Psi_0x}{a^2} \biggr) \biggr]^2 \biggr\} \, . </math>

It also means that,

<math>~\zeta_z </math>	<math>~=</math>	<math>~ \frac{\partial }{\partial x}\biggl[\biggl( \frac{\rho}{\rho_c}\biggr)^{q-2} \biggl(\frac{2q \Psi_0 x}{\rho_c a^2} \biggr) \biggr] + \frac{\partial }{\partial y} \biggl[\biggl( \frac{\rho}{\rho_c}\biggr)^{q-2} \biggl(\frac{2q \Psi_0 y}{\rho_c b^2} \biggr) \biggr] </math>
	<math>~=</math>	<math>~ \biggl( \frac{\rho}{\rho_c}\biggr)^{q-2} \biggl(\frac{2q \Psi_0 }{\rho_c a^2} \biggr) ~+ ~ \biggl(\frac{2q \Psi_0 x}{\rho_c a^2} \biggr) \frac{\partial }{\partial x} \biggl( \frac{\rho}{\rho_c}\biggr)^{q-2} ~+ ~ \biggl( \frac{\rho}{\rho_c}\biggr)^{q-2} \biggl(\frac{2q \Psi_0 }{\rho_c b^2} \biggr) ~+ ~ \biggl(\frac{2q \Psi_0 y}{\rho_c b^2} \biggr) \frac{\partial }{\partial y} \biggl( \frac{\rho}{\rho_c}\biggr)^{q-2} </math>
	<math>~=</math>	<math>~ \biggl( \frac{\rho}{\rho_c}\biggr)^{q-2} \biggl(\frac{2q \Psi_0 }{\rho_c } \biggr) \biggl[ \frac{1}{a^2} + \frac{1}{b^2}\biggr] ~+ ~ \biggl(\frac{2q \Psi_0 x}{\rho_c a^2} \biggr) \biggl[ (q-2) \rho_c^{2-q} \rho^{q-3} \frac{\partial \rho}{\partial x}\biggr] ~+ ~ \biggl(\frac{2q \Psi_0 y}{\rho_c b^2} \biggr) \biggl[ (q-2) \rho_c^{2-q} \rho^{q-3} \frac{\partial \rho}{\partial y}\biggr] </math>
	<math>~=</math>	<math>~ \biggl( \frac{\rho}{\rho_c}\biggr)^{q-2} \biggl(\frac{2q \Psi_0 }{\rho_c } \biggr) \biggl[ \frac{1}{a^2} + \frac{1}{b^2}\biggr] ~- ~ \biggl(\frac{2q \Psi_0 x}{\rho_c^2 a^2} \biggr) \biggl[ (q-2) \biggl(\frac{\rho}{\rho_c}\biggr)^{q-3} \frac{2\rho_c x}{a^2}\biggr] ~- ~ \biggl(\frac{2q \Psi_0 y}{\rho_c^2 b^2} \biggr) \biggl[ (q-2) \biggl(\frac{\rho}{\rho_c}\biggr)^{q-3} \frac{2\rho_c y}{b^2}\biggr] </math>
	<math>~=</math>	<math>~\frac{q\Psi_0}{\rho_c} \biggl\{ 2\biggl( \frac{\rho}{\rho_c}\biggr)^{q-2} \biggl[ \frac{1}{a^2} + \frac{1}{b^2}\biggr] ~- ~ \biggl(\frac{2 x}{ a^2} \biggr) \biggl[ (q-2) \biggl(\frac{\rho}{\rho_c}\biggr)^{q-3} \frac{2x}{a^2}\biggr] ~- ~ \biggl(\frac{2 y}{ b^2} \biggr) \biggl[ (q-2) \biggl(\frac{\rho}{\rho_c}\biggr)^{q-3} \frac{2 y}{b^2}\biggr] \biggr\} </math>
	<math>~=</math>	<math>~ \frac{q\Psi_0}{\rho_c} \biggl(\frac{\rho}{\rho_c}\biggr)^{q-3} \biggl[ 2\biggl( \frac{\rho}{\rho_c}\biggr) \biggl( \frac{1}{a^2} + \frac{1}{b^2}\biggr) ~- ~ 4(q-2) \biggl(\frac{x^2}{ a^4} \biggr) ~- ~ 4(q-2) \biggl(\frac{y^2}{ b^4} \biggr) \biggr] </math>

Hence,

<math>~(\vec\zeta + 2\vec\Omega) \times \bold{u}</math>	<math>~=</math>	<math>~ \biggl[ \boldsymbol{\hat\imath} \frac{\partial \Psi}{\partial x} + \boldsymbol{\hat\jmath} \frac{\partial \Psi}{\partial y} \biggr] \biggl\{ \frac{ 2\Omega_f}{\rho_c}\biggl( \frac{\rho}{\rho_c}\biggr)^{-1} + \biggl( \frac{\rho}{\rho_c}\biggr)^{-1} \frac{ \zeta_z }{\rho_c} \biggr\} </math>
	<math>~=</math>	<math>~-~ \biggl[ \boldsymbol{\hat\imath} \biggl( \frac{\rho}{\rho_c}\biggr)^{q-1} \biggl(\frac{2q \Psi_0 x}{a^2} \biggr) + \boldsymbol{\hat\jmath} \biggl( \frac{\rho}{\rho_c}\biggr)^{q-1} \biggl(\frac{2q \Psi_0 y}{b^2} \biggr) \biggr] \biggl\{ \frac{ 2\Omega_f}{\rho_c}\biggl( \frac{\rho}{\rho_c}\biggr)^{-1} </math>
		<math>~ + \frac{q\Psi_0}{\rho_c^2} \biggl(\frac{\rho}{\rho_c}\biggr)^{q-4} \biggl[ 2\biggl( \frac{\rho}{\rho_c}\biggr) \biggl( \frac{1}{a^2} + \frac{1}{b^2}\biggr) ~- ~ 4(q-2) \biggl(\frac{x^2}{ a^4} \biggr) ~- ~ 4(q-2) \biggl(\frac{y^2}{ b^4} \biggr) \biggr] \biggr\} </math>
	<math>~=</math>	<math>~ -~2q\Psi_0 \biggl[ \boldsymbol{\hat\imath} \biggl(\frac{x}{a^2} \biggr) + \boldsymbol{\hat\jmath} \biggl(\frac{y}{b^2} \biggr) \biggr] \biggl( \frac{\rho}{\rho_c}\biggr)^{q-2} \biggl\{ \frac{ 2\Omega_f}{\rho_c} </math>
		<math>~ + \frac{q\Psi_0}{\rho_c^2} \biggl(\frac{\rho}{\rho_c}\biggr)^{q-3} \biggl[ 2\biggl( \frac{\rho}{\rho_c}\biggr) \biggl( \frac{1}{a^2} + \frac{1}{b^2}\biggr) ~- ~ 4(q-2) \biggl(\frac{x^2}{ a^4} \biggr) ~- ~ 4(q-2) \biggl(\frac{y^2}{ b^4} \biggr) \biggr] \biggr\} \, . </math>

Exponent q = 2

Notice that, if <math>~q = 2</math>,

<math>~\biggl[ (\vec\zeta + 2\vec\Omega) \times \bold{u} \biggr]_{q=2}</math>

<math>~ -~\frac{8\Psi_0}{\rho_c}\biggl\{ \Omega_f + \frac{2\Psi_0}{\rho_c} \biggl( \frac{1}{a^2} + \frac{1}{b^2}\biggr) \biggr\} \biggl[ \boldsymbol{\hat\imath} \biggl(\frac{x}{a^2} \biggr) + \boldsymbol{\hat\jmath} \biggl(\frac{y}{b^2} \biggr) \biggr] \, . </math>

Now, if we choose a function,

<math>~ \frac{D_1}{2} \Psi^{1 / 2} = \frac{D_1}{2} \Psi_0^{1 / 2} \biggl(\frac{\rho}{\rho_c}\biggr) = \frac{D_1}{2} \Psi_0^{1 / 2} \biggl[1 - \biggl( \frac{x^2}{a^2} + \frac{y^2}{b^2} \biggr) \biggr] \, ,</math>

we obtain,

<math>~\nabla F_B(\Psi)</math>	<math>~=</math>	<math>~ -~\frac{D_1}{2} \Psi_0^{1 / 2} \cdot \nabla \biggl( \frac{x^2}{a^2} + \frac{y^2}{b^2} \biggr) </math>
	<math>~=</math>	<math>~ -~D_1 \Psi_0^{1 / 2} \biggl[ \boldsymbol{\hat\imath} \biggl(\frac{x}{a^2} \biggr) + \boldsymbol{\hat\jmath} \biggl(\frac{y}{b^2} \biggr) \biggr] \, . </math>

This is consistent with the elliptic PDE constraint if,

<math>~\frac{8\Psi_0^{1 / 2}}{\rho_c}\biggl\{ \Omega_f + \frac{2\Psi_0}{\rho_c} \biggl( \frac{1}{a^2} + \frac{1}{b^2}\biggr) \biggr\} \, .</math>

Also if <math>~q = 2</math>, we have,

<math>~\bold{u}</math>	<math>~=</math>	<math>~ \frac{1}{\rho} \biggl\{ \boldsymbol{\hat\imath} \frac{\partial \Psi}{\partial y} - \boldsymbol{\hat\jmath} \frac{\partial \Psi}{\partial x} \biggr\} </math>
	<math>~=</math>	<math>~ \frac{1}{\rho} \biggl\{ -\boldsymbol{\hat\imath} \biggl( \frac{\rho}{\rho_c}\biggr) \biggl(\frac{2q \Psi_0 y}{b^2} \biggr) + \boldsymbol{\hat\jmath} \biggl( \frac{\rho}{\rho_c}\biggr) \biggl(\frac{2q \Psi_0 x}{a^2} \biggr) \biggr\} </math>
	<math>~=</math>	<math>~ \biggl\{ -\boldsymbol{\hat\imath} \biggl(\frac{4 \Psi_0 y}{\rho_c b^2} \biggr) + \boldsymbol{\hat\jmath} \biggl(\frac{4 \Psi_0 x}{\rho_c a^2} \biggr) \biggr\} </math>
<math>~\Rightarrow ~~~ \bold{u}\cdot \bold{u}</math>	<math>~=</math>	<math>~ \biggl(\frac{4 \Psi_0 y}{\rho_c b^2} \biggr)^2 + \biggl(\frac{4 \Psi_0 x}{\rho_c a^2} \biggr)^2 </math>
	<math>~=</math>	<math>~ \biggl( \frac{4 \Psi_0}{\rho_c}\biggr)^2 \biggl[ \frac{x^2}{a^4} + \frac{y^2}{b^4} \biggr] \, . </math>

Keep in mind that, as discussed above, we are trying to satisfy the scalar Bernoulli relation,

<math>~C_B - H - \Phi_\mathrm{grav}</math>	<math>~=</math>	<math>~ F_B(\Psi) + \frac{1}{2}u^2 - \frac{1}{2}\Omega_f^2 (x^2 + y^2) </math>
	<math>~=</math>	<math>~ \frac{4\Psi_0}{\rho_c}\biggl\{ \Omega_f + \frac{2\Psi_0}{\rho_c} \biggl( \frac{1}{a^2} + \frac{1}{b^2}\biggr) \biggr\}\biggl[1 - \biggl( \frac{x^2}{a^2} + \frac{y^2}{b^2} \biggr) \biggr] + \frac{1}{2}\biggl( \frac{4 \Psi_0}{\rho_c}\biggr)^2 \biggl[ \frac{x^2}{a^4} + \frac{y^2}{b^4} \biggr] - \frac{1}{2}\Omega_f^2 (x^2 + y^2) \, . </math>

The right-hand-side of this expression does not appear to be rich enough to balance the gravitational potential (on the left-hand-side) which, as detailed above, contains <math>~x^2 y^2</math> and <math>~x^4</math> and <math>~y^4</math> terms.

Exponent q = 3

Alternatively, if <math>~q = 3</math>,

<math>~\biggl[ (\vec\zeta + 2\vec\Omega) \times \bold{u} \biggr]_{q=3}</math>	<math>~=</math>	<math>~-~6\Psi_0 \biggl[ \boldsymbol{\hat\imath} \biggl(\frac{x}{a^2} \biggr) + \boldsymbol{\hat\jmath} \biggl(\frac{y}{b^2} \biggr) \biggr] \biggl( \frac{\rho}{\rho_c}\biggr) \biggl\{ \frac{ 2\Omega_f}{\rho_c} + \frac{3\Psi_0}{\rho_c^2} \biggl[ 2\biggl( \frac{\rho}{\rho_c}\biggr) \biggl( \frac{1}{a^2} + \frac{1}{b^2}\biggr) ~- ~ 4 \biggl(\frac{x^2}{ a^4} \biggr) ~- ~ 4 \biggl(\frac{y^2}{ b^4} \biggr) \biggr] \biggr\} </math>
	<math>~=</math>	<math>~-~ 6\Psi_0 \biggl\{ \frac{ 2\Omega_f}{\rho_c} \biggl( \frac{\rho}{\rho_c}\biggr) - \frac{12\Psi_0}{\rho_c^2} \biggl[ \frac{x^2}{ a^4} + \frac{y^2}{ b^4} \biggr]\biggl( \frac{\rho}{\rho_c}\biggr) + \frac{6\Psi_0}{\rho_c^2} \biggl( \frac{1}{a^2} + \frac{1}{b^2}\biggr) \biggl( \frac{\rho}{\rho_c}\biggr)^2 \biggr\} \boldsymbol{\hat{f}} \, , </math>

where,

<math>~\boldsymbol{\hat{f}}</math>

<math>~\equiv</math>

<math>~ \biggl[ \boldsymbol{\hat\imath} \biggl(\frac{x}{a^2} \biggr) + \boldsymbol{\hat\jmath} \biggl(\frac{y}{b^2} \biggr) \biggr] \, . </math>

If we choose a function,

<math>~F_B(\Psi)</math>	<math>~=</math>	<math>~ D_1\Psi^{2/3} + D_2 \Psi </math>
	<math>~=</math>	<math>~ D_1\Psi_0^{2/3} \biggl(\frac{\rho}{\rho_c}\biggr)^2 + D_2 \Psi_0 \biggl(\frac{\rho}{\rho_c}\biggr)^3 \, , </math>

we obtain,

<math>~\nabla F_B(\Psi)</math>	<math>~=</math>	<math>~ -~\biggl[ 2D_1 \Psi_0^{2/3} \biggl(\frac{\rho}{\rho_c}\biggr) + 3D_2\Psi_0 \biggl( \frac{\rho}{\rho_c}\biggr)^2\biggr] \cdot \nabla \biggl( \frac{x^2}{a^2} + \frac{y^2}{b^2} \biggr) </math>
	<math>~=</math>	<math>~ -~\biggl\{ 2D_1 \Psi_0^{2/3} \biggl[1 - \biggl(\frac{x^2}{a^2} + \frac{y^2}{b^2}\biggr) \biggr] + 3D_2\Psi_0 \biggl[1 - \biggl(\frac{x^2}{a^2} + \frac{y^2}{b^2}\biggr) \biggr]^2 \biggr\} \cdot \biggl[ \boldsymbol{\hat\imath} \biggl(\frac{2x}{a^2} \biggr) + \boldsymbol{\hat\jmath} \biggl(\frac{2y}{b^2} \biggr) \biggr] </math>
	<math>~=</math>	<math>~ -~\biggl\{ 2D_1 \Psi_0^{2/3} \biggl[1 - \biggl(\frac{x^2}{a^2} + \frac{y^2}{b^2}\biggr) \biggr] + 3D_2\Psi_0 \biggl[1 - 2\biggl(\frac{x^2}{a^2} + \frac{y^2}{b^2}\biggr) + \biggl(\frac{x^2}{a^2} + \frac{y^2}{b^2}\biggr)^2 \biggr] \biggr\} \cdot \biggl[ \boldsymbol{\hat\imath} \biggl(\frac{2x}{a^2} \biggr) + \boldsymbol{\hat\jmath} \biggl(\frac{2y}{b^2} \biggr) \biggr] </math>
	<math>~=</math>	<math>~ -~\biggl\{ \biggl(2D_1 \Psi_0^{2/3} + 3D_2\Psi_0 \biggr) - (2D_1 \Psi_0^{2/3} +6D_2\Psi_0)\biggl(\frac{x^2}{a^2} + \frac{y^2}{b^2}\biggr) + 3D_2\Psi_0 \biggl(\frac{x^2}{a^2} + \frac{y^2}{b^2}\biggr)^2 \biggr\} \boldsymbol{\hat{f}} </math>
	<math>~=</math>	<math>~ -~\biggl\{ \biggl(2D_1 \Psi_0^{2/3} + 3D_2\Psi_0 \biggr) - (2D_1 \Psi_0^{2/3} +6D_2\Psi_0)\biggl(\frac{x^2}{a^2} + \frac{y^2}{b^2}\biggr) + 3D_2\Psi_0 \biggl[ \frac{x^4}{a^4} + 2\biggl( \frac{x^2 y^2}{a^2 b^2} \biggr) + \frac{y^4}{b^4}\biggr] \biggr\} \boldsymbol{\hat{f}} \, . </math>

Let's reorganize and expand the terms in both of these expressions in order to ascertain whether or not they can be matched. First …

<math>~\biggl[ (\vec\zeta + 2\vec\Omega) \times \bold{u} \biggr]_{q=3}</math>	<math>~=</math>	<math>~-~ 6\Psi_0 \biggl\{ \frac{ 2\Omega_f}{\rho_c} \biggl[ 1 - \biggl( \frac{x^2}{a^2} + \frac{y^2}{b^2} \biggr) \biggr] - \frac{12\Psi_0}{\rho_c^2} \biggl[ \frac{x^2}{ a^4} + \frac{y^2}{ b^4} \biggr] \biggl[ 1 - \biggl( \frac{x^2}{a^2} + \frac{y^2}{b^2} \biggr) \biggr] </math>
		<math>~ + \frac{6\Psi_0}{\rho_c^2} \biggl( \frac{1}{a^2} + \frac{1}{b^2}\biggr) \biggl[ 1 - \biggl( \frac{x^2}{a^2} + \frac{y^2}{b^2} \biggr) \biggr]^2 \biggr\} \boldsymbol{\hat{f}} \, , </math>
	<math>~=</math>	<math>~-~ 6\Psi_0 \biggl\{ \frac{ 2\Omega_f}{\rho_c} - \frac{ 2\Omega_f}{\rho_c} \biggl( \frac{x^2}{a^2} + \frac{y^2}{b^2} \biggr) - \frac{12\Psi_0}{\rho_c^2} \biggl[ \frac{x^2}{ a^4} + \frac{y^2}{ b^4} \biggr] + \frac{12\Psi_0}{\rho_c^2} \biggl[ \frac{x^2}{ a^4} + \frac{y^2}{ b^4} \biggr] \biggl( \frac{x^2}{a^2} + \frac{y^2}{b^2} \biggr) </math>
		<math>~ + \frac{6\Psi_0}{\rho_c^2} \biggl( \frac{1}{a^2} + \frac{1}{b^2}\biggr) \biggl[ 1 - 2\biggl( \frac{x^2}{a^2} + \frac{y^2}{b^2} \biggr) + \biggl( \frac{x^2}{a^2} + \frac{y^2}{b^2} \biggr)^2 \biggr] \biggr\} \boldsymbol{\hat{f}} \, , </math>
	<math>~=</math>	<math>~ -~ 6\Psi_0 \biggl\{ \biggl[ \frac{ 2\Omega_f}{\rho_c} +\frac{6\Psi_0}{\rho_c^2} \cdot \frac{(a^2 + b^2)}{a^2b^2} \biggr] - \biggl[ \frac{ 2\Omega_f}{\rho_c} + \frac{12\Psi_0}{\rho_c^2}\biggl( \frac{1}{a^2} + \frac{1}{b^2}\biggr)\biggr]\biggl( \frac{x^2}{a^2} + \frac{y^2}{b^2} \biggr) </math>
		<math>~ - \frac{12\Psi_0}{\rho_c^2} \biggl( \frac{x^2}{ a^4} + \frac{y^2}{ b^4} \biggr) + \frac{12\Psi_0}{\rho_c^2} \biggl( \frac{x^4}{ a^6} + \frac{x^2 y^2}{ a^4 b^2} + \frac{x^2 y^2}{ a^2 b^4} + \frac{y^4}{ b^6}\biggr) + \frac{6\Psi_0}{\rho_c^2}\biggl( \frac{1}{a^2} + \frac{1}{b^2}\biggr)\biggl( \frac{x^4}{a^4} + \frac{2x^2 y^2}{a^2 b^2}+ \frac{y^4}{b^4} \biggr) \biggr\} \boldsymbol{\hat{f}} \, , </math>

In order for the zeroth-order terms to match, we need,

<math>~2D_1 \Psi_0^{2/3} + 3D_2\Psi_0 </math>	<math>~=</math>	<math>~ 6\Psi_0 \biggl\{ \frac{ 2\Omega_f}{\rho_c} +\frac{6\Psi_0}{\rho_c^2} \cdot \frac{(a^2 + b^2)}{a^2b^2} \biggr\} </math>
<math>~\Rightarrow ~~~ 2D_1 \Psi_0^{2/3} </math>	<math>~=</math>	<math>~ 3\Psi_0 \biggl\{ \frac{ 4\Omega_f \rho_c a^2 b^2}{\rho_c^2 a^2 b^2} +\frac{12\Psi_0(a^2 + b^2)}{\rho_c^2 a^2 b^2} \biggr\} - 3D_2\Psi_0 </math>
	<math>~=</math>	<math>~ 3\Psi_0 \biggl\{ \frac{ 4\Omega_f \rho_c a^2 b^2}{\rho_c^2 a^2 b^2} +\frac{12\Psi_0(a^2 + b^2)}{\rho_c^2 a^2 b^2} -D_2\biggr\} \, . </math>

Then we also need,

<math>~ 6D_2\Psi_0</math>	<math>~=</math>	<math>~6\Psi_0 \biggl[ \frac{ 2\Omega_f}{\rho_c} + \frac{12\Psi_0}{\rho_c^2}\biggl( \frac{1}{a^2} + \frac{1}{b^2}\biggr)\biggr] - 2D_1 \Psi_0^{2/3} </math>
	<math>~=</math>	<math>~6\Psi_0 \biggl[ \frac{ 2\Omega_f}{\rho_c} + \frac{12\Psi_0}{\rho_c^2}\biggl( \frac{1}{a^2} + \frac{1}{b^2}\biggr)\biggr] - 3\Psi_0 \biggl\{ \frac{ 4\Omega_f \rho_c a^2 b^2}{\rho_c^2 a^2 b^2} +\frac{12\Psi_0(a^2 + b^2)}{\rho_c^2 a^2 b^2} -D_2\biggr\} </math>
<math>~\Rightarrow~~~ 3D_2\Psi_0</math>	<math>~=</math>	<math>~6\Psi_0 \biggl[ \frac{ 2\Omega_f \rho_c a^2b^2}{\rho_c^2 a^2 b^2} + \frac{12\Psi_0}{\rho_c^2}\frac{(a^2 + b^2)}{a^2b^2} \biggr] - 6\Psi_0 \biggl[ \frac{ 2\Omega_f \rho_c a^2 b^2}{\rho_c^2 a^2 b^2} +\frac{6\Psi_0(a^2 + b^2)}{\rho_c^2 a^2 b^2} \biggr] </math>
	<math>~=</math>	<math>~6\Psi_0 \biggl[ \frac{6\Psi_0(a^2 + b^2)}{\rho_c^2a^2b^2}\biggr] </math>
<math>~\Rightarrow~~~ D_2</math>	<math>~=</math>	<math>~\frac{12\Psi_0(a^2 + b^2)}{\rho_c^2a^2b^2} \, . </math>

Finally, we need,

<math>~3D_2\Psi_0 \biggl[ \frac{x^4}{a^4} + 2\biggl( \frac{x^2 y^2}{a^2 b^2} \biggr) + \frac{y^4}{b^4}\biggr] \cdot \frac{1}{6\Psi_0} </math>	<math>~=</math>	<math>~ - \frac{12\Psi_0}{\rho_c^2} \biggl( \frac{x^2}{ a^4} + \frac{y^2}{ b^4} \biggr) + \frac{12\Psi_0}{\rho_c^2} \biggl( \frac{x^4}{ a^6} + \frac{x^2 y^2}{ a^4 b^2} + \frac{x^2 y^2}{ a^2 b^4} + \frac{y^4}{ b^6}\biggr) + \frac{6\Psi_0}{\rho_c^2}\biggl( \frac{1}{a^2} + \frac{1}{b^2}\biggr)\biggl( \frac{x^4}{a^4} + \frac{2x^2 y^2}{a^2 b^2}+ \frac{y^4}{b^4} \biggr) </math>
<math>~\Rightarrow ~~~ \frac{(a^2 + b^2)}{a^2b^2} \biggl[ \frac{x^4}{a^4} + 2\biggl( \frac{x^2 y^2}{a^2 b^2} \biggr) + \frac{y^4}{b^4}\biggr] </math>	<math>~=</math>	<math>~ - 2 \biggl( \frac{x^2}{ a^4} + \frac{y^2}{ b^4} \biggr) + 2 \biggl( \frac{x^4}{ a^6} + \frac{x^2 y^2}{ a^4 b^2} + \frac{x^2 y^2}{ a^2 b^4} + \frac{y^4}{ b^6}\biggr) + \frac{(a^2 + b^2)}{a^2b^2} \biggl( \frac{x^4}{a^4} + \frac{2x^2 y^2}{a^2 b^2}+ \frac{y^4}{b^4} \biggr) </math>
<math>~\Rightarrow ~~~\biggl( \frac{x^2}{ a^4} + \frac{y^2}{ b^4} \biggr) </math>	<math>~=</math>	<math>~ \biggl( \frac{x^4}{ a^6} + \frac{x^2 y^2}{ a^4 b^2} + \frac{x^2 y^2}{ a^2 b^4} + \frac{y^4}{ b^6}\biggr) </math>
<math>~\Rightarrow ~~~\frac{x^2}{ a^4} + \frac{y^2}{ b^4} </math>	<math>~=</math>	<math>~ \frac{x^4}{ a^6} + \frac{y^4}{ b^6} + \frac{x^2 y^2(a^2 + b^2) }{a^4 b^4} </math>
<math>~\Rightarrow ~~~ x^4 b^6 + y^4 a^6 </math>	<math>~=</math>	<math>~ a^2 b^2 [x^2 b^4 - (a^2+b^2)x^2 y^2 + y^2a^4] </math>
	<math>~=</math>	<math>~ a^2 b^2 [x^4 b^4 - (a^2+b^2)x^2 y^2 + y^4a^4] +a^2b^2[x^2b^4 + y^2a^4] -a^2b^2[x^4b^4 + y^4a^4] </math>

Keeping in mind that,

<math>~\Psi_0 \biggl( \frac{\rho}{\rho_c} \biggr)^q \, ,</math>

and that, after setting <math>~q = 3</math>, we have chosen,

<math>~F_B(\Psi)</math>	<math>~=</math>	<math>~ D_1\Psi^{2/3} + D_2 \Psi </math>
	<math>~=</math>	<math>~ D_1\Psi_0^{2/3} \biggl(\frac{\rho}{\rho_c}\biggr)^2 + D_2 \Psi_0 \biggl(\frac{\rho}{\rho_c}\biggr)^3 \, , </math>

let's try again.

<math>~\biggl[ (\vec\zeta + 2\vec\Omega) \times \bold{u} \biggr]_{q=3} - \nabla F_B(\Psi)</math>	<math>~=</math>	<math>~-~ 6\Psi_0 \biggl\{ \frac{ 2\Omega_f}{\rho_c} \biggl( \frac{\rho}{\rho_c}\biggr) - \frac{12\Psi_0}{\rho_c^2} \biggl[ \frac{x^2}{ a^4} + \frac{y^2}{ b^4} \biggr]\biggl( \frac{\rho}{\rho_c}\biggr) + \frac{6\Psi_0}{\rho_c^2} \biggl( \frac{1}{a^2} + \frac{1}{b^2}\biggr) \biggl( \frac{\rho}{\rho_c}\biggr)^2 \biggr\} \cdot \boldsymbol{\hat{f}} </math>
		<math>~ +~\biggl[ 2D_1 \Psi_0^{2/3} \biggl(\frac{\rho}{\rho_c}\biggr) + 3D_2\Psi_0 \biggl( \frac{\rho}{\rho_c}\biggr)^2\biggr] \cdot \nabla \biggl( \frac{x^2}{a^2} + \frac{y^2}{b^2} \biggr) </math>
	<math>~=</math>	<math>~-~ 3\Psi_0 \biggl\{ \frac{ 2\Omega_f}{\rho_c} - \frac{12\Psi_0}{\rho_c^2} \biggl[ \frac{x^2}{ a^4} + \frac{y^2}{ b^4} \biggr] + \frac{6\Psi_0}{\rho_c^2} \biggl( \frac{1}{a^2} + \frac{1}{b^2}\biggr) \biggl( \frac{\rho}{\rho_c}\biggr) \biggr\} \cdot 2\biggl(\frac{\rho}{\rho_c}\biggr)\boldsymbol{\hat{f}} </math>
		<math>~ +~\biggl[ 2D_1 \Psi_0^{2/3} + 3D_2\Psi_0 \biggl( \frac{\rho}{\rho_c}\biggr)\biggr] \cdot 2\biggl(\frac{\rho}{\rho_c}\biggr)\boldsymbol{\hat{f}} \, . </math>

Now, set …

<math>~2D_1 \Psi_0^{2/3}</math>	<math>~=</math>	<math>~\frac{6\Psi_0 \Omega_f}{\rho_c}</math>
<math>~\Rightarrow ~~~ D_1 </math>	<math>~=</math>	<math>~\frac{3\Psi_0^{1 / 3} \Omega_f}{\rho_c} \, ;</math>

and, set …

<math>~3D_2 \Psi_0</math>	<math>~=</math>	<math>~ \frac{18 \Psi_0^2(a^2 + b^2)}{\rho_c^2 a^2 b^2} </math>
<math>~\Rightarrow~~~ D_2 </math>	<math>~=</math>	<math>~ \frac{6 \Psi_0(a^2 + b^2)}{\rho_c^2 a^2 b^2} \, . </math>

We then have,

<math>~\biggl[ (\vec\zeta + 2\vec\Omega) \times \bold{u} \biggr]_{q=3} - \nabla F_B(\Psi)</math>	<math>~=</math>	<math>~-~ 3\Psi_0 \biggl\{ - \frac{12\Psi_0}{\rho_c^2} \biggl[ \frac{x^2}{ a^4} + \frac{y^2}{ b^4} \biggr] \biggr\} \cdot 2\biggl(\frac{\rho}{\rho_c}\biggr)\boldsymbol{\hat{f}} </math>
	<math>~=</math>	<math>~ \frac{72\Psi_0^2}{\rho_c^2} \biggl[ \frac{x^2}{ a^4} + \frac{y^2}{ b^4} \biggr] \cdot \biggl(\frac{\rho}{\rho_c}\biggr)\boldsymbol{\hat{f}} \, . </math>
	<math>~=</math>	<math>~ (\bold{u}\cdot \bold{u}) \nabla \biggl(\frac{\rho}{\rho_c}\biggr) \, . </math>

VERY INTERESTING! (29 September 2020)

Exponent q = 4

<math>~\biggl[ (\vec\zeta + 2\vec\Omega) \times \bold{u} \biggr]_{q=4}</math>

<math>~ -~8\Psi_0 \biggl[ \boldsymbol{\hat\imath} \biggl(\frac{x}{a^2} \biggr) + \boldsymbol{\hat\jmath} \biggl(\frac{y}{b^2} \biggr) \biggr] \biggl( \frac{\rho}{\rho_c}\biggr)^{2} \biggl\{ \frac{ 2\Omega_f}{\rho_c} + \frac{8\Psi_0}{\rho_c^2} \biggl(\frac{\rho}{\rho_c}\biggr) \biggl[ \biggl( \frac{\rho}{\rho_c}\biggr) \biggl( \frac{1}{a^2} + \frac{1}{b^2}\biggr) ~- ~ 4 \biggl(\frac{x^2}{ a^4} \biggr) ~- ~ 4 \biggl(\frac{y^2}{ b^4} \biggr) \biggr] \biggr\} \, . </math>

Trial #4

We begin with the,

Euler Equation
written in terms of the Vorticity and
as viewed from a Rotating Reference Frame

<math>\frac{\partial \bold{u}}{\partial t} + (\boldsymbol\zeta+2{\vec\Omega}_f) \times {\bold{u}}= - \frac{1}{\rho} \nabla P - \nabla \biggl[\Phi + \frac{1}{2}u^2 - \frac{1}{2}|{\vec{\Omega}}_f \times \vec{x}|^2 \biggr]</math> .

Next, we rewrite this expression to incorporate the following three realizations:

For a barotropic fluid, the term involving the pressure gradient can be replaced with a term involving the enthalpy via the relation, <math>~\nabla H = \nabla P/\rho</math>.
The expression for the centrifugal potential can be rewritten as, <math>~\tfrac{1}{2}|\vec\Omega_f \times \vec{x}|^2 = \tfrac{1}{2}\Omega_f^2 (x^2 + y^2)</math>.
In steady state, <math>~\partial \bold{u}/\partial t = 0</math>.

This means that,

<math>~ (\boldsymbol\zeta+2{\vec\Omega}_f) \times {\bold{u}}</math>

<math>- \nabla \biggl[H + \Phi_\mathrm{grav} + \frac{1}{2}u^2 - \frac{1}{2}\Omega_f^2 (x^2 + y^2) \biggr] \, .</math>

If the term on the left-hand-side of this equation can be expressed in terms of the gradient of a scalar function, then it can be readily grouped with all the other terms on the right-hand-side, which already are in the gradient form.

Building on the insight that we have gained from the above examination of systems for which the exponent, q = 3, let's change the <math>~\tfrac{1}{2}\nabla u^2</math> term on the RHS to <math>~\tfrac{1}{2}\nabla (\rho u)^2</math> then reexamine the LHS.

<math>~ [ (\boldsymbol\zeta+2{\vec\Omega}_f) \times {\bold{u}} ]_{q=3}</math>	<math>~=</math>	<math>- \nabla \biggl[H + \Phi_\mathrm{grav} - \frac{1}{2}\Omega_f^2 (x^2 + y^2) \biggr] - \nabla \biggl[\frac{1}{2}u^2\biggr]</math>
	<math>~=</math>	<math> - \nabla \biggl[H + \Phi_\mathrm{grav} - \frac{1}{2}\Omega_f^2 (x^2 + y^2) \biggr] - \biggl(\frac{\rho}{\rho_c}\biggr)^{-2} \biggl\{ \biggl(\frac{\rho}{\rho_c}\biggr)^{2} \nabla \biggl[\frac{1}{2}u^2\biggr]\biggr\} </math>
	<math>~=</math>	<math> - \nabla \biggl[H + \Phi_\mathrm{grav} - \frac{1}{2}\Omega_f^2 (x^2 + y^2) \biggr] - \biggl(\frac{\rho}{\rho_c}\biggr)^{-2} \biggl\{ \nabla \biggl[\frac{1}{2}\biggl(\frac{\rho}{\rho_c}\biggr)^{2} u^2\biggr] - \frac{1}{2}u^2\nabla \biggl[\biggl(\frac{\rho}{\rho_c}\biggr)^{2}\biggr] \biggr\} </math>
<math>~\Rightarrow~~~ [ (\boldsymbol\zeta+2{\vec\Omega}_f) \times {\bold{u}} ]_{q=3} - \nabla F_B</math>	<math>~=</math>	<math> - \nabla \biggl[H + \Phi_\mathrm{grav} + F_B - \frac{1}{2}\Omega_f^2 (x^2 + y^2) \biggr] - \biggl(\frac{\rho}{\rho_c}\biggr)^{-2} \biggl\{ \nabla \biggl[\frac{1}{2}\biggl(\frac{\rho}{\rho_c}\biggr)^{2} u^2\biggr] - \frac{1}{2}u^2\nabla \biggl[\biggl(\frac{\rho}{\rho_c}\biggr)^{2}\biggr] \biggr\} \, . </math>

Now, rewriting the LHS gives,

<math>~[ (\boldsymbol\zeta+2{\vec\Omega}_f) \times {\bold{u}} ]_{q=3} - \nabla F_B</math>	<math>~=</math>	<math>~-~ 3\Psi_0 \biggl\{ \frac{ 2\Omega_f}{\rho_c} - \frac{12\Psi_0}{\rho_c^2} \biggl[ \frac{x^2}{ a^4} + \frac{y^2}{ b^4} \biggr] + \frac{6\Psi_0}{\rho_c^2} \biggl( \frac{1}{a^2} + \frac{1}{b^2}\biggr) \biggl( \frac{\rho}{\rho_c}\biggr) \biggr\} \cdot 2\biggl( \frac{\rho}{\rho_c}\biggr)\boldsymbol{\hat{f}} </math>
		<math>~ +~\biggl\{ 2D_1 \Psi_0^{2/3} + 3D_2\Psi_0 \biggl( \frac{\rho}{\rho_c}\biggr) \biggr\} \cdot 2\biggl( \frac{\rho}{\rho_c}\biggr)\boldsymbol{\hat{f}} </math>
	<math>~=</math>	<math>~3\Psi_0 \biggl\{ \frac{12\Psi_0}{\rho_c^2} \biggl[ \frac{x^2}{ a^4} + \frac{y^2}{ b^4} \biggr] \biggr\} \cdot 2\biggl( \frac{\rho}{\rho_c}\biggr)\boldsymbol{\hat{f}} </math>
	<math>~=</math>	<math>~-~ \frac{36\Psi_0^2}{\rho_c^2} \biggl[ \frac{x^2}{ a^4} + \frac{y^2}{ b^4} \biggr] \cdot \nabla \biggl(\frac{\rho}{\rho_c}\biggr)^2 \, , </math>

where we have set,

<math>~D_1 = \frac{3\Psi_0^{1 / 3} \Omega_f}{\rho_c}</math>

and,

Notice that when the exponent, <math>~q=3</math>, we have,

<math>~[ \bold{u}\cdot \bold{u} ]_{q=3} </math>	<math>~=</math>	<math>~ \frac{1}{\rho^2} \biggl\{ \biggl[ \biggl( \frac{\rho}{\rho_c}\biggr)^{q-1} \biggl(\frac{2q \Psi_0 y}{b^2} \biggr) \biggr]^2 + \biggl[ \biggl( \frac{\rho}{\rho_c}\biggr)^{q-1} \biggl(\frac{2q \Psi_0x}{a^2} \biggr) \biggr]^2 \biggr\}_{q=3} </math>
	<math>~=</math>	<math>~ \frac{1}{\rho^2} \biggl\{ \biggl[ \biggl( \frac{\rho}{\rho_c}\biggr)^{2} \biggl(\frac{6 \Psi_0 y}{b^2} \biggr) \biggr]^2 + \biggl[ \biggl( \frac{\rho}{\rho_c}\biggr)^{2} \biggl(\frac{6 \Psi_0x}{a^2} \biggr) \biggr]^2 \biggr\} </math>
	<math>~=</math>	<math>~\frac{36\Psi_0^2}{\rho_c^2} \biggl( \frac{\rho}{\rho_c}\biggr)^2 \biggl[ \frac{x^2}{a^4} + \frac{ y^2}{b^4} \biggr] \, . </math>

Hence,

<math>~[ (\boldsymbol\zeta+2{\vec\Omega}_f) \times {\bold{u}} ]_{q=3} - \nabla F_B</math>

<math>~-~ \biggl( \frac{\rho}{\rho_c}\biggr)^{-2} u^2 \cdot \nabla \biggl(\frac{\rho}{\rho_c}\biggr)^2 \, . </math>

Trial#5

Let's return to the above-mentioned,

Eulerian Representation
of the Euler Equation
as viewed from a Rotating Reference Frame

<math>\frac{\partial \bold{u}}{\partial t} + (\bold{u}\cdot \nabla) \bold{u} = - \frac{1}{\rho} \nabla P - \nabla \Phi_\mathrm{grav}

- {\vec{\Omega}}_f \times ({\vec{\Omega}}_f \times \vec{x}) - 2{\vec{\Omega}}_f \times \bold{u} \, .</math>

In steady state, this can be rewritten as,

<math>~(\bold{u}\cdot \nabla) \bold{u} + 2{\vec{\Omega}}_f \times \bold{u}</math>

<math>~ - \nabla \biggl[H + \Phi_\mathrm{grav} - \frac{1}{2}|{\vec{\Omega}}_f \times \vec{x}|^2 \biggr] \, . </math>

Let's focus on the left-hand-side, which is expressed entirely in terms of the rotating-frame velocity, <math>~\bold{u}</math>, and the (constant) angular frequency of rotation of the coordinate frame, <math>~\Omega_f</math>. Rewriting the LHS, we have,

LHS: <math>~2{\vec{\Omega}}_f \times \bold{u} + (\bold{u}\cdot \nabla) \bold{u}</math>	<math>~=</math>	<math>~ 2\Omega_f \biggl[ \boldsymbol{\hat\jmath} u_x - \boldsymbol{\hat\imath} u_y\biggr] + \biggl[ u_x \frac{\partial}{\partial x} + u_y \frac{\partial}{\partial y} \biggr] \biggl[\boldsymbol{\hat\imath} u_x + \boldsymbol{\hat\jmath} u_y \biggr] </math>
	<math>~=</math>	<math>~ 2\Omega_f \biggl[ \boldsymbol{\hat\jmath} u_x - \boldsymbol{\hat\imath} u_y\biggr] + \boldsymbol{\hat\imath}\biggl[ u_x \frac{\partial u_x}{\partial x} + u_y \frac{\partial u_x}{\partial y} \biggr] + \boldsymbol{\hat\jmath} \biggl[ u_x \frac{\partial u_y}{\partial x} + u_y \frac{\partial u_y}{\partial y} \biggr] </math>
	<math>~=</math>	<math>~ 2\Omega_f \biggl[ \boldsymbol{\hat\jmath} u_x - \boldsymbol{\hat\imath} u_y\biggr] + \biggl[\boldsymbol{\hat\imath} u_x \frac{\partial u_x}{\partial x} + \boldsymbol{\hat\jmath} u_y \frac{\partial u_y}{\partial y} \biggr] + \biggl[ \boldsymbol{\hat\imath}u_y \frac{\partial u_x}{\partial y} + \boldsymbol{\hat\jmath} u_x \frac{\partial u_y}{\partial x} \biggr] </math>
	<math>~=</math>	<math>~ 2\Omega_f \biggl[ \boldsymbol{\hat\jmath} u_x - \boldsymbol{\hat\imath} u_y\biggr] + \frac{1}{2}\biggl\{ \biggl[\boldsymbol{\hat\imath} \frac{\partial u_x^2}{\partial x} + \boldsymbol{\hat\jmath} \frac{\partial u_y^2}{\partial y} \biggr] + \biggl[\boldsymbol{\hat\imath} \frac{\partial u_y^2}{\partial x} + \boldsymbol{\hat\jmath} \frac{\partial u_x^2}{\partial y} \biggr] \biggr\} + \biggl[ \boldsymbol{\hat\imath}u_y \frac{\partial u_x}{\partial y} + \boldsymbol{\hat\jmath} u_x \frac{\partial u_y}{\partial x} \biggr] - \frac{1}{2}\biggl[\boldsymbol{\hat\imath} \frac{\partial u_y^2}{\partial x} + \boldsymbol{\hat\jmath} \frac{\partial u_x^2}{\partial y} \biggr] </math>
	<math>~=</math>	<math>~ 2\Omega_f \biggl[ \boldsymbol{\hat\jmath} u_x - \boldsymbol{\hat\imath} u_y\biggr] + \frac{1}{2}\nabla u^2 + \boldsymbol{\hat\imath} u_y \biggl[ \frac{\partial u_x}{\partial y} - \frac{\partial u_y}{\partial x} \biggr] + \boldsymbol{\hat\jmath} u_x \biggl[ \frac{\partial u_y}{\partial x} - \frac{\partial u_x}{\partial y} \biggr] \, . </math>

Next, to the extent possible, let's express the LHS in terms of the dimensionless mass density,

<math>~\sigma \equiv \frac{\rho}{\rho_c}</math>

<math>~1 - \biggl( \frac{x^2}{a^2} + \frac{y^2}{b^2}\biggr) \, .</math>

We will assume that the stream-function,

<math>~\Psi_0 \sigma^q \, ,</math>

in which case,

<math>~\rho_c \sigma \bold{u}</math>	<math>~=</math>	<math>~\boldsymbol{\hat\imath} \frac{\partial\Psi}{\partial y} - \boldsymbol{\hat\jmath} \frac{\partial\Psi}{\partial x}</math>
<math>~\Rightarrow~~~ \rho_c \bold{u}</math>	<math>~=</math>	<math>~\frac{q}{(q-1)} \biggl[ \boldsymbol{\hat\imath} \frac{\partial\sigma^{q-1}}{\partial y} - \boldsymbol{\hat\jmath} \frac{\partial\sigma^{q-1}}{\partial x} \biggr] \, .</math>

That is,

<math>~\frac{q}{(q-1)} \biggl( \frac{\partial\sigma^{q-1}}{\partial y} \biggr) </math>

and,

<math>~- \frac{q}{(q-1)} \biggl( \frac{\partial\sigma^{q-1}}{\partial x} \biggr) \, .</math>

The first term on the LHS becomes,

<math>~ 2\Omega_f \biggl[ \boldsymbol{\hat\jmath} u_x - \boldsymbol{\hat\imath} u_y\biggr] </math>	<math>~=</math>	<math>~ \frac{2\Omega_f}{\rho_c} \cdot \frac{q}{q-1} \biggl[ \boldsymbol{\hat\jmath} \biggl( \frac{\partial\sigma^{q-1}}{\partial y} \biggr) + \boldsymbol{\hat\imath} \biggl( \frac{\partial\sigma^{q-1}}{\partial x} \biggr)\biggr] </math>
	<math>~=</math>	<math>~ \frac{2\Omega_f}{\rho_c} \cdot \frac{q}{q-1} \nabla (\sigma^{q-1}) \, . </math>

The third term on the LHS becomes,

<math>~ \boldsymbol{\hat\imath} u_y \biggl[ \frac{\partial u_x}{\partial y} - \frac{\partial u_y}{\partial x} \biggr] </math>	<math>~=</math>	<math>~ \boldsymbol{\hat\imath} \biggl( \frac{u_y}{\rho_c} \biggr) \biggl\{ \frac{\partial }{\partial y} \biggl[ \frac{q}{(q-1)} \biggl( \frac{\partial\sigma^{q-1}}{\partial y} \biggr) \biggr] + \frac{\partial }{\partial x} \biggl[ \frac{q}{(q-1)} \biggl( \frac{\partial\sigma^{q-1}}{\partial x} \biggr) \biggr] \biggr\} </math>
	<math>~=</math>	<math>~ \frac{q}{(q-1)} \cdot \boldsymbol{\hat\imath} \biggl( \frac{u_y}{\rho_c} \biggr) \biggl\{ \nabla^2 \sigma^{q-1} \biggr\} </math>
	<math>~=</math>	<math>~ -~\biggl[ \frac{q}{\rho_c(q-1)}\biggr]^2 \biggl\{ \nabla^2 \sigma^{q-1} \cdot \boldsymbol{\hat\imath} \biggl( \frac{\partial\sigma^{q-1}}{\partial x} \biggr) \biggr\} </math>

Similarly, the fourth term on the LHS becomes,

<math>~ \boldsymbol{\hat\jmath} u_x \biggl[ \frac{\partial u_y}{\partial x} - \frac{\partial u_x}{\partial y} \biggr] </math>	<math>~=</math>	<math>~ -~\boldsymbol{\hat\jmath} \biggl( \frac{u_x}{\rho_c} \biggr) \biggl\{ \frac{\partial }{\partial x} \biggl[ \frac{q}{(q-1)} \biggl( \frac{\partial\sigma^{q-1}}{\partial x} \biggr) \biggr] + \frac{\partial }{\partial y} \biggl[ \frac{q}{(q-1)} \biggl( \frac{\partial\sigma^{q-1}}{\partial y} \biggr) \biggr] \biggr\} </math>
	<math>~=</math>	<math>~ -~ \frac{q}{(q-1)} \cdot \boldsymbol{\hat\jmath} \biggl( \frac{u_x}{\rho_c} \biggr) \biggl\{ \nabla^2 \sigma^{q-1} \biggr\} </math>
	<math>~=</math>	<math>~ -~ \biggl[ \frac{q}{\rho_c(q-1)}\biggr]^2\biggl\{ \nabla^2 \sigma^{q-1} \cdot \boldsymbol{\hat\jmath} \biggl( \frac{\partial \sigma^{q-1}}{\partial y} \biggr) \biggr\} \, . </math>

Note that,

<math>~\frac{\partial \sigma^{q-1}}{\partial x_i}</math>	<math>~=</math>	<math>~ (q-1)\sigma^{q-2} \frac{\partial \sigma}{\partial x_i} = - (q-1)\sigma^{q-2} \biggl( \frac{2x_i}{a_i^2} \biggr) </math>
<math>~\Rightarrow ~~~ (\rho_c u)^2 = \rho_c^2(u_x^2 + u_y^2)</math>	<math>~=</math>	<math>~ \frac{q^2}{(q-1)^2}\biggl\{ \biggl[ - (q-1)\sigma^{q-2} \biggl( \frac{2x}{a^2} \biggr) \biggr]^2 + \biggl[ - (q-1)\sigma^{q-2} \biggl( \frac{2y}{b^2} \biggr) \biggr]^2 \biggr\} </math>
	<math>~=</math>	<math>~ 4q^2 \sigma^{2(q-2)} \biggl[\frac{x^2}{a^4} +\frac{y^2}{b^4} \biggr] \, . </math>

Hence,

<math>~\nabla^2 \sigma^{q-1}</math>	<math>~=</math>	<math>~ \frac{\partial}{\partial x} \biggl[ \frac{\partial \sigma^{q-1}}{\partial x}\biggr] + \frac{\partial}{\partial y} \biggl[ \frac{\partial \sigma^{q-1}}{\partial y}\biggr] </math>
	<math>~=</math>	<math>~ - (q-1)\biggl\{ \frac{\partial}{\partial x} \biggl[ \sigma^{q-2} \biggl( \frac{2x}{a^2} \biggr)\biggr] + \frac{\partial}{\partial y} \biggl[ \sigma^{q-2} \biggl( \frac{2y}{b^2} \biggr)\biggr] \biggr\} </math>
	<math>~=</math>	<math>~ - (q-1)\biggl\{ \sigma^{q-2} \frac{\partial}{\partial x} \biggl[ \biggl( \frac{2x}{a^2} \biggr)\biggr] + \biggl( \frac{2x}{a^2} \biggr)\frac{\partial}{\partial x} \biggl[ \sigma^{q-2} \biggr] + \sigma^{q-2}\frac{\partial}{\partial y} \biggl[ \biggl( \frac{2y}{b^2} \biggr)\biggr] + \biggl( \frac{2y}{b^2} \biggr)\frac{\partial}{\partial y} \biggl[ \sigma^{q-2} \biggr] \biggr\} </math>
	<math>~=</math>	<math>~ - (q-1)\biggl\{ \biggl[ \frac{2}{a^2} + \frac{2}{b^2} \biggr] \sigma^{q-2} + \frac{2x(q-2)}{a^2} \biggl[ \sigma^{q-3} \frac{\partial \sigma}{\partial x}\biggr] + \frac{2y(q-2)}{b^2} \biggl[ \sigma^{q-3} \frac{\partial \sigma}{\partial y}\biggr] \biggr\} </math>
	<math>~=</math>	<math>~ - (q-1) \biggl\{ \biggl[ \frac{2}{a^2} + \frac{2}{b^2} \biggr] \sigma^{q-2} - 2 (q-2)\sigma^{q-3} \biggl[ \frac{2x^2}{a^4} +\frac{2y^2}{b^4} \biggr] \biggr\} \, . </math>

So, when they are added together, the third and fourth terms give,

<math>~ \boldsymbol{\hat\imath} u_y \biggl[ \frac{\partial u_x}{\partial y} - \frac{\partial u_y}{\partial x} \biggr] + \boldsymbol{\hat\jmath} u_x \biggl[ \frac{\partial u_y}{\partial x} - \frac{\partial u_x}{\partial y} \biggr] </math>	<math>~=</math>	<math>~ -~\biggl[ \frac{q}{\rho_c(q-1)}\biggr]^2 \nabla^2 \sigma^{q-1} \biggl\{ \boldsymbol{\hat\imath} \biggl( \frac{\partial\sigma^{q-1}}{\partial x} \biggr) +~ \boldsymbol{\hat\jmath} \biggl( \frac{\partial \sigma^{q-1}}{\partial y} \biggr) \biggr\} </math>
	<math>~=</math>	<math>~ \biggl[ \frac{2q^2}{\rho_c^2(q-1)}\biggr] \biggl\{ \biggl[ \frac{1}{a^2} + \frac{1}{b^2} \biggr] \sigma^{q-2} - 2 (q-2)\sigma^{q-3} \biggl[ \frac{x^2}{a^4} +\frac{y^2}{b^4} \biggr] \biggr\} \biggl\{ \boldsymbol{\hat\imath} \biggl( \frac{\partial\sigma^{q-1}}{\partial x} \biggr) +~ \boldsymbol{\hat\jmath} \biggl( \frac{\partial \sigma^{q-1}}{\partial y} \biggr) \biggr\} </math>
	<math>~=</math>	<math>~ \biggl[ \frac{2q^2}{\rho_c^2(q-1)}\biggr]\biggl[ \frac{1}{a^2} + \frac{1}{b^2} \biggr] \biggl\{ \sigma^{q-2} \biggr\} \biggl\{ \boldsymbol{\hat\imath} \biggl( \frac{\partial\sigma^{q-1}}{\partial x} \biggr) +~ \boldsymbol{\hat\jmath} \biggl( \frac{\partial \sigma^{q-1}}{\partial y} \biggr) \biggr\} </math>
		<math>~ -~\biggl[ \frac{4(q-2)q^2}{\rho_c^2(q-1)}\biggr] \biggl\{ \biggl[\frac{(\rho_c u)^2}{4q^2}\biggr] \sigma^{1-q} \biggr\}\biggl\{ \boldsymbol{\hat\imath} \biggl( \frac{\partial\sigma^{q-1}}{\partial x} \biggr) +~ \boldsymbol{\hat\jmath} \biggl( \frac{\partial \sigma^{q-1}}{\partial y} \biggr) \biggr\} </math>
	<math>~=</math>	<math>~ \frac{2q^2}{\rho_c^2(2q-3)} \biggl[ \frac{1}{a^2} + \frac{1}{b^2} \biggr] \nabla \sigma^{2q-3} -~\frac{(q-2)(\rho_c u)^2}{\rho_c^2} \nabla\ln(\sigma) \, . </math>

Hence,

LHS: <math>~2{\vec{\Omega}}_f \times \bold{u} + (\bold{u}\cdot \nabla) \bold{u}</math>

<math>~ \frac{2q\Omega_f}{\rho_c(q-1)} \nabla (\sigma^{q-1}) + \frac{1}{2}\nabla u^2 + \frac{2q^2}{\rho_c^2(2q-3)} \biggl[ \frac{1}{a^2} + \frac{1}{b^2} \biggr] \nabla \sigma^{2q-3} -~(q-2)u^2 \nabla\ln(\sigma) \, . </math>

Note that,

<math>~ \frac{1}{2}\nabla u^2 -~(q-2)u^2 \nabla\ln(\sigma) </math>	<math>~=</math>	<math>~ \frac{1}{2} \biggl\{ \nabla u^2 -~u^2 \nabla\ln[\sigma^{2(q-2)}] \biggr\} </math>
	<math>~=</math>	<math>~ \frac{u^2}{2} \biggl\{ \nabla \ln [u^2] -~ \nabla\ln[\sigma^{2(q-2)}] \biggr\} </math>
	<math>~=</math>	<math>~ \frac{u^2}{2} \biggl\{ \nabla \ln \biggl[ \frac{u^2}{ \sigma^{2(q-2)} } \biggr] \biggr\} </math>
	<math>~=</math>	<math>~ \frac{u^2}{2} \biggl\{ \nabla \ln \biggl[ \biggl( \frac{2q}{ \rho_c } \biggr)^2 \biggl( \frac{x^2}{a^4} + \frac{y^2}{b^4} \biggr)\biggr] \biggr\} </math>
	<math>~=</math>	<math>~ \frac{2q^2}{\rho_c^2} \sigma^{2(q-2)} \biggl[\frac{x^2}{a^4} +\frac{y^2}{b^4} \biggr] \biggl\{ \nabla \ln \biggl[ \biggl( \frac{2q}{ \rho_c } \biggr)^2 \biggl( \frac{x^2}{a^4} + \frac{y^2}{b^4} \biggr)\biggr] \biggr\} </math>
	<math>~=</math>	<math>~ \frac{\sigma^{2(q-2)}}{2} \biggl\{ \nabla \biggl[ \biggl( \frac{2q}{ \rho_c } \biggr)^2 \biggl( \frac{x^2}{a^4} + \frac{y^2}{b^4} \biggr)\biggr] \biggr\} \, . </math>

Or,

<math>~ \frac{1}{2}\nabla u^2 +~(2-q)u^2 \nabla\ln(\sigma) </math>	<math>~=</math>	<math>~ \frac{1}{2} \biggl\{ \nabla u^2 +~u^2 \nabla\ln[\sigma^{2(2-q)}] \biggr\} </math>
	<math>~=</math>	<math>~ \biggl( \frac{1}{2}\biggr)\sigma^{-2(2-q)} \biggl\{ \sigma^{2(2-q)}\nabla u^2 +~u^2 \nabla[\sigma^{2(2-q)}] \biggr\} </math>
	<math>~=</math>	<math>~ \biggl( \frac{1}{2}\biggr)\sigma^{2(q-2)} \biggl\{ \nabla[u^2 \sigma^{2(2-q)}] \biggr\} \, . </math>

Exponent q = 2

LHS: <math>~[2{\vec{\Omega}}_f \times \bold{u} + (\bold{u}\cdot \nabla) \bold{u}]_{q=2}</math>

<math>~ \frac{4\Omega_f}{\rho_c} \nabla \sigma + \frac{1}{2}\nabla u^2 + \frac{8}{\rho_c^2} \biggl[ \frac{1}{a^2} + \frac{1}{b^2} \biggr] \nabla \sigma \, . </math>

Exponent q = 3

LHS: <math>~[2{\vec{\Omega}}_f \times \bold{u} + (\bold{u}\cdot \nabla) \bold{u} ]_{q=3}</math>	<math>~=</math>	<math>~ \frac{3\Omega_f}{\rho_c} \nabla \sigma^{2} + \frac{1}{2}\nabla u^2 + \frac{6}{\rho_c^2} \biggl[ \frac{1}{a^2} + \frac{1}{b^2} \biggr] \nabla \sigma^{3} -~u^2 \nabla\ln(\sigma) \, . </math>
	<math>~=</math>	<math>~ \frac{3\Omega_f}{\rho_c} \nabla \sigma^{2} + \frac{6}{\rho_c^2} \biggl[ \frac{1}{a^2} + \frac{1}{b^2} \biggr] \nabla \sigma^{3} + u^2\nabla \ln\biggl(\frac{u}{\sigma}\biggr) \, . </math>

Trial #6

LHS: <math>~2{\vec{\Omega}}_f \times \bold{u} + (\bold{u}\cdot \nabla) \bold{u}</math>

<math>~ 2\Omega_f \biggl[ \boldsymbol{\hat\jmath} u_x - \boldsymbol{\hat\imath} u_y\biggr] + \frac{1}{2}\nabla u^2 + \boldsymbol{\hat\imath} u_y \biggl[ \frac{\partial u_x}{\partial y} - \frac{\partial u_y}{\partial x} \biggr] + \boldsymbol{\hat\jmath} u_x \biggl[ \frac{\partial u_y}{\partial x} - \frac{\partial u_x}{\partial y} \biggr] \, . </math>

We will assume that the stream-function,

<math>~\Psi</math>	<math>~=</math>	<math>~\Psi_0 (\alpha - \sigma )^q </math>
<math>~\Rightarrow ~~~ \frac{\partial\Psi}{\partial x_i}</math>	<math>~=</math>	<math>~-~\Psi_0 q(\alpha - \sigma )^{q-1} \frac{\partial\sigma}{\partial x_i} \, ,</math>

in which case,

<math>~\rho_c \sigma \bold{u}</math>	<math>~=</math>	<math>~\boldsymbol{\hat\imath} \frac{\partial\Psi}{\partial y} - \boldsymbol{\hat\jmath} \frac{\partial\Psi}{\partial x}</math>
	<math>~=</math>	<math>~ -~\boldsymbol{\hat\imath} \Psi_0 q(\alpha - \sigma )^{q-1} \frac{\partial\sigma}{\partial y} + \boldsymbol{\hat\jmath}~\Psi_0 q(\alpha - \sigma )^{q-1} \frac{\partial\sigma}{\partial x} </math>
<math>~\Rightarrow~~~ \rho_c \bold{u}</math>	<math>~=</math>	<math>~ -~\boldsymbol{\hat\imath} \biggl[ \frac{\Psi_0 q(\alpha - \sigma )^{q-1} }{\sigma} \biggr]\frac{\partial\sigma}{\partial y} + \boldsymbol{\hat\jmath}~\biggl[ \frac{\Psi_0 q(\alpha - \sigma )^{q-1}}{\sigma}\biggr] \frac{\partial\sigma}{\partial x} \, .</math>

That is,

<math>~-~\biggl[ \frac{\Psi_0 q(\alpha - \sigma )^{q-1} }{\sigma} \biggr]\frac{2y}{b^2} </math>

and,

<math>~ \biggl[ \frac{\Psi_0 q(\alpha - \sigma )^{q-1}}{\sigma}\biggr] \frac{2x}{a^2} \, .</math>

The first term on the LHS becomes,

<math>~ 2\Omega_f \biggl[ \boldsymbol{\hat\jmath} u_x - \boldsymbol{\hat\imath} u_y\biggr] </math>	<math>~=</math>	<math>~ \frac{2\Omega_f}{\rho_c} \biggl[- \frac{\Psi_0 q(\alpha - \sigma )^{q-1}}{\sigma}\biggr] \biggl[\boldsymbol{\hat\jmath} \biggl( \frac{2y}{b^2} \biggr) + \boldsymbol{\hat\imath} \biggl( \frac{2x}{a^2} \biggr)\biggr] </math>
	<math>~=</math>	<math>~ + \frac{2\Omega_f}{\rho_c} \biggl[ \frac{\Psi_0 q(\alpha - \sigma )^{q-1}}{\sigma}\biggr] \nabla \sigma \, . </math>

Trial #7

Uncluttered Setup

Let's simply look at the vortensity expression as defined in Part II, above, namely,

<math>~\equiv</math>

<math>~\frac{(2\Omega_f + \zeta_z)}{\rho_c\sigma} \, ,</math>

and recognize that we are ultimately interested in the function, <math>~F_B(\Psi)</math>, defined such that,

We start with the expression for the z-component of the vorticity,

<math>~\zeta_z </math>	<math>~=</math>	<math>~ -~\frac{\partial }{\partial x}\biggl[\frac{1}{\rho} \frac{\partial \Psi}{\partial x} \biggr] - \frac{\partial }{\partial y} \biggl[\frac{1}{\rho} \frac{\partial \Psi}{\partial y} \biggr] </math>
	<math>~=</math>	<math>~ -~\biggl\{ \frac{\partial }{\partial x}\biggl[ \biggl( \frac{\Psi_0}{\rho_c \sigma} \biggr) \frac{\partial \sigma^q}{\partial x} \biggr] + \frac{\partial }{\partial y} \biggl[ \biggl( \frac{\Psi_0}{\rho_c\sigma} \biggr) \frac{\partial \sigma^q}{\partial y} \biggr] \biggr\} \, . </math>

Next, appreciate that,

<math>~\frac{1}{\sigma} \frac{\partial \sigma^q}{\partial x_i}</math>

<math>~ q \sigma^{q-2} \frac{\partial \sigma}{\partial x_i} </math>

<math>~ \biggl( \frac{q}{q-1} \biggr) \frac{\partial \sigma^{q-1}}{\partial x_i} \, . </math>

Hence,

<math>~\zeta_z </math>	<math>~=</math>	<math>~ -~ \frac{q\Psi_0}{\rho_c (q-1) } \biggl\{ \frac{\partial }{\partial x}\biggl[ \frac{\partial \sigma^{q-1}}{\partial x} \biggr] + \frac{\partial }{\partial y} \biggl[ \frac{\partial \sigma^{q-1}}{\partial y} \biggr] \biggr\} </math>
	<math>~=</math>	<math>~ -~ \biggl[ \frac{q\Psi_0}{\rho_c (q-1) } \biggr] \nabla^2\sigma^{q-1} \, . </math>

So, the vortensity is,

<math>~g(\sigma) = \frac{(2\Omega_f + \zeta_z )}{\rho_c\sigma}</math>

<math>~ \biggl(\frac{2\Omega_f}{\rho_c}\biggr)\sigma^{-1} - \biggl[ \frac{q\Psi_0}{\rho_c^2 (q-1) } \biggr] \sigma^{-1} \nabla^2\sigma^{q-1} \, . </math>

Let's switch to the stream-function via the assumed relation,

<math>~\sigma</math>	<math>~=</math>	<math>~\biggl( \frac{\Psi}{\Psi_0}\biggr)^{1/q} \, ;</math>
<math>~\Rightarrow ~~~ g(\Psi) </math>	<math>~=</math>	<math>~ \biggl(\frac{2\Omega_f}{\rho_c}\biggr)\biggl( \frac{\Psi}{\Psi_0}\biggr)^{-1/q} - \biggl[ \frac{q\Psi_0}{\rho_c^2 (q-1) } \biggr] \biggl( \frac{\Psi}{\Psi_0}\biggr)^{-1/q} \nabla^2\biggl( \frac{\Psi}{\Psi_0}\biggr)^{(q-1)/q} \, . </math>

This expression gives the vortensity in what appears to be the desired form — that is, expressed strictly in terms of the stream function, <math>~\Psi</math> — for a wide range of values of the exponent, <math>~q</math>. [CAUTION: the <math>~\nabla^2</math> operator is an exception.] It is not yet (13 October 2020) clear to me how — or if — the second term on the right-hand-side of this expression can be integrated to give <math>~F(\Psi)</math>. But the first term can be obtained from,

<math>~F_{1^\mathrm{st}}(\Psi)</math>	<math>~=</math>	<math>~ -\biggl(\frac{2\Omega_f}{\rho_c}\biggr) \frac{q}{(q-1)} \biggl( \frac{\Psi}{\Psi_0}\biggr)^{(q-1)/q} </math>
<math>~\Rightarrow ~~~ \frac{dF_{1^\mathrm{st}}(\Psi)}{d\Psi}</math>	<math>~=</math>	<math>~ -\biggl(\frac{2\Omega_f}{\rho_c}\biggr) \biggl( \frac{\Psi}{\Psi_0}\biggr)^{-1/q} \, . </math>

T5 Coordinates

Let's evaluate the <math>~\nabla^2</math> operator by expressing it and its argument in terms of T5 Coordinates. Note that,

<math>~\lambda_1</math>

<math>~\equiv</math>

<math>~(x^2 + \kappa^2 y^2)^{1 / 2} \, ,</math>

where, <math>~\kappa \equiv a/b</math>. The specified density distribution is, therefore,

<math>~\sigma \equiv \frac{\rho}{\rho_c}</math>

<math>~1 - \frac{\lambda_1^2}{a^2} \, ,</math>

and the stream-function is,

<math>~\sigma^q = \biggl[ 1 - \frac{\lambda_1^2}{a^2} \biggr]^q \, .</math>

The relevant T5-Coordinate System Laplacian is,

<math>~\nabla^2 \mathfrak{F}</math>

<math>~ \biggl[ \frac{1}{\lambda_1^2 \ell^2} \biggr] \biggl[ \frac{\partial^2 \mathfrak{F}}{\partial \lambda_1^2}\biggr] - \biggl[ \frac{1}{\lambda_1^3 \ell^2} \biggr] \frac{\partial \mathfrak{F}}{\partial \lambda_1} + \biggl[ \frac{(1 + \kappa^2)}{\lambda_1 } \biggr] \frac{\partial \mathfrak{F}}{\partial \lambda_1} \, , </math>

where,

<math>~\ell \equiv (x^2 + \kappa^4 y^2)^{- 1 / 2} \, ,</math>

and in the present context,

<math>~\mathfrak{F}</math>

<math>~ \biggl(\frac{\Psi}{\Psi_0}\biggr)^{(q-1)/q} = \biggl[ 1 - \frac{\lambda_1^2}{a^2} \biggr]^{q-1} \, . </math>

Hence,

<math>~\frac{\partial \mathfrak{F}}{\partial \lambda_1}</math>	<math>~=</math>	<math>~ - \frac{2(q-1)\lambda_1}{a^2} \biggl[ 1 - \frac{\lambda_1^2}{a^2} \biggr]^{q-2} = - \biggl[ \frac{2(q-1)\lambda_1}{a^2}\biggr]\biggl( \frac{\Psi}{\Psi_0} \biggr)^{(q-2)/q} \, ; </math>
<math>~\frac{\partial^2 \mathfrak{F}}{\partial \lambda_1^2}</math>	<math>~=</math>	<math>~ - \frac{2(q-1)}{a^2} \biggl[ 1 - \frac{\lambda_1^2}{a^2} \biggr]^{q-2} + \frac{2(q-1)\lambda_1}{a^2} \biggl[ 1 - \frac{\lambda_1^2}{a^2} \biggr]^{q-3} \biggl[ \frac{2(q-2) \lambda_1 }{ a^2 } \biggr] </math>
	<math>~=</math>	<math>~ - \frac{2(q-1)}{a^2} \biggl( \frac{\Psi}{\Psi_0} \biggr)^{(q-2)/q} + \frac{4(q-1)(q-2)\lambda_1^2}{a^4} \biggl( \frac{\Psi}{\Psi_0} \biggr)^{(q-3)/q} \, . </math>

And,

<math>~\nabla^2 \mathfrak{F}</math>	<math>~=</math>	<math>~ \biggl[ \frac{1}{\lambda_1^2 \ell^2} \biggr] \biggl\{ - \frac{2(q-1)}{a^2} \biggl( \frac{\Psi}{\Psi_0} \biggr)^{(q-2)/q} + \frac{4(q-1)(q-2)\lambda_1^2}{a^4} \biggl( \frac{\Psi}{\Psi_0} \biggr)^{(q-3)/q} \biggr\} </math>
		<math>~ - \biggl[ \frac{1}{\lambda_1^3 \ell^2} \biggr] \biggl\{ - \biggl[ \frac{2(q-1)\lambda_1}{a^2}\biggr]\biggl( \frac{\Psi}{\Psi_0} \biggr)^{(q-2)/q} \biggr\} + \biggl[ \frac{(1 + \kappa^2)}{\lambda_1 } \biggr] \biggl\{ - \biggl[ \frac{2(q-1)\lambda_1}{a^2}\biggr]\biggl( \frac{\Psi}{\Psi_0} \biggr)^{(q-2)/q} \biggr\} </math>
	<math>~=</math>	<math>~ \frac{4(q-1)(q-2)}{a^4\ell^2} \biggl( \frac{\Psi}{\Psi_0} \biggr)^{(q-3)/q} - \biggl[ \frac{2(q-1)(1 + \kappa^2)}{a^2}\biggr]\biggl( \frac{\Psi}{\Psi_0} \biggr)^{(q-2)/q} \, . </math>

All Together

Putting this all together, we obtain,

<math>~ g(\Psi) </math>	<math>~=</math>	<math>~ \biggl(\frac{2\Omega_f}{\rho_c}\biggr)\biggl( \frac{\Psi}{\Psi_0}\biggr)^{-1/q} - \biggl[ \frac{q\Psi_0}{\rho_c^2 (q-1) } \biggr] \biggl( \frac{\Psi}{\Psi_0}\biggr)^{-1/q} \biggl\{ \frac{4(q-1)(q-2)}{a^4\ell^2} \biggl( \frac{\Psi}{\Psi_0} \biggr)^{(q-3)/q} - \biggl[ \frac{2(q-1)(1 + \kappa^2)}{a^2}\biggr]\biggl( \frac{\Psi}{\Psi_0} \biggr)^{(q-2)/q} \biggr\} </math>
	<math>~=</math>	<math>~ \biggl(\frac{2\Omega_f}{\rho_c}\biggr)\biggl( \frac{\Psi}{\Psi_0}\biggr)^{-1/q} + \biggl[ \frac{2q(1 + \kappa^2)\Psi_0}{\rho_c^2 a^2} \biggr] \biggl( \frac{\Psi}{\Psi_0} \biggr)^{(q-3)/q} - \biggl[ \frac{4 q(q-2) \Psi_0}{\rho_c^2 a^4} \biggr] (x^2 + \kappa^4 y^2) \biggl( \frac{\Psi}{\Psi_0} \biggr)^{(q-4)/q} \, . </math>

Note that, for the specific example case of <math>~\kappa^2 = 2</math>,

<math>~\biggl[ \frac{1}{\ell^2} \biggr]_{\kappa^2=2} = (x^2 + 4y^2)</math>

<math>~ x^2 \Lambda = \frac{\Lambda \lambda_2^2}{2} \biggl[ \Lambda - 1 \biggr] = 2\lambda_1^2 \biggl[ \frac{\Lambda}{\Lambda + 1} \biggr] \, . </math>

where,

<math>~\Lambda</math>

<math>~\equiv</math>

<math>~\biggl[ 1 + \frac{4\lambda_1^2}{\lambda_2^2} \biggr]^{1 / 2} \, .</math>

User:Tohline/ThreeDimensionalConfigurations/Challenges

Contents

Challenges Constructing Ellipsoidal-Like Configurations

Riemann S-type Ellipsoids

2^nd-Order TVE Expressions

Ou's (2006) Detailed Force Balance

Hybrid Scheme

Compressible Structures

Part I

Part II

Part III

Equation of Continuity

Related Useful Expressions

Euler Equation

Strategy

Compare

Incompressible

Compressible

Trial #1

Trial #2

Trial #3

Exponent q = 2

Exponent q = 3

Exponent q = 4

Trial #4

Trial#5

Exponent q = 2

Exponent q = 3

Trial #6

Trial #7

Uncluttered Setup

T5 Coordinates

All Together

See Also

Navigation menu

Search

User:Tohline/ThreeDimensionalConfigurations/Challenges

Challenges Constructing Ellipsoidal-Like Configurations

Riemann S-type Ellipsoids

2nd-Order TVE Expressions

Ou's (2006) Detailed Force Balance

Hybrid Scheme

Compressible Structures

Part I

Part II

Part III

Equation of Continuity

Related Useful Expressions

Euler Equation

Strategy

Compare

Incompressible

Compressible

Trial #1

Trial #2

Trial #3

Exponent q = 2

Exponent q = 3

Exponent q = 4

Trial #4

Trial#5

Exponent q = 2

Exponent q = 3

Trial #6

Trial #7

Uncluttered Setup

T5 Coordinates

All Together

See Also

Navigation menu

Search

2^nd-Order TVE Expressions