Graphical Representation of the Energy Surface of a Triaxial Rotor

Questions
,
1

Energy surface - Mathematica

I am working with a triaxial rigid rotor, which is a 3D object with no symmetry axis, so the moments of inertia about the principal axes are all different.
It is well known from textbooks that the energy (rotational) of a triaxial rotor is given by the formula:

E=A_1I_1^2+A_2I_2^2+A_3I_3^2\ ,

where the terms A_k=\frac{1}{2\mathcal{I}_k} are the inverse of the moments of inertia \mathcal{I}_k, and I_k represent the components of the total angular momentum \mathbf{I}.
My problem involves plotting the energy surface of the triaxial rotor, which should be a triaxial ellipsoid, deformed accordingly to the values of the inertia parameters A_k.

I chose to work in spherical coordinates, so a representation of the angular momentum vector using its components will be something like this:

I_1\equiv x_1 =r\sin\theta\cos\varphi \\ I_2\equiv x_2=r\sin\theta\sin\varphi \\ I_3\equiv x_3=r\cos\theta

where r\in[0,I] (can be any number but no larger than the maximum spin value I) , \theta\in[0,\pi] (polar angle) and \varphi\in[0,2\pi] (azimuthal angle).

The following arbitrary values for the inertia factors A3=1 , A_1=6, A_2=3, and the maximum spin I=2 are given. Note that the smallest inertia factor is the one corresponding to the 3-axis (so it means that the 3-axis has the largest moment of inertia). Therefore, choosing x_3 as the quantization axis is the right choice.

Now, for the graphical representation of the energy ellipsoid E, I was thinking of doing a SphericalPlot3D, since I work in spherical coordinates. Fixing r to a value of I=2, and then trying to plot E=f(r=I,\theta,\varphi), should (?) give me an ellipsoid. However, this is not what I obtain. The 3d object has some kind of neck. I don’t understand why that is.

This is the code:

x1[r_, theta_, fi_] := r Sin[theta] Cos[fi]; 
x2[r_, theta_, fi_] := r Sin[theta] Sin[fi]; 
x3[r_, theta_, fi_] := r Cos[theta];

EnSphere[a1_, a2_, a3_, r_, theta_, fi_] := 
  a1*x1[r, theta, fi]^2 + a2*x2[r, theta, fi]^2 + 
   a3*x3[r, theta, fi]^2;

SphericalPlot3D[
 En3[A1, A2, A3, 1, \[Theta], \[Phi]], {\[Theta], 
  0, \[Pi]}, {\[Phi], 0, 2\[Pi]}, PlotStyle -> Blue]

This is what I obtain

ellipsoid2.png
ellipsoid2.png474×263 31.9 KB

This is what I expected to obtain

ellipsoid.jpg
ellipsoid.jpg1414×808 166 KB

Why is this happening? And how can I solve this? Is there a way to obtain the 3D energy surface of the triaxial rigid rotor starting from the values of the inertial parameters A_k?

Thank you in advance!

2

Like this perhaps?

ContourPlot3D[ a1 x^2 + a2 y^2 + a3 z^2 /. {a1 -> 6, a2 -> 3, a3 -> 1}, 
{x, -1, 1}, {y, -1, 1}, {z, -1.3, 1.3}, Contours -> 1.5]
3

Well, you used cartesian coordinates. If I switch to the components x1, x2, and x3 which are expressed in terms of spherical coordinates, then the result I get is different.

ContourPlot3D[
 a1 x1[r, \[Theta], \[Phi]]^2 + a2 x2[r, \[Theta], \[Phi]]^2 + 
   a3 x3[r, \[Theta], \[Phi]]^2 /. {a1 -> 6, a2 -> 3, a3 -> 1}, {r, 0,
   2}, {\[Theta], 0, \[Pi]}, {\[Phi], 0, 2 \[Pi]}, Contours -> 1.5]
4

SphericalPlot3D works as a parametric plot, not as an implicit plot. You must first solve with respect to r and then feed the result to SphericalPlot3D:

Block[{a1 = 1, a2 = 2, a3 = 3},
 sol = Solve[EnSphere[a1, a2, a3, r, \[Theta], \[Phi]] == 1, r];
 rAsFunctionOFThetaFi = r /. sol[[2]];
 SphericalPlot3D[
  Evaluate[rAsFunctionOFThetaFi], {\[Theta], 0, \[Pi]}, {\[Phi], 0, 
   2 \[Pi]}, PlotStyle -> Blue]]
5

Well, I think we have two totally different representations.
The contourplot asks for all momenta {x,y,z} with a given energie (1.6 here , Gianluca looked at energy = 1 )

ContourPlot3D[ 1.6 == 6 x^2 + 3 y^2 + z^2, {x, -1, 1}, {y, -1, 1}, {z, -1.3, 1.3}]

(or

RegionPlot3D[ 0 < 6 x^2 + 3 y^2 + z^2 < 1.6, {x, -1, 1}, {y, -1, 1}, {z, -1, 1}]

)

But the momenta giving this energy don’t have all the length 1. E.g. chose a momentum { .2, .3, z }
and ask, for which z the energy 1.6 is achieved . You will get a solution, but this has not the length 1.

Clear[x, y, z]
vec = {.2, .3, z} /. Solve[1.6 == 6 x^2 + 3 y^2 + z^2 /. x -> .2 /. y -> .3, z][[1, 1]]
vec.vec

If on the other hand you ask for the energies produced by a momentum of length 1 you have that these energies are angle-dependent.

Plot[6 x1[1, Pi/2, \[Phi]]^2 + 3 x2[1, Pi/2, \[Phi]]^2 +   1 x3[1, Pi/2, \[Phi]]^2, {\[Phi], 0, 2 \[Pi]}]
Plot[6 x1[1, th, 0]^2 + 3 x2[1, th, 0]^2 + 1 x3[1, th, 0]^2, {th, 0,  Pi}]

This explains the first figure in blue: it is the energy given by momenta of length 1 pointing in different directions.