Something that comes recurrently when you work with quantum fields is, given any operator that consists of products of powers of annihilation \(a\) and creation \(a^\dagger\) Bose operators, such as, e.g.,
\[a^{\dagger3}a^2a^{\dagger3}a^2a^{\dagger}a\]
compute its normal order[1]:
\[:a^{\dagger3}a^2a^{\dagger3}a^2a^{\dagger}a:\] You can do this tediously by iterating the fundamental relation:
\[aa^\dagger=1+a^\dagger a\,.\]
It's some algebra thing. Soon you get to know rules like:
\[a^na^\dagger=na^{n-1}+a^\dagger a^n\,,\]
which, written as a commutator, shows that \(a^\dagger\) does some sort of derivative (a similar relation holds for \(a^{\dagger m}\)):
\[[a^n,a^\dagger]=na^{n-1}\,.\]
And there is something deep here, but it's not my point. My point is, when you get to compute something like the above, how do you do, technically? I've always taken as granted that it's an obvious problem for a computer, just iterating. It doesn't appear so trivial to put in the computer, though, that it would take more time to program rather than compute it yourself (stupidly, by iterations or using rules as above or even more powerful one such as the one for \(a^ma^{\dagger^n}\) which is a sum over normally ordered terms of reducing orders with some products of binomials).
Well, I got tired of iterating this lazy reasoning and put the thing in the computer at last. It's not very complicated, at least if you have a grisp of λ-calculus, such as lisp or Mathematica. It took me a day, the larger part of it because I got stuck on how to set apart the concept of list in Mathematica with the list as a defining element of my objects. This was solved by deferring application of my functions so that I could flatten and deal with normal lists without hampering on the "objects" themselves.
My fundamental "object" is a list of two elements:
\[\{\{\ldots, i_3,i_2,i_1\},\alpha\}\,.\]
Here we're going to play with "objects" like they do in set theory when they define numbers, where, e.g.,
3 = {0,1,2} = {0, {0}, {0, {0}}} = {{ }, {{ }}, {{ }, {{ }}}}
with 0 the emptyset \(\emptyset=\{\}\), or, the way surreal numbers are defined.
In my list, the first element is another list, of integers, and the second is a coefficient. The list of integers contains the powers of operators in the sequence:
\[\ldots a^{\dagger}aa^{\dagger}aa^{\dagger}aa^{\dagger}aa^{\dagger}a\,.\]
that is formed by products of \(a^{\dagger}a\), so the list can be of any size. The list above defines: (the coefficient comes as a coefficient!):
\[\alpha\ldots a^{i_3}a^{\dagger i_2}a^{i_1}\,.\]
So the operator we started with above is, for the computer:
\[a^{\dagger3}a^2a^{\dagger3}a^2a^{\dagger}a\rightarrow\{\{3,2,3,2,1,1\},1\}\,.\]
A list of these objects is a sum of them, so:
\[1+a^\dagger a\rightarrow\{\{\{0,0\},1\},\{\{1,1\},1\}\}\,.\]
The task for the programmer is to set up a code that, with the above definition, turns any list into a normally ordered one. A normally ordered list is simple enough, it's a list of lists whose first element is at most a doublet.
This module[2], let us call it NormalOrder, would act as follows: (first line is the input, second line is the result)
NormalOrder[{{{1, 1, 0}, 1}}]
{{{0, 0}, 1}, {{1, 1}, 1}}
This is \(aa^\dagger=1+a^\dagger a\) (you will recognize the output as the previous expression and the lhs as our fundamental relation).
Another well known result:
NormalOrder[{{{7, 1, 0}, 1}}]
{{{0, 6}, 7}, {{1, 7}, 1}}
which is another particular case we've already seen above. More difficult, although still known through a closed form expression (that I didn't give explicitly):
NormalOrder[{{{7, 7, 0}, 1}}]
{{{0, 0}, 5040}, {{1, 1}, 35280}, {{2, 2}, 52920}, {{3, 3},
29400}, {{4, 4}, 7350}, {{5, 5}, 882}, {{6, 6}, 49}, {{7, 7}, 1}}
Or the problem I started with:
NormalOrder[{{{3, 2, 3, 2, 1, 1}, 1}}]
{{{4, 2}, 12}, {{5, 3}, 24}, {{6, 4}, 10}, {{7, 5}, 1}}
which, if I didn't blunder, brings us to\[:a^{\dagger3}a^2a^{\dagger3}a^2a^{\dagger}a:=12a^{\dagger4}a^2+24a^{\dagger5}a^3+10a^{\dagger6}a^4+a^{\dagger7}a^5\,,\]
which doesn't look obviously incorrect (same signature of exponents is retained). I didn't check carefully the code yet. In fact I won't even give it for now, only because it's a bit late and that I will leave it as an exercise for the reader, for the time being.
I am not sure, yet, until I study it further, how much is accomplished here. Not much more than iterating, but universality comes from this, so it's hardly a point. In principle, we're not very far off here of starting from any quantum Hamiltonian or Liouvillian and derive completely automatically all the equations of motions for any observable, to any desired order. It's just to do the same for anticommutation rules and build more complex objects that bring us from the vector space structure of here to that of tensor product (I've already got the symbolic code that sets up and solve the equation of motions, but it needs someone went through the algebra). Some difficulties would arise for fermions but for commuting operators, like qubits, this is maybe merely another day's work. I'll leave it there for tonight.