Defining the Moyal Product in Mathematica

What would be the best way to define the Moyal Product, $f(x,p)\star g(x,p)$, in Mathematica?

$f(x,p)\star g(x,p)$ may be written as $f\star g=\sum\limits_{n=0}^{\infty}\frac{1}{n!}\left(\frac{i\hbar}{2}\right)^n\Pi^n(f,g)$ where $\Pi^n=\sum\limits_{k=0}^n(-1)^k{n\choose k}\frac{\partial^nf}{\partial p^k\partial x^{n-k}}\frac{\partial^ng}{\partial x^k\partial p^{n-k}}$

I'm also wondering how I would then be able to define $f(x,p)\tilde{\star}g(x,p)$ as $f\star(1+\alpha\star_1+\alpha^2\star_2)g$ with $\star_1=cp^3\overleftarrow{\partial_p}\overrightarrow{\partial_p}$, $\star_2=c^2p^6\overleftarrow{\partial^2_p}\overrightarrow{\partial^2_p}$, and c being a constant.

Thanks.

$\Large{\rm{\bf{Edit}}}$

Here's how I tried to define the Moyal product:

Star[f_, g_] =
Sum[(I*h/2)^n*
Sum[(-1)^k*Binomial[n, k]*D[f[x, p], {p, k}, {x, n - k}]*
D[g[x, p], {x, k}, {p, n - k}], {k, 0, n}], {n, 0, Infinity}]


which output

$\sum _{n=0}^{\infty } 2^{-n} (i h)^n \sum _{k=0}^n (-1)^k \binom{n}{k} f^{(n-k,k)}(x,p) g^{(k,n-k)}(x,p)$

I know that $x\star p=xp+\frac{i\hbar}{2}$, so I tried Star[x,p], but got the following output:

$\sum _{n=0}^{\infty } 2^{-n} (i h)^n \sum _{k=0}^n (-1)^k \binom{n}{k} p^{(k,n-k)}(x,p) x^{(n-k,k)}(x,p)$

I think my problem is where I used f and g in the function of Star, but I'm not sure how to fix it. Once I figure out how to define $f\star g$ in Mathematica, I think I'll be able to define $f\tilde{\star} g$.

$\Large{\rm{\bf{Edit}}}$

I've included the $\frac{1}{n!}$ in the definition of the Moyal Product.

• It seems quite straightforward. What difficulties are you finding? Jul 4, 2014 at 16:38
• You may want to look at the Notation package. Jul 4, 2014 at 16:49
• In hindsight, I should have attached my first attempt at defining the function. I've edited my question. Jul 4, 2014 at 17:01
• Apparently the problem is related to a limit calculation. If you use star[f_, g_] := Sum[(I*h/2)^n* Sum[(-1)^k*Binomial[n, k]*D[f, {p, k}, {x, n - k}]* D[g, {x, k}, {p, n - k}], {k, 0, n}], {n, 0, 10}] it works as expected Jul 4, 2014 at 17:31
• There's also the fact that with your current definition you'd have to call Star[Function[{x, p}, x], Function[{x, p}, p]] instead of Star[x, p].
– user484
Jul 8, 2014 at 5:17

Edit

The definition of the Moyal product in the original question was missing a factor 1/n! under the sum. I used the textbook definition instead. A reference for this definition, without any unnecessary technicalities, is here: Quantum Mechanics in Phase Space (arXiv), see the appendix. Thanks to Rahul Narain for spotting the discrepancy between my textbook definition and the original form of the question.

In the following, I also set the constant h from the question equal to 1.

First part of the solution: polynomials

The best implementation depends on the function space that you're interested in. Based on the examples for f and g in your question, (f=x and g=p for example), I am going to assume first that the function space consists of the polynomials in the variables x and p. At the end, I'll address how to define the star product for more general functions.

Given two polynomials, we already know that for any two functions the series will reduce to a finite sum. We just have to make sure that we do the sum to a large enough number of terms. Instead of implementing the definition you give in the question, however, I decided to use an equivalent but different definition, in which the operators can be ordered so that they are to the left of both functions.

The way to do this is simply to write the "dyadic" derivatives appearing in the original question as follows:

$$f(x, p)\frac{\overleftarrow{\partial}}{\partial x}\frac{\overrightarrow{\partial}}{\partial p} g(x,p)= \left[\frac{\partial^2}{\partial x_1\partial p} f(x_1, p_1) g(x,p) \right]_{x_1\to x, p_1\to p}$$

Then the series expansion can be written in terms of the (normal) operator exponential. In the answer to the lined question, I defined this operator using Fold, but this required knowing the expansion order beforehand.

In our case, the assumption of polynomial functions f and g alone doesn't tell us the maximum number of terms that will be needed. Instead, we should define the operator exponential such that its recursive formula automatically terminates when the sum no longer changes (this cutoff is guaranteed to exist for polynomials).

This is why I implement another version of the operatorExp function below, using FixedPoint to sum up all terms in the usual exponential series until the result no longer changes. The construction is a little uglier than using the linked approach based on Fold because I have to keep track of the powers "manually" in a list (containing three elements, the first of which is the partial sum) that is passed on in each recursion step.

Using the first equation, the Moyal product can be written as

$$f\star g=\lim_{x'\to x}\lim_{p'\to p}\sum_{\nu=0}^{\infty}\frac{1}{\nu!}\left(\frac{i}{2}\right)^{\nu}\left(\frac{\partial}{\partial x'}\frac{\partial}{\partial p}-\frac{\partial}{\partial p'}\frac{\partial}{\partial x}\right)^{\nu}f(x',p')\, g(x,p)$$

which is an exponential:

$$f\star g=\lim_{x'\to x}\lim_{p'\to p}\exp\!\left[\frac{i}{2}\left(\frac{\partial}{\partial x'}\frac{\partial}{\partial p}-\frac{\partial}{\partial p'}\frac{\partial}{\partial x}\right)\right]\,f(x',p')\, g(x,p)$$

After defining the exponential in such a way that I don't have to evaluate factorials and binomials at every step, I define the exponent (poissonOp, in the square brackets above) as an operator acting on the four variables $x_1, p_1, x, p$ introduced above.

Finally, I define the star product, but only for polynomials.

Clear[operatorExp];
operatorExp[dop_, n_: 100][f_] := First@FixedPoint[
{#[[1]] + dop[#[[2]]], dop[#[[2]]]/#[[3]], #[[3]] + 1} &, {f, f, 2},
n, SameTest -> (PossibleZeroQ[#[[1]] - #2[[1]]] &)]

Clear[poissonOp];
poissonOp[x1_, p1_, x_, p_] :=
Function[f, I/2 (D[f, x1, p] - D[f, p1, x])];

Clear[star];
star[f_?(PolynomialQ[#, {x, p}] &), g_?(PolynomialQ[#, {x, p}] &)] :=
Module[{x1, p1},
operatorExp[
poissonOp[x1, p1, x, p]][(f /. {x :> x1, p :> p1}) g] /. {x1 :>
x, p1 :> p}
]


Here are some tests:

star[x, p]

(* ==> I/2 + p x *)

star[x, p] - star[p, x]

(* ==> I *)

star[x, p^2]

(* ==> I p + p^2 x *)

star[x^6 p^2, p^7 x^4]

(*
==> (945 p x^2)/4 - 2205/2 I p^2 x^3 - (2205 p^3 x^4)/4 -
2205/2 I p^4 x^5 - 840 p^5 x^6 + 42 I p^6 x^7 - (153 p^7 x^8)/2 +
17 I p^8 x^9 + p^9 x^10
*)


In principle, this can be applied to arbitrary polynomials. In practice, I put in a "failsafe" to prevent hanging, by adding a third argument to FixedPoint in operatorExp, specifying the maximum number of terms to which the recursion is pushed. The default value is n=100, but it could in principle be made arbitrarily large. Again, this assumes that you know your function space to be polynomials.

If, on the other hand, the functions are not polynomials, then you could still use the Fourier transform approach I also mention in the answer linked above. This will also be made easier by rewriting the left and right acting derivatives as I did here.

Edit: More general functions

Elaborating on the Fourier-transform approach, you should also be able to use the following definition (limited only by Mathematica's ability to do the required transforms):

star[f_, g_] := FourierTransform[
InverseFourierTransform[f, {x, p}, {kf, lf}]
InverseFourierTransform[g, {x, p}, {kg, lg}]
Exp[I/2 (-kf lg + kg lf)],
{kf, lf, kg, lg},
{x1, p1, x, p}
] /. {x1 -> x, p1 -> p}


This uses the fact that the derivative operators turn into numbers when applied to the Fourier basis functions. Here is an example using Gaussians of different widths:

star[Exp[-x^2 - p^2], Exp[-x^2/2 - p^2/2]]

(* ==> 2/3 E^(1/3 (-3 p^2 + (p - I x)^2 - 3 x^2 + (I p + x)^2)) *)


The Fourier method also works for polynomials, but is slower in that case because it requires some pretty complicated manipulations with Dirac delta functions (which Mathematica does correctly, though). For speed reasons, both definitions of star can be used side by side. The correct choice of method will be made based on f and g.

Edit: brute-force fix for pathological functions

Here is a version of the Fourier approach that attempts to deal with cases where the required multi-dimensional Fourier transform fails. Since this was motivated by Karan's comment regarding a pathological function with non-compact support, I don't claim that it's the unique and best way to deal with such cases. But one can argue that re-arranging the order of integration variables is a good first try when an integration fails:

star[f_, g_] := If[
# /. {kg, lg, kf, lf} -> {kf, lf, kg, lg}, #] &[
FourierTransform[
InverseFourierTransform[
f, {x, p}, {kf, lf}] InverseFourierTransform[
g, {x, p}, {kg, lg}] Exp[I/2 (-kf lg + kg lf)], {kg, lg, kf,
lf}, {x1, p1, x, p}]] /. {x1 -> x, p1 -> p}


An example where this is needed would be a pathological function such as f = 1/x^2. In Mathematica 8.0.4, there is no problem doing the required Fourier transforms in star[1/x^2,p] or star[p,1/x^2], independently of the order, provided you first set \$Assumptions=x>0. But in versions 9 and 10, you will need the last definition of star. What happens is that the FourierTransform initially remains unevaluated for star[p,1/x^2] (Mathemtatica tries for a while but then gives up). In that case, the If statement takes the unevaluated result and switches the order of the Fourier integration variables. The justification for this is that it produces a result which differs from the opposite order star[1/x^2,p] in the expected way, so it preserves the desired skew symmetry of the Moyal bracket.

• Using the Fourier transform code, I tried to compute star[1/x,p] with Assumptions x>0 and it gives the correct answer but with star[1/x^2,p] I get FourierTransform[-2 I E^(1/2 I (kg lf - kf lg)) kf \[Pi]^2 DiracDelta[kg] DiracDelta[lf] Sign[kf] Derivative[1][DiracDelta][lg], {kf, lf, kg, lg}, {x, p, x, p}] whereas star[p,1/x^2] gives the correct answer. How can this be fixed? Jul 7, 2014 at 22:18
• @Karan Indeed, I checked and your observation is a bug in Mathematica version 9.0.0 which has been fixed in version 9.0.1.
– Jens
Jul 8, 2014 at 4:44
• Sorry, by "final test" I meant the final one of your "Here are some tests" block, namely star[x^6 p^2, p^7 x^4]. I also should have mentioned that I took the sum up to 100 terms and the result didn't change, so I assumed that it had converged. As you stated, for any two polynomials the series will reduce to a finite sum, and that is what the above result seems to be.
– user484
Jul 8, 2014 at 6:12
• @RahulNarain Oh, I didn't even notice that the definition of the star product in the question is wrong. The textbook definition is the exponential I wrote down. If you believe my version (which you can find in the literature), then the definition in the question is missing a factor 1/n! under the sum. Also, I of course set h=1. So the discrepancy won't show up unless higher terms n>1 are involved in the sum. I can add some references, but I haven't heard from the OP at all so I don't know if it's worth it. Thanks for checking!
– Jens
Jul 8, 2014 at 7:03
• @Karan OK - ConditionalExpressions are of course also a warning for you to stop and think about the validity of what you're doing... basically the order of integrations seems to matter, and boundary terms could be cropping up in partial integrations. I think that goes beyond what the question asked for. But I could post a brute-force fix for this issue anyway.
– Jens
Jul 8, 2014 at 18:18
star[f_, g_] :=
Sum[(I \[HBar]/2)^n 1/
n! Nest[ D[#, p, x1] - D[#, x, p1] &, (f /. {x -> x1, p -> p1})*g,
n], {n, 0, Plus @@ Exponent[f, {x, p}]}] /. {x1 -> x, p1 -> p}

AbsoluteTiming @ star[x^6 p^2, p^7 x^4]

{0.00201501,
p^9 x^10 + 17 I p^8 x^9 \[HBar] - 153/2 p^7 x^8 \[HBar]^2 +
42 I p^6 x^7 \[HBar]^3 - 840 p^5 x^6 \[HBar]^4 -
2205/2 I p^4 x^5 \[HBar]^5 - 2205/4 p^3 x^4 \[HBar]^6 -
2205/2 I p^2 x^3 \[HBar]^7 + 945/4 p x^2 \[HBar]^8}

• Hi F. Kampas, welcome to SE! Do you think you could include an example of your code? Thanks! Apr 21, 2018 at 21:33
• AbsoluteTiming @ star[x^6 p^2, p^7 x^4] {0.00201501, p^9 x^10 + 17 I p^8 x^9 [HBar] - 153/2 p^7 x^8 [HBar]^2 + 42 I p^6 x^7 [HBar]^3 - 840 p^5 x^6 [HBar]^4 - 2205/2 I p^4 x^5 [HBar]^5 - 2205/4 p^3 x^4 [HBar]^6 - 2205/2 I p^2 x^3 [HBar]^7 + 945/4 p x^2 [HBar]^8} Apr 23, 2018 at 0:10