# Operator operation

There are two basic operators r={x,y,z} and p=-I*h*{D[#,x] , D[#,y], D[#,z]} &. The angular momentum operator is defined as $\hat{\textbf L}=\hat{\textbf r}\times\hat{\textbf p}$. So how can I get L from p and r in mathematica? The first thing that comes to mind, l=Cross[r,p], won't work.

P.S. how could I define p in a better, more MMA-styled way?

Assuming that you have decided to call the coordinates x, y, and z throughout, you could do this, building on the definition you already have (I shortened your momentum definition, and in addition, you could also use $\hbar$ instead of h by entering ESChbESC):

r = {x, y, z};
p= -I h D[#,{{x,y,z}}]&
L = Function[{f}, r\[Cross]p[f]];
L[f[x, y, z]]


$$\left\{i h z f^{(0,1,0)}(x,y,z)-i h y f^{(0,0,1)}(x,y,z),\\i h x f^{(0,0,1)}(x,y,z)-i h z f^{(1,0,0)}(x,y,z),\\i h y f^{(1,0,0)}(x,y,z)-i h x f^{(0,1,0)}(x,y,z)\right\}$$

By using Function, the angular momentum operator executes the cross product only when it is provided with the argument of the function. This makes the Cross product into an operator on the space of functions $f(x,y,z)$, with the understanding that the differentiation variables for the momentum are always $x, y, z$.

Now the vector operator itself is less useful than its components, so you may also want this:

{Lx, Ly, Lz} =
Map[Function[{f}, L[f][[#]]] &, {1, 2, 3}];


It allows you to check the commutation relations:

Simplify[
Lx[Ly[ψ[r]]] - Ly[Lx[ψ[r]]] == I h Lz[ψ[r]]]

(* ==> True *)


For additional ideas, e.g. on how to construct linear combinations of operators, see my answer to the question Having the derivative be an operator.

Some more angular momentum calculations

With this starting point for the angular momentum vector operator, we can also define the squared angular momentum L2, and the lowering operator lMinus:

L2 = Function[{f}, Lx[Lx[f]] + Ly[Ly[f]] + Lz[Lz[f]]];

lMinus = Function[{f}, Lx[f] - I Ly[f]];


Now test some properties of the eigenfunctions of angular momentum. The maximal state of an orbital angular momentum ladder is the spherical harmonic $Y_l^l$, which in Cartesian coordinates is proportional to $(x+ i y)^l$. Check that this is an eigenfunction of the z component and square of angular momentum:

First[
Solve[
Simplify[Lz[(x + I y)^l] == eigenvalue (x + I y)^l], eigenvalue]
]
*)

(* ==> {eigenvalue -> h l} *)

First[
Solve[
Simplify[L2[(x + I y)^l] == eigenvalue (x + I y)^l], eigenvalue]
]

(* ==> {eigenvalue -> h^2 l (1 + l)} *)


Indeed, the eigenvalues are as expected.

Next, create a lower z angular momentum by applying the lowering operator to the maximal state:

First[
Solve[
Simplify[Lz[lMinus[(x + I y)^l]] == eigenvalue lMinus[(x + I y)^l]],
eigenvalue]
]

(* ==> {eigenvalue -> h (-1 + l)} *)

First[
Solve[
Simplify[L2[lMinus[(x + I y)^l]] == eigenvalue lMinus[(x + I y)^l]],
eigenvalue]
]

(* ==> {eigenvalue -> h^2 l (1 + l)} *)


The z component has been lowered, but the state is still an eigenstate of L2 with the same eigenvalue as before.

Remark on many-particle systems

If you want to use the angular momentum operator on a system in which there is not just one particle with coordinates {x, y, z}, then it becomes necessary to generalize your definitions so that the operator knows which of the particles is supposed to be acted on. That can be done by giving all the above definitions an additional set of arguments to name the variables with respect to which the momentum operator does its differentiations:

Clear[p, L];

p[{x_, y_, z_}] = Function[{f}, -I h D[f, {{x, y, z}}]];
L[{x_, y_, z_}] = Function[{f}, {x, y, z}\[Cross]p[{x, y, z}][f]];

L[{x1, y1, z1}][f[x, y, z]]

(* ==> {0, 0, 0} *)

L[{x1, y1, z1}][f[x1, y1, z1]]


$$\left\{i h \text{z1} f^{(0,1,0)}(\text{x1},\text{y1},\text{z1})-i h \text{y1} f^{(0,0,1)}(\text{x1},\text{y1},\text{z1}),\\i h \text{x1} f^{(0,0,1)}(\text{x1},\text{y1},\text{z1})-i h \text{z1} f^{(1,0,0)}(\text{x1},\text{y1},\text{z1}),\\i h \text{y1} f^{(1,0,0)}(\text{x1},\text{y1},\text{z1})-i h \text{x1} f^{(0,1,0)}(\text{x1},\text{y1},\text{z1})\right\}$$

Here I defined p new, so the old definitions are erased at this point. The consequence is that L now gives zero angular momentum in the first case where the function has the wrong coordinates. Then you can later form products of functions with different coordinates and get the correct results for the total angular momentum, if you define a new function

Lt = Function[{f},L[{x1,y1,z1}][f]+L[{x,y,z}][f]];


for the total angular momentum.

To test the total angular momentum operator, let's define its components and its square, as above for the single angular momentum:

{Ltx, Lty, Ltz} = Map[Function[{f}, Lt[f][[#]]] &, {1, 2, 3}];

Lt2 = Function[{f}, Ltx[Ltx[f]] + Lty[Lty[f]] + Ltz[Ltz[f]]];

First[
Simplify[
Solve[
Lt2[(x + I y)^l (x1 + I y1)^l1] ==
eigenvalue (x + I y)^l (x1 + I y1)^l1, eigenvalue]
]
]

(* ==> {eigenvalue -> h^2 (l + l^2 + l1 + 2 l l1 + l1^2)} *)

Simplify[(eigenvalue /. %) == h^2 (l + l1) (l + l1 + 1)]

(* ==> True *)


Here I verified that the product of maximal wave functions for two particles with coordinates x, y, z and x1, y1, z1, respectively, and angular momenta l, and l1, is also a total-angular momentum eigenstate with maximal angular momentum quantum number l + l1.