How to change coordinates of a differential operator?

Question

I'm doing a basic quantum mechanics problem and am trying to learn how to do it in Mathematica. Any help would be much appreciated.

$\vec{L} = \vec{x} \times \vec{p}$

where $\vec{x}$ has components:

$\begin{align*} x &= r \sin \theta \cos \phi\\ y &= r \sin \theta \sin \phi\\ z &= r \cos \theta \end{align*}$

and $\vec{p}$ has components:

$p_x = - i \hbar \dfrac{\partial}{\partial x}$

etc.

Show that $L_z = -i \hbar \dfrac{\partial}{\partial \phi}$

Answer 1

Here is the Mathematica proof. I'll leave out the prefactor $\hbar/i$ for simplicity. Also, in case this is a homework problem, I decided not to add too many comments to the code. Instead I'll let you figure it out. The basic idea is to do cross products and gradients in spherical coordinates. The calculation shown here actually gives you a way to calculate all angular momentum components (not just the $z$ component).

In terms of notation, I use the spherical coordinates r, θ, φ and call their unit vectors er, eθ, eφ.

Needs["VectorAnalysis`"]

Clear[r, θ, φ]

SetCoordinates[Spherical[r, θ, φ]]

(* ==> Spherical[r, θ, φ] *)

{er, eθ, eφ} = 
 Transpose[JacobianMatrix[{r, θ, φ}]]/
  ScaleFactors[{r, θ, φ}]

$\left( \begin{array}{ccc} \cos (\phi ) \sin (\theta ) & \sin (\theta ) \sin (\phi ) & \cos (\theta ) \\ \cos (\theta ) \cos (\phi ) & \cos (\theta ) \sin (\phi ) & -\sin (\theta ) \\ -\sin (\phi ) & \cos (\phi ) & 0 \\ \end{array} \right)$

l = 
 FullSimplify@
  Cross[r er, 
   Transpose[{er, eθ, eφ}].Grad[f[θ, φ]]]

$\left( \begin{array}{c} \cot (\theta ) (-\cos (\phi )) f^{(0,1)}(\theta ,\phi )-\sin (\phi ) f^{(1,0)}(\theta ,\phi )\\ \cos (\phi ) f^{(1,0)}(\theta ,\phi )-\cot (\theta ) \sin (\phi ) f^{(0,1)}(\theta ,\phi )\\ f^{(0,1)}(\theta ,\phi) \end{array} \right)$

This is a list of Cartesian components (because the Cross product needed to be done in Cartesian coordinates) - the last entry in the list l is the z component of the angular momentum.

Here f is a test function to which I apply the operator. $f^{(0,1)}(\theta ,\phi)$ is the derivative with respect to the azimuthal coordinate.

