Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I have painstakingly derived the vector-spherical harmonics $\mathbf{V}_{J,\,M}^\ell(\theta, \phi)$, which are the generalization of ordinary spherical harmonics $Y_\ell^m(\theta, \phi)$ to vector fields. But now, I would like to visualize them.

The vector-spherical harmonics takes three integers ($\ell$, $J$, $M$), and yields a 3D vector field on the surface of a sphere ($\theta$, $\phi$). The integers are subject to the constraints: $J\geq0$, $\ell\geq0$, $|J-\ell|\leq 1$, $|M|\leq J$. The function VectorSphericalHarmonicV below generates a 3-component complex vector.

Clear[ϵ]; (*Polarization vector*)
ϵ[λ_] = Switch[λ,
   -1, {1, -I, 0}/Sqrt[2],
    0, {0, 0, 1},
    1, {1, I, 0}/Sqrt[2]   ];

VectorSphericalHarmonicV[ℓ_, J_, M_, θ_, ϕ_] /; 
  J >= 0 && ℓ >= 0 && Abs[J - ℓ] <= 1 && Abs[M] <= J :=
    If[Abs[M - λ] <= ℓ, ClebschGordan[{ℓ, M - λ}, {1, λ}, {J, M}], 0]*
    SphericalHarmonicY[ℓ, M - λ, θ, ϕ]*ϵ[λ], {λ, -1, 1}   ]

VectorSphericalHarmonicV[1, 1, -1, θ, ϕ]
VectorSphericalHarmonicV[1, 1, 0, θ, ϕ]
VectorSphericalHarmonicV[1, 1, 1, θ, ϕ]

Honestly, not sure what the best visual representation is. But, I am thinking a good one is to display a bunch of vectors on the surface of a sphere in the spirit of this one:

enter image description here

A sticking point for me is how do I tell Mathematica to display vectors ONLY an the surface? Also, because these vectors are complex, maybe the vectors can be colored to indicate the phases?

Perhaps you have a better way to represent these functions?

share|improve this question
I don't think this is exactly what you want so I'm not writing it as an answer but this may help: vec = Re[{{ Sin[\[Theta]],Cos[\[Theta]] Sin[\[Phi]],Cos[\[Theta]]Cos[\[Phi]]},({(1 + #[[1]]) Sin[\[Theta]], (1+#[[2]])Cos[\[Theta]] Sin[\[Phi]], (1+#[[3]])Cos[\[Theta]]Cos[\[Phi]]}& /@{VectorSphericalHarmonicV[1,1,1,\[Theta],\[Phi]]})[[1]]}] then define a table of vectors: tab = Flatten[Table[Arrow[vec],{\[Theta],0,2\[Pi],2 \[Pi]/30},{\[Phi],0,\[Pi],\[Pi]/30}]] then plot it: Graphics3D[{Arrowheads[Small], #} & /@ tab] –  Jonathan Shock May 12 '13 at 6:16
it seems you could use: Mollweide maps subject to a colour coding the module? –  chris May 12 '13 at 7:07

2 Answers 2

up vote 25 down vote accepted

Edit: I added more explanations below, because this visualization method is quite different from conventional vector plots

For just this purpose I had at some point invented the following visualization technique. I'll reproduce your definition first. It defines a complex vector field on the surface of a unit sphere.

Clear[\[Epsilon]];(*Polarization vector*)\[Epsilon][\[Lambda]_] = 
 Switch[\[Lambda], -1, {1, -I, 0}/Sqrt[2], 0, {0, 0, 1}, 
  1, {1, I, 0}/Sqrt[2]];

VectorSphericalHarmonicV[\[ScriptL]_, J_, M_, \[Theta]_, \[Phi]_] /; 
  J >= 0 && \[ScriptL] >= 0 && Abs[J - \[ScriptL]] <= 1 && 
   Abs[M] <= J := 
 Sum[If[Abs[M - \[Lambda]] <= \[ScriptL], 
    ClebschGordan[{\[ScriptL], M - \[Lambda]}, {1, \[Lambda]}, {J, 
      M}], 0]*SphericalHarmonicY[\[ScriptL], 
    M - \[Lambda], \[Theta], \[Phi]]*\[Epsilon][\[Lambda]], \
{\[Lambda], -1, 1}]

What makes complex-valued 3D vector fields special?

The question is specifically about an example of a complex vector field in three dimensions, as it's commonly encountered in electromagnetism.

For a real-valued three-dimensional vector field, one often uses arrows attached to a set of points to get a visualization that contains all the information about the field. For a complex-valued field, one can naively do the same thing by plotting only the real part (or the imaginary part) of the vectors. However, by doing this, we lose exactly half of the information that's contained in the complex vector field: out of its six real numbers at every space point, we use only three to make a vector.

The method below does not throw away any information while creating a single 3D representation of the field. It can replace conventional vector plots on a three-dimensional domain when your vector field is complex.

Now comes my main function. It visualizes this field by plotting its polarization ellipses. It's an interesting fact that the field, when multiplied by a common complex phase factor (in electrodynamics that woud be time evolution), locally lies in a plane and sweeps out an ellipse (as a function of that phase factor). I calculate the axes of this ellipse at every point in a discrete grid (here chosen to span part of a spherical surface). The ellipses are plotted as little disks with the appropriate orientation. The color represents the relative phase of the field vectors.

ellipseListPlot[f_, pointList_, scale_: .1, thickness_: .01, 
  colorfunction_: ColorData["Rainbow"]] := 
 Module[{m = N[f @@@ pointList], \[Epsilon]},
  \[Epsilon] = Arg[Map[Dot[#, #] &, m]]/2;
     GeometricTransformation[{{colorfunction[(1 + Cos[#3])/2], 
         Line[{{-scale, 0, 0}, {scale, 0, 0}}], 
         Line[{{0, -scale, 0}, {0, scale, 0}}], 
         Cylinder[{{0, 0, 0}, {0, 0, thickness}}, 
          scale]}, {colorfunction[(1 - Cos[#3])/2], 
         Cylinder[{{0, 0, -thickness}, {0, 0, 0}}, scale]}}, 
          Append[#1, Normalize[Cross[#1[[1]], #1[[2]]]]]], #2}
        ]] &,
     {Map[{Re[#], Im[#]} &, m Exp[-I \[Epsilon]]], 
      pointList, \[Epsilon]}]
    }, Lighting -> "Neutral"


Clear[r, \[Theta], \[Phi]]

grid = N[With[{kr = 5},
     Table[CoordinatesToCartesian[{kr, \[Theta], \[Phi]}, 
       Spherical], {\[Theta], Pi/40, Pi/2, Pi/40}, {\[Phi], 
       0, Pi/2, Pi/20}], 1]]];

 VectorSphericalHarmonicV[2, 1, -1, 
   Sequence @@ Rest[CoordinatesFromCartesian[{##}, Spherical]]] &,
 grid, .5


This type of visualization contains practically all the information about the complex field at the given points. In this example, you see that the polarization is locally linear in some places and circular in others, and accordingly the ellipses degenerate to lines or circles.

The function ellipseListPlot can take an arbitrary list of points as the argument grid, so you can also plot the disks (ellipses) in three dimensions inside the sphere, e.g. This becomes interesting if you add on a radial dependence (spherical Bessel functions, say).

The optional arguments scale and thickness define the overall size of the polarization ellipses. The last optional argument is the ColorFunction for the relative phase.

The upshot of the ellipse representation is that you could imagine a movie (see the example below) where you take the real part of the field vector $\vec{v}$ at every point in space, but only after multiplying it by a phase factor $\exp(i \alpha)$. Now let $\alpha$ vary from $0$ to $2\pi$ and record where the vector $\text{Re}[\vec{v}\,\exp(i \alpha)]$ points. This will describe an elliptical trajectory as a function of $\alpha$, and it is these ellipses that I'm plotting. The ellipses have an orientation in three dimensions, just like arrows, but their axis ratios and size encode additional information.

The color represents the angle with respect to the major axis of each ellipse at which this vector is found when $\alpha = 0$ (I called it the relative phase above).

If you count how many real numbers one needs to uniquely specify the size and orientation of the ellipse at each point, together with the color information, it comes out to be six, just as many as are contained in the original vector field. So this representation preserves all the information about the complex vector field, in contrast to a simple vector plot of the real or imaginary part.

A use case for this ellipse representation would be Figure 3 of this paper.

Edit: Detailed explanation

The ellipses are drawn using AffineTransform, and I explained the specific method in this answer to How to draw an ellipse arc in 3D.

To show more clearly what the ellipses have to do with the actual complex vector field, here is a slightly modified version of the ellipseListPlot function where I just added a black line to each ellipse, pointing in the direction of Re[f], the real part of the first argument:

ellipseListPlot[f_, pointList_, scale_: .1, thickness_: .01, 
  colorfunction_: ColorData["Rainbow"]] := 
 Module[{m = N[f @@@ pointList], \[Epsilon]}, \[Epsilon] = 
   Arg[Map[Dot[#, #] &, m]]/2;
  Graphics3D[{EdgeForm[], MapThread[
       {Thick, Black, Line[{#2, #2 + scale Re[#4]}]},
          colorfunction[(1 + Cos[#3])/2],
          Line[{{-scale, 0, 0}, {scale, 0, 0}}], 
          Line[{{0, -scale, 0}, {0, scale, 0}}],
          Cylinder[{{0, 0, 0}, {0, 0, thickness}}, scale]
          colorfunction[(1 - Cos[#3])/2], 
          Cylinder[{{0, 0, -thickness}, {0, 0, 0}}, scale]
           Append[#1, Normalize[Cross @@ #1]]], #2}]]
       } &, {Map[{Re[#], Im[#]} &, m Exp[-I \[Epsilon]]], 
      pointList, \[Epsilon], m}]}, Lighting -> "Neutral"]]

plots = Table[
     Exp[I α] VectorSphericalHarmonicV[2, 1, -1, 
        Sequence @@ 
         Rest[CoordinatesFromCartesian[{##}, Spherical]]] &, grid, .5],
    ViewPoint -> {1.5, 1, 1}
   {α, Pi/24, 2 Pi, Pi/24}

plots1 = 
  Map[ImageResize[Image[#, ImageResolution -> 300], 350] &, plots];

Export["vsh1.gif", plots1, "DisplayDurations" -> .05, 
 AnimationRepetitions -> Infinity]

movie pointers

The animation was created by plotting not just VectorSphericalHarmonicV but Exp[I α] VectorSphericalHarmonicV and varying $\alpha$ to create a list plots. To get better anti-aliasing, I used a higher ImageResolution and scaled it down again, before exporting as GIF.

Now the rotating black lines are like hands of a clock, and they trace out the ellipses. Equal pointer position with respect to the major axis show up as equal color. There is one unavoidable ambiguity in the colors which happens when the role of the major axis switches from one to the other principal axis of the ellipse (going from "prolate" to "oblate"). In that circular limiting case the color also changes discontinuously (unless you choose a color scheme that is cyclic across that transition).

Mathematical details

To explain the math behind finding the axes of the ellipses, let's call the function in the last movie $\vec{w}\equiv \vec{f} \exp(i \alpha)$. Define its real and imaginary parts as

$$\vec{q} \equiv \text{Re}[\vec{w}],\qquad \vec{p} \equiv \text{Im}[\vec{w}] $$

The pointer in the movie shows $\vec{q}$. To determine at what value of $\alpha$ this pointer equals one of the principal axes, you can look for the extrema of the length of $\vec{q}$. That is, we require

$$ 0 = \frac{d}{d\alpha}\left(\vec{q}\right)^2 $$

Expressing the real part as the sum of $\vec{w}$ and its complex conjugate, this leads to

$$ 0 = \frac{d}{d\alpha}\left(\vec{f}\exp(i\alpha)+\vec{f}^*\exp(-i\alpha)\right)^2 $$

The cross terms are independent of $\alpha$ and therefore the condition simplifies to

$$ 0 = \frac{d}{d\alpha}\left[\left(\vec{f}\exp(i\alpha)\right)^2+\left(\vec{f}^*\exp(-i\alpha)\right)^2\right] $$

Exactly the same condition also applies to the extrema of $\vec{p}$ because its square differs only in the sign of the cross terms (which make no contribution). This means the condition for the extremal phase angle $\alpha$, at which both real and imaginary parts of $\vec{w}$ must be principal axis vectors, becomes

$$ \frac{\left(\vec{f}\right)^2}{\left(\vec{f}^*\right)^2} = \exp(-2i\alpha) $$

Introduce modulus and phase on the left-hand side:

$$ \left(\vec{f}\right)^2 = \left|\left(\vec{f}\right)^2\right| \exp(i \arg(\left(\vec{f}\right)^2)) $$

and the moduli cancel out to yield the solution

$$ \alpha = -\frac{1}{2} \arg(\left(\vec{f}\right)^2) $$

This is the quantity I call $\epsilon$ in the code above. Once you have this angle, the ellipse axes are obtained by extracting the real and imaginary parts of $\vec{w}=\vec{f} \exp(i \alpha)$.

For another, alternative derivation of the principal axes of a complex 3D vector field, see the book by J.F. Nye: Natural Focusing and Fine Structure of Light: Caustics and Wave Dislocations.

share|improve this answer
This is fantastic. In v9, there is the built-in function CoordinateTransform that supersedes CoordinatesToCartesian and (I believe) the VectorAnalysis package as a whole. Soon, I will be multiplying with the BesselJ function, just as you predicted, and will need your help when I do! –  QuantumDot May 12 '13 at 22:53
Jens, have you seen this, by any chance? –  Guess who it is. May 13 '13 at 2:53
@J.M. Thanks for the link. But I'm sad to say that my institution doesn't have access to that journal. Anyway, the math behind this is not too hard - I could add it tomorrow. My favorite reference for this would be Natural Focusing and Fine Structure of Light: Caustics and Wave Dislocations. I haven't been able to find a reference where these ellipses were used for visualization in 3D space, the way I do here. –  Jens May 13 '13 at 3:23
I brought it up since I was about ready to post an answer based on that paper, and then saw that you've already posted something. If my reading of that paper's right, your method's a practical realization of the proposal in that paper. (Sadly I cannot give another +1.) –  Guess who it is. May 13 '13 at 3:27
@QuantumDot Thanks for the bounty! I'll get back to the field line issues in the next day or so... –  Jens May 21 '13 at 15:48

Following this question you can define:

invmollweide[{x_, y_}] := 
With[{theta = ArcSin[y]}, {Pi (x)/(2 Cos[theta]), 
ArcSin[(2 theta + Sin[2 theta])/Pi]}];

fc[phi_] := Block[{theta}, If[Abs[phi] == Pi/2, phi, theta /. 
 FindRoot[2 theta + Sin[2 theta] == Pi Sin[phi], {theta, phi}]]];
cart[{lambda_, phi_}] := 
With[{theta = fc[phi]}, {2/Pi*lambda Cos[theta], Sin[theta]}]


 gr2 = VectorPlot[{-1 - Sin[x]^2 + Sin[3 y] + Cos[y]^2, 
 1 + Sin[2 x] - Cos[y]^2}, {x, -Pi, Pi}, {y, -Pi/2, Pi/2}, 
 VectorScale -> Small, VectorPoints -> 50, 
 VectorColorFunction -> Function[{x, y, vx, vy, n}, 
 ColorData["Heat"][Sin[x + y] Cos[x]^2]],
 AspectRatio -> 1/2, Frame -> False];
 gr2 = gr2 /. Arrow[pts_] :> Arrow[(cart /@ pts)] /. 
 Point[pts_] :> Point[cart[pts]] // 
 Show[#, PlotRange -> {{-2, 2}, {-1, 1}}] &


where the colour code could code the third component of your spin Ylms while the vector codes the first 2?

One could mix this answer and Jen's if the complex phase needs to be coded too?

share|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.