thanks for your quick answer Edgardo!
As I figured, this is so much nicer (and shorter) than how i was thinking of doing it. Also very helpful.
Small caveat is that the differential operator works slightly differently than the ave function that i was imagining because:
"ave(beta^2)" does not give "beta^2" (i.e last term), but it's trivial to just manually set it to 1 after the rest of the relevant ave(ABC) terms get transformed to the f-functions.
Also, thanks for sending the link to your worksheet... indeed that looks tangentially related.
This thing i'm doing, is essentially a mean-field expansion, which is then cut off at some nth order (i.e. meaning the nth order cumulant is approximated as zero). This is actually extremely widely used in various areas of physics, and comes up in many-many papers. One often wants to understand some dynamics of a large size system, for which solving the full Shrodinger (or more realistically Lindblad master) equation is not feasible.
Perhaps at some stage, when you're bored, you might consider adding this kind of functionality to the Physics package (probably an hour's work for you - the hardest parts, i.e handling, sorting, noncommuting operators - you've already done).
The idea is basically:
1) write down the eoms for a given system. For coherent system, that's just the standard Heisenberg eoms, and there is trivial correction if one has Lindblad (noise describing) terms.
2) average all the operators
Key point is that the eoms might not be closed (meaning a diff equation for, say, <AB> might depend on, say, <ABC>), hence:
3) do the n-th order cummulant approximation; i.e., write the all nth order averages in terms of terms only containing (n-1)-th order averages.
There is currently a way to do this up to 4th order in mathematica, using the Quantum package (not updated since 2011, but works with recent versions).
There is also a recent paper about a Julia package that does it (with some limitations)
In the past, I'be been doing this mainly by hand in some simpler cases, but more recently have switched to using the mathematica package above to generate relevant equations of motion, then slightly rewriting them so that maple can easily import via FromMma()... but indeed would be nice to do it directly in maple (hence this post). So will try to get some crude version going based on your answer.
I have another couple of simple questions about functions of operators and commutation, as well as about defining a visual representations for objects (i.e. can i have, say, f([A,B,C]) get printed as "<ABC>"), but maybe will ask those in their own thread at some stage.