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?

  • 2
    $\begingroup$ I’m voting to close this question because it should be asked on e.g. Mathematics SE. It is about general mathematical concept, not about quantum computing. $\endgroup$ Commented Dec 25, 2021 at 15:53
  • 4
    $\begingroup$ @MartinVesely I totally disagree. Mathematics.SE will not be the best place to ask about software on how to simplify expressions involving creation and annihilation operators. I'm voting to leave this question open, but if it does get closed, I'd recommend the asker to ask it on MMSE with the "software-recommendations" tag, or on Physics.SE. $\endgroup$ Commented Dec 25, 2021 at 21:35
  • $\begingroup$ Effectively you want some code that can manipulate noncommutative polynomials. For example, this can be done in sympy. $\endgroup$
    – Rammus
    Commented Dec 26, 2021 at 16:03

1 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,
        {l : nc[__], {first_, second_}},
        Sequence @@ {
            DirectedEdge[l, firstIfList@first],
            DirectedEdge[l, firstIfList@second]

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:

enter image description here

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

enter image description here

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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.