The spherical components of the angular momentum (which you didn't ask for, but I mentioned in an earlier version of this answer) are obtained as follows:

FullSimplify[{er, eθ , eφ}.l]

$\left(0,-\csc (\theta ) f^{(0,1)}(\theta ,\varphi ),f^{(1,0)}(\theta ,\varphi )\right)$

Edit for version 9

Although the above solution still works if you load the VectorAnalysis package, there are compatibility warnings now because the package has become obsolete with Mathematica version 9. So for this new version it's better to extract the coordinate transformation data in a different way:

{er, eθ, eφ} = 
  CoordinateTransformData[{"Spherical" -> "Cartesian"}, 
   "OrthonormalBasisRotation", {r, θ, φ}];
l = FullSimplify@
  Cross[r er, 
   Transpose[{er, eθ, eφ}].Grad[
     f[θ, φ], {r, θ, φ}, 
     "Spherical"]]

The result is the same as above. Typing the lengthy "OrthonormalBasisRotation" expression to define the orthonormal basis vectors of the spherical coordinate system takes some getting used to. It's almost equally easy to re-derive them from first principles - but I wanted to show where the obsolete JacobianMatrix and ScaleFactors information is hidden now.

Answer 2

I wrote some quite general solution, that uses built in functions for coordinate transforming. I hope that it will be helpful.

\[ScriptCapitalD]Transform[expr_,$\[ScriptCapitalD]Functions_,$\[ScriptCapitalD]Variables_,$\[ScriptCapitalD]NewVariables_,system1_->system2_,assumptions_:True]:=
    Module[{inner,\[ScriptCapitalD]},
    \[ScriptCapitalD][sym_/;MemberQ[$\[ScriptCapitalD]Functions,sym],var_/;MemberQ[$\[ScriptCapitalD]Variables,var]]:=Sum[\[ScriptCapitalD][sym,{inner}]\[ScriptCapitalD][inner,{var}],{inner,$\[ScriptCapitalD]NewVariables}];

\[ScriptCapitalD][expr1_+expr2_,var_]:=\[ScriptCapitalD][expr1,var]+\[ScriptCapitalD][expr2,var];

\[ScriptCapitalD][\[ScriptCapitalD][vars1__]\[ScriptCapitalD][vars2__],var_]:=\[ScriptCapitalD][\[ScriptCapitalD][vars1],var]\[ScriptCapitalD][vars2]+\[ScriptCapitalD][vars1]\[ScriptCapitalD][\[ScriptCapitalD][vars2],var];

\[ScriptCapitalD][\[ScriptCapitalD][sym_/;MemberQ[$\[ScriptCapitalD]Functions,sym],var1_List],var2_]:=Sum[\[ScriptCapitalD][sym,Sort[var1~Join~{inner}]]\[ScriptCapitalD][inner,{var2}],{inner,$\[ScriptCapitalD]NewVariables}];

\[ScriptCapitalD][\[ScriptCapitalD][sym_/;MemberQ[$\[ScriptCapitalD]NewVariables,sym],var1_List],var2_]:=\[ScriptCapitalD][sym,Sort[var1~Join~{var2}]];

\[ScriptCapitalD][Derivative[ders__][f_][vars__]]:=Fold[\[ScriptCapitalD][#1,#2]&,f,Flatten[MapThread[Table[#1,{#2}]&,{{vars},{ders}}]]];

\[ScriptCapitalD][someExpr_]:=someExpr/.Derivative[ders__][f_][vars__]:>\[ScriptCapitalD][Derivative[ders][f][vars]];

coordinateTransformRules=Thread[$\[ScriptCapitalD]NewVariables->CoordinateTransform[ system1->system2,$\[ScriptCapitalD]Variables]];

inverseCoordinateTransformRules=Thread[$\[ScriptCapitalD]Variables->CoordinateTransform[ system2->system1,$\[ScriptCapitalD]NewVariables]];

\[ScriptCapitalD]InverseRules={HoldPattern[\[ScriptCapitalD][sym_/;MemberQ[$\[ScriptCapitalD]NewVariables,sym],vars_]]:>With[{function=sym/.coordinateTransformRules},D[function,Sequence@@vars]],

HoldPattern[\[ScriptCapitalD][sym_/;MemberQ[$\[ScriptCapitalD]Functions,sym],vars_]]:>D[sym@@$\[ScriptCapitalD]NewVariables,Sequence@@vars]};

FullSimplify[Expand[\[ScriptCapitalD][expr]]/.\[ScriptCapitalD]InverseRules/.inverseCoordinateTransformRules,assumptions]]

Some examples:

\[ScriptCapitalD]Transform[D[u[x, y], x, x] + D[u[x, y], y, y], {u}, {x, y}, {r, \[Phi]}, "Cartesian" -> "Polar", {r > 0, 0 < \[Phi] < 2 Pi}]

\[ScriptCapitalD]Transform[D[u[x, y, z], x, x] + D[u[x, y, z], y, y] + D[u[x, y, z], z, z], {u}, {x, y, z}, {r, \[Theta], \[Phi]}, "Cartesian" -> "Spherical", {r > 0, 0 < \[Theta] < Pi, -Pi < \[Phi] < Pi}]

\[ScriptCapitalD]Transform[D[T[A, B], A] + 2 D[T[A, B], A, B] + 1/x D[T[A, B], B, B], {T}, {A,  B}, {\[CapitalPsi], \[CapitalOmega]}, "Cartesian" -> "PlanarParabolic", {\[CapitalPsi] \[Element] Reals, 0 <= \[CapitalOmega], (\[CapitalPsi] != 0 || \[CapitalOmega] != 0)}]