Defining the Moyal Product in Mathematica

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[
    Head[#] === FourierTransform,
    # /. {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.