# Why is the square of a differential operator not evaluated automatically?

I defined an operator:

py[f_[x, y]] := (I*h*D[f[x, y], y] - x) // Expand


(not important, but it represents a simplified component of a quantum momentum in a magnetic field). This works:

py[f[x, y]]


yielding

This does not:

py[py[f[x, y]]]


Why?

• The pattern matches only if the argument of py is a function with arguments that are exactly x and y. – QuantumDot Mar 24 '16 at 11:32
• @QuantumDot So, how should I extend the definition, so that my case would also be accounted for? – Alexei Boulbitch Mar 24 '16 at 12:44
• I would do something like this: py = (I h D[#, y] - x &) – chuy Mar 24 '16 at 14:00
• I need a little more info: are py[x] and py[y] (without f in them) supposed to vanish? If so, the construction by @chuy won't work – QuantumDot Mar 24 '16 at 14:01
• true, what result do you want for py[py[f[x,y]]? – chuy Mar 24 '16 at 14:02

## 2 Answers

Your definition of the momentum operator is incorrect because as a linear operator it needs to have the function f appearing linearly in all terms. This includes the x in your definition.

To fix this, I would define the operator as follows (using ψ instead of f):

py = Function[ψ, Expand[-I D[ψ, y] + ay[x, y] ψ]];

ay[x_, y_] := x

py@py@ψ[x, y]


$$x^2 \psi (x,y)-2 i x \,\psi ^{(0,1)}(x,y)-\psi ^{(0,2)}(x,y)$$

Here, I added a minus sign to your derivative to make it consistent with quantum mechanics, and left the definition of the vector potential unspecified in py. I define it only after defining py. Note that ay[x, y] ψ appears, not just ay[x, y].

Also, since you have implicitly decided to use x and y as the global differentiation variables, it's simply not necessary to explicitly include them as a pattern to be matched by the operator. This is why I define the operator as a Function where the wave function appears only as a single variable ψ.

As you can see, the square of the operator is now evaluated correctly using py@py.

Update

Based on Jens' answer, I am adapting the code. The idea is the same, but it calculates the derivative correctly:

ClearAll@py
SetAttributes[py, HoldFirst]
py[f_[x_, y_], order_] := Expand@Nest[I h D[#, y] - x # &, f[x, y], order]


Usage:

ClearAll@f
py[f[x, y], 2]
(* x^2*f[x, y] - 2*I*h*x*Derivative[0, 1][f][x, y] - h^2*Derivative[0, 2][f][x, y] *)

f[x_, y_] := x^3 y^3
py[f[x, y], 2]
(* -6 h^2 x^3 y - 6 I h x^4 y^2 + x^5 y^3 *)


Original Post

Here. Instead of using f_[x, y] as the pattern for the argument of py, just use a general expression and take the derivative with respect to it, while putting in the variables explicitly as arguments.

ClearAll@py
py[expr_, {x_, y_}] := (I*h*D[expr, y] - x) // Expand
py[expr_, {x_, {y_, order_}}] := Nest[py[#, {x, y}] &, expr, order]


Usage:

py[f[x, y], {x, y}]
py[py[f[x, y], {x, y}], {x, y}]
(* -x + I h Derivative[0, 1][f][x, y] *)
(* -x - h^2 Derivative[0, 2][f][x, y] *)


or

py[f[x, y], {x, {y, 1}}]
py[f[x, y], {x, {y, 2}}]
py[%, {x, {y, 1}}]
(* -x + I h Derivative[0, 1][f][x, y] *)
(* -x - h^2 Derivative[0, 2][f][x, y] *)
(* -x - I h^3 Derivative[0, 3][f][x, y] *)


Alternatively, if you don't want to have to put the variables in explicitly, you could do this. This does not allow you to call py on the expression after you've evaluated it, but that's why we build in a certain number of derivatives to begin with.

ClearAll@py
SetAttributes[py, HoldFirst]
py[f_[x_, y_], order_] := Nest[I h D[#, y] - x &, f[x, y], order]


Usage:

f[x_, y_] := x^3 y^3
py[f[x, y], 3]
(* -x - 6 I h^3 x^3 *)

ClearAll@f
py[f[x, y], 3]
(* -x - I h^3 Derivative[0, 3][f][x, y] *)