Is there a software to simplify commutators such as $[{\hat a}^\dagger{\hat a}^\dagger {\hat a},{\hat a}^\dagger\hat a \hat a]$?

Question

Sometimes I need to simplify some commutators in the form $[{\hat a}^\dagger{\hat a}^\dagger {\hat a},{\hat a}^\dagger\hat a \hat a]$.

Doing it by hand is exhausting and boring. Is there software (Matlab, Python, Mathematica, etc) that can handle such commutators?

Answer 1

I actually wrote a short script a bunch of years ago to do something similar in Mathematica. The code was written to visualise how the different integer coefficients you get in the final expressions come from the commutator algebra. Also, it does not make any assumption on what the actual commutator look like, meaning it was written to take expressions with products of two non-commuting symbols and rewrite them in terms of nested commutator etc. However, a slight modification of that code should work for the simpler case of $[a,a^\dagger]=1$.

The code defining the necessary functions is:

<< MaTeX`

distribute[args_] := (args //. {
    HoldPattern[nc[l___, Plus[m__], r___]] :> Total[nc[l, #, r] & /@ {m}],
    nc[l___, c_*nc[m__], r___] :> c nc[l, m, r],
    nc[l___, nc[m__], r___] :> nc[l, m, r],
    nc[-a_, b_] :> -nc[a, b],
    nc[a_, -b_] :> -nc[a, b],
    nc[nc[l : __], r_] :> nc[l, r], nc[l_, nc[r : __]] :> nc[l, r],
    nc[a_] :> a,
    nc[a__, 1] :> nc[a],
    nc[1, a__] :> nc[a]
});

comm[a_, b_] := nc[a, b] - nc[b, a];

groupPowers[args_] := args //. {
    nc[l___, a_, a_, r___] :> nc[l, a^2, r],
    nc[l___, Power[a_, n_], a_, r___] :> nc[l, a^(n + 1), r],
    nc[l___, a_, Power[a_, n_], r___] :> nc[l, a^(n + 1), r]
};

comm[a_, b_] := nc[a, b] - nc[b, a];

makeCsToBrackets[expr_] := expr //. {
    c[a_, b_] :> "[" <> ToString @ makeCsToBrackets @ a <> "," <> ToString @ makeCsToBrackets @ b <> "]",
    Power[a_, n_] :> ToString @ makeCsToBrackets @ a <> "^" <> ToString @ n
};

beautify[expr_] := makeCsToBrackets@expr //. {
    nc[args__] :> MaTeX[StringJoin @@ (ToString /@ {args}), Magnification -> 1.5],
    s_String :> MaTeX[s, Magnification -> 1.5]
};

singleStepExpand[expr_] := expr /. {
     nc[l___, a, b, r___] :> {nc[l, b, a, r], nc[l, c[a, b], r]},
     nc[l___, a, m : c[__], r___] :> {nc[l, m, a, r], nc[l, c[a, m], r]},
     nc[l___, m : c[__], b, r___] :> {nc[l, b, m, r], nc[l, c[m, b], r]}
 } // Map@distribute;

stepExpand[expr_] := singleStepExpand@expr // If[Length@# > 1, stepExpand /@ #, #] &;

stepExpandFullStory[expr_] := singleStepExpand@expr // (
    If[Head @ # === List,
        Append[{expr}, stepExpandFullStory /@ #],
        #
    ] &
);

firstIfList[expr_] := If[Head @ expr === List, First @ expr, expr];

nestedListToListOfEdges[expr_] := Cases[expr,
    RuleDelayed[
        {l : nc[__], {first_, second_}},
        Sequence @@ {
            DirectedEdge[l, firstIfList@first],
            DirectedEdge[l, firstIfList@second]
        }
    ],
    All
];

edgesToGraphWithNiceLabels[edges_] := Graph[edges,
    VertexLabels -> Map[
        # -> beautify @ groupPowers @ # &,
        DeleteDuplicates @ Flatten[edges, Infinity, DirectedEdge]
    ],
    GraphLayout -> "LayeredDigraphEmbedding"
];

Note that the only external library I'm using is MaTeX, which is used to pretty-print the output results in LaTeX. If you don't care about that, you don't need it (also, in most cases it will take more time to generate the LaTeX outputs than to compute the actual expressions, for if you care about performance, definitely remove the pretty-printing parts). If you want to use this repeatedly in a larger piece of software, I'd also recommend putting this code into a package to properly scope the functions.

An usage example would be:

Note that here nc stands for non-commutative product, and a and b are generic operators/elements of the algebra.

To instead assume the commutator relation $[a,a^\dagger]=1$, and obtain as a final expression a sum of terms, we can slightly tune the above code defining the following functions:

singleStepExpand2[expr_] := expr /. {
     nc[l___, a, b, r___] :> Plus[nc[l, b, a, r], nc[l, r]]
 } // Map@distribute;
stepExpand2[expr_] := singleStepExpand2 @ expr // If[Length @ # > 1, stepExpand2 /@ #, #] &;

Again, this is using a notation with $[a,b]=1$, that is, $b\equiv a^\dagger$. For example, to handle the case in the question, we can use

comm[nc[b, b, a], nc[b, a, a]] // stepExpand2 // groupPowers // beautify

I should also point out that, because this was written with in mind a generic algebra of two operators, with no assumptions on the relations between the different nested commutators, that are certainly more efficient ways of doing this for the specific case at hand. Also, I haven't really tested this thoroughly, so let me know if you find examples in which it produces the wrong result.

There's also a bunch posts on mathematica.SE related to implementing noncommutative algebras. You might want to have a look at e.g. https://mathematica.stackexchange.com/q/95616/27539, and links therein.