Sympy: How to get simplified commutators using the second quantization module?

Question

So I want to use the fact that [b,bd]=1, where [] is the commutator, to get some commutators of more complicated expressions using sympy instead of doing it by hand, but instead, I get huge expressions that contain the commutator but it is not replaced for 1, here's the code

from sympy import *
from sympy.physics.secondquant import * 
comm1=simplify(Commutator(B(0),Bd(0)).doit())
comm1

the output in this case is 1 , it corresponds to the [b,bd]=1 case, but if I input a more complicated expression such as

w1,w2,g=symbols('w1 w1 g')
H=w1*B(0)*Bd(0)+w2*B(1)*Bd(1)+g*Bd(0)*B(1)+conjugate(g)*Bd(1)*B(0)

comm2=simplify(Commutator(H,B(0)))
print(simplify(comm2))

I get

-g*(AnnihilateBoson(0)*CreateBoson(0)*AnnihilateBoson(1) - CreateBoson(0)*AnnihilateBoson(1)*AnnihilateBoson(0)) - w1*(AnnihilateBoson(0)*AnnihilateBoson(1)*CreateBoson(1) - AnnihilateBoson(1)*CreateBoson(1)*AnnihilateBoson(0)) + w1*(AnnihilateBoson(0)*CreateBoson(0)*AnnihilateBoson(0) - AnnihilateBoson(0)**2*CreateBoson(0)) - conjugate(g)*(AnnihilateBoson(0)*CreateBoson(1)*AnnihilateBoson(0) - CreateBoson(1)*AnnihilateBoson(0)**2)

Which clearly would be simplified quite a lot if [b,bd]=1 was substituted, Does anyone know how to this ? or could anyone point me to another tool capable of doing this?

Answer 1

The key trick is always to use the commutation relation to move either 'a' or 'a^\dagger' to the left, until you cannot do that anymore.

I don't have a good Sympy answer, but since you asked about other tools, here's a shameless plug about how you do it in Cadabra (https://cadabra.science) (which uses Sympy for various things, though not this particular computation). First setup the two sets of creation/annihilation operators using:

{a_{0}, ad_{0}}::NonCommuting;
{a_{1}, ad_{1}}::NonCommuting;
{a_{0}, ad_{0}, a_{1}, ad_{1}}::SortOrder.

They'll print nicer with

\bar{#}::Accent;
ad_{n?}::LaTeXForm("a^\dagger",n?,"").

Your Hamiltonian:

H:= w_{1} a_{0} ad_{0} + w_{2} a_{1} ad_{1} + g ad_{0} a_{1} + \bar{g} ad_{1} a_{1};

The commutator you want to compute:

ex:= @(H) a_{0} - a_{0} @(H);

Just expanding this (without simplification using the [a,ad]=1 commutator) is done with

distribute(ex);
sort_product(ex);

where the 2nd line moves operators with different subscripts through each other, but keeps the order of operators with the same subscripts. Applying the commutator until the expression no longer changes:

converge(ex):
    substitute(ex, $a_{n?} ad_{n?} = ad_{n?} a_{n?} + 1$)
    distribute(ex)
;

to finally give '-a_0 w_1 - a_1 g'.