# Exponential of a Differential Operator

In Mathematica, is it possible to exponentiate a differential operator such that the operator will act on a function, $f(x,p)$? Specifically, I wondering if I can get Mathematica to do this:

$\exp(c\frac{d}{dp})f(x,p)=(1+c\frac{d}{dp}+\frac{c^2}{2}\frac{d^2}{dp^2}+\cdots)f(x,p)=(f(x,p)+c\frac{df}{dp}+\frac{c^2}{2}\frac{d^2f}{dp^2}+\cdots)$

where $c\in\Re$.

Thanks.

• How do you expect mathematica to treat the ... bit? One approach would be Series[f[x, p + c], {c, 0, 2}] but that's because I know what the resummed series looks like.
– acl
Jun 25, 2014 at 15:37
• I'd like Mathematica to treat the ... as the continuation of the exponential function, namely $\frac{c^n}{n!}\frac{d^nf}{dp^n}$ for the nth term. Jun 25, 2014 at 17:24
• Yes, I appreciate that. What I'm asking is what would you expect mathematica to do where you wrote ...?
– acl
Jun 25, 2014 at 20:32

There are two possible approaches:

(1)

To map an operator exponential onto a MatrixExp, one would have to assume that there exists an orthogonal basis of the function space in which the target functions f[x] can be expanded. Then construct a matrix representing the operator in that basis. For example, the momentum operator $p$ is diagonal in the Fourier basis, and therefore we can do the exponential very easily in that basis. But we may not always be able to determine the expansion coefficients in a chosen basis in closed form. As a result, you end up with matrices corresponding to a truncated expansion of the function f[x]. Exact results can be obtained in this approach if your function space is finite-dimensional and the operator doesn't throw you out of that space.

(2)

Instead of truncating the expansion of f[x] in some basis, you can truncate the series for the exponential at some order of the expansion. Truncation in either approach genrally can't be avoided, but this approach is very common since the exponential series is absolutely convergent over the complex numbers.

However, the usual form of the series expansion for the exponential requires you to calculate factorials at every order. This is inefficient from a computational point of view.

Instead, the exponential of an operator can be conveniently (and efficiently) defined by using a recursive formula for the exponential. This is implemented here:

Clear[f, operatorExp];
operatorExp[dop_, n_] := Function[f, Fold[f + dop[#1]/#2 &, f, Reverse@Range[n]]]


This one-liner defines the exponential of dop up to a number n of terms. The result is itself an operator, taking as its argument a function of the form f[x]. Here I'll illustrate what this does to a general operator dop:

Clear[dop,x];
operatorExp[dop, 4][f[x]]


$\text{dop}\left(\frac{1}{2} \text{dop}\left(\frac{1}{3} \text{dop}\left(\frac{1}{4} \text{dop}(f(x))+f(x)\right)+f(x)\right)+f(x)\right)+f(x)$

Now try it for a well-known example (the Taylor series), and then for an operator with non-commuting terms. To do this, I use Function to define the operators dOp and dOp2. They assume that the independent variable of the differential operation is x.

dOp = Function[f, D[f, x]];
operatorExp[dOp, 3][f[x]]


$f(x) + f'(x)+\frac{1}{2} \left(\frac{1}{3} f^{(3)}(x)+f''(x)\right)$

dOp2 = Function[f, x + D[f, x]];
operatorExp[dOp2, 3][f[x]]


$x + f(x) + f'(x)+\frac{1}{2} \left(\frac{1}{3} f^{(3)}(x)+f''(x)+1\right)$

Applied to the specific problem in the question

Here I use the momentum as the differentiation variable in the operator:

Clear[c,p];
dOp3 = Function[f, c D[f, p]];
Expand[operatorExp[dOp3, 3][f[x, p]]]


$\frac{1}{6} c^3 f^{(0,3)}(x,p)+\frac{1}{2} c^2 f^{(0,2)}(x,p)+c f^{(0,1)}(x,p)+f(x,p)$

By using Expand, I turned the recursively factored form back into something that looks like the usual Taylor expansion.

• Perfect! Your second approach is exactly what I needed. Jun 27, 2014 at 17:20
• Is the result for $operatorExp[dOp2, 3][f[x]]$ correct? $x$ and $D[f, x]$ are non-commuting operators. Aug 19, 2017 at 13:16
• @wondering Yes, it's correct. It doesn't matter whether or not the terms in dop commute, as long as dop is a linear operator.
– Jens
Aug 19, 2017 at 15:26
• But if I calculate the first three terms of the exponential expansion, i.e. $\Exp[x+\partial] f=(1+(x+\partial_x)+(x+\partial_x)^2+(x+\partial_x))^3+\dots)f(x)$, then I don't get the result above. How comes that you get terms like $x$ or $1/2$ in the series, if $f(x)$ has to appear in every term as a derivative of some order ($0,1,2,...$)? Aug 23, 2017 at 14:35
• @wondering You're simply not doing it right,
– Jens
Aug 23, 2017 at 16:08

may be something like this:

Plus @@ (c^#/#! D[f[x, p], {p, #}] & /@ Range[0, 5])


or

Sum[c^i/i! D[f[x, p], {p, i}], {i, 0, 5}]


Here is a trick for making a series expansion of a function of a single operator — i.e. the series expansion has no non-commutation issues — and then correctly manipulating powers of the operator, where the operator may have non-commuting internal parts.

Expand the exponential Exp[op] of an operator op — 3 terms are sufficient for illustrative purposes.

Series[Exp[op], {op, 0, 3}] // Normal

(* 1 + op + op^2/2 + op^3/6 *)


Convert the leading 1 into an identity operator id.

% // ReplacePart[#, 1 -> id] &

(* id + op + op^2/2 + op^3/6 *)


Operate on a state s.

(%)[s] // Through

(* id[s] + op[s] + ((op^2)/2)[s] + ((op^3)/6)[s] *)


Convert the powers of operators into the equivalent composition of operators. This trick sets things up so that, when op is expanded out in terms of its internal parts, any non-commuting parts in op are handled correctly.

% /. ((u_: 1) op^n_)[v_] :> u (Composition @@ Table[op, {n}])[v]

(* id[s] + op[s] + 1/2 op[op[s]] + 1/6 op[op[op[s]]] *)


Now use coordinate representation s[x] (i.e. a wave function) for the state s, and (for example) use (# + I D[#, x]) &) for the operator op .

% /. s -> s[x] /. {id -> (# &), op -> ((# + I D[#, x]) &)} // Expand

(* (8*s[x])/3 + (5/2)*I*Derivative[s][x] - Derivative[s][x] - (1/6)*I*Derivative[s][x] *)


This approach is very general — i.e. the (# + D[#, x])& is merely an example of the sort of operator you can use.

One way would be to get the series by treating the operator as a variable, then replace the variable with operator. Let say, you want the series to order 9.

n = 9;
op = Series[Exp[c dp], {dp, 0, n}] // Normal;
op /. {dp^m_ -> D[#, {p, m}]} &@f[x, p]


Here I first treat the operator as variable dp first and then replace it with the operator once I have my series. And you have your answer. You can use it for any function of operator (I mean in general, not those tricky function with singularity or even worse complexity!)

For any arbitrary form of operator

I can still use the same method, but the power will be taken care by Nest which will allow the repeated application of the operator.

opd[x_]=(D[#,x]#+x #)&; (*The operator*)
func[op_]:=Exp[c op]; (*The func. of operator*)
f[x_,p_]=y[x ,p]; (*Target function*)

n=1; (*order of series expansion*)
opr[op_]=op Series[func[op],{op,0,n}]//Normal;
out[g_,x_,p_]=(opr[op])/.op^m_.->Nest[opd[x],g[x,p],m-1]//Quiet;

out[f,x,p]//Expand


The output is $y[x,p]+c x y[x,p]+c y[x,p] y^{(1,0)}[x,p]$.

Now you can use any form or function of operator.

Regarding this I like to mention one thing. While making the series opr[op] I multiply it with op so that the power series starts from op instead of 1. That is why I use the Nest m-1 times for op^m. Because 1 is not interpreted as op^0 and in final result it shows 1 instead of y[x,p]. I don't know if there is any better way to handle it. I post this question little ago. I would be glad to hear comments about that which may also improve this answer.

Better Approach

opd[x_]=(D[#,x]#+x #)&;(*The operator*)
func[op_]:=Exp[c op];(*The func.of operator*)
f[x_,p_]=y[x,p];(*Target function*)
n=1;(*order of series expansion*)
opr[op_]= Series[func[op],op,0,n}]//Normal;
out[g_,x_,p_]=Sum[Coefficient[opr[op],op,m]Nest[opd[x],g[x,p],m],{m,0,n}];

out[f,x,p]//Expand

• What if the operator has parts that don't commute with each other? eg $\partial_x + x$? Does this work correctly then?
– acl
Jun 25, 2014 at 17:14
• If, instead of c dp, would it still work with $c*p^3*(dp)^2+g*p^2*dp$ where $g$ is also real and $(dp)^2$ corresponds to $\frac{d^2}{dp^2}$? Jun 25, 2014 at 17:21
• @acl and user85503 : I am afraid this method will work only for commutating elements. For noncommutating operators the series expansion may not be valid. However in that case you can define the series expansion in terms of commutators and then do the replacement. Or you can replace the (operator)^m using Jens' second approach above. Jun 25, 2014 at 22:51
• yes, it's a bit more work if the operator has noncommuting parts
– acl
Jun 26, 2014 at 0:07
• Thanks Sumit. Your answer combined with Jens regarding non-commuting operators is precisely what I wanted to do. Jun 27, 2014 at 17:19