BosonNormalOrder is a Mathematica code I wrote to normal-order bosonic operators, e.g., like the following:
$$aa^{\dagger k}a^la^\dagger=a^{\dagger(k+1)}a^{l+1}+(k+l+1)a^{\dagger k}a^l+kla^{\dagger(k-1)}a^{l-1}\,.$$
The implementation relies on what I call a "weighted sequence" (wseq), which is of the type
{w, {k, l, m, n}}
where w is the weight and k, l, m, n, etc. (the list can be of any size) are the exponents of a product of annihilation and creation operators, always starting with an annihilation one on the right (normal order), so, e.g.,
{1, {1, 1}}
corresponds to $\ud{a}a$ while
{2, {1, 1, 0}}
is $2a\ud{a}$. One can make more complex sequences in this way, e.g.,
{-1, {1, 2, 3, 1}}
is $-\ud{a}a^2\ud{a}^3a$. Various such weighted sequences in a list are understood as a sum of them, e.g.
{{-1, {1, 2, 3, 1}}, {2, {1, 1, 0}}}
is $-\ud{a}a^2\ud{a}^3a+2a\ud{a}$. Etc. What the module does is, given an expression in any order, it returns it normally-ordered version, e.g.:
BosonNormalOrder[{{1, {1, 1, 0}}}]
{{1, {0}}, {1, {1, 1}}}
For simplicity, if there is only one correlator only, instead of passing the full wseq, one can use:
BosonNormalOrderCorrelator[{1, 1, 0}]
{{1, {0}}, {1, {1, 1}}}
This is $a^2\ud{a}^2=2+4\ud{a}a+\ud{a}^2a^2$:
BosonNormalOrderCorrelator[{2, 2, 0}]
{{2, {0}}, {4, {1, 1}}, {1, {2, 2}}}
More details on the code and notebook can be foudn in this blog post.