Bug in group_by_operators

[akhmerov]

Hello sympt team!

While implementing second quantization in pymablock, @iarayad and I believe that we saw a bug in your implementation of group_by_operators, which influences the example https://qcode-uni-a.github.io/sympt/examples/I1_Transmon_resonator/ or any other Hamiltonian where $H_0$ isn't linear in boson number. In those cases, perturbation theory generates terms like $a^\dagger (a^\dagger a + x)^{-1} a$. group_by_operators will process this term as $a^\dagger a: (a^\dagger a + x)^{-1}$, which is incorrect because the energy denominator does not commute with $a$. Because group_by_operators is used in all examples as a way to examine the results, it seems like a serious issue.

To be able to manipulate expressions with energy denominators, where normal ordering does not work, and also to be able to use sympy simplification routines with operator expressions, we developed a number ordering scheme. Please find its implementation over here:

https://github.com/quantum-tinkerer/pymablock/blob/745a68291e5ae5c1c6c4a1bbfad73c70cf26a20c/pymablock/second_quantization.py#L64-L327

[OtiDioti]

Thanks for flagging this issue. We looked into it and found that if you use an expression like $a^\dagger \left(a^\dagger a + x\right)^{-1} a$, then group_by_operators incorrectly groups it as you correctly commented.

That said, this doesn’t affect sympt’s output because sympt always reorders the terms—using commutation relations—to place any negative power of bosonic operators on the far left of any "positive power expression" of bosonic operators. This is handled via the MulGroup class which naturally keeps the two types of expressions separate. In other words, the problematic expression never actually appears when sympt processes a Hamiltonian.

Again, thank you for providing this insight. We will start working on this, so that users will be able to use group_by_operators not only for Hamiltonian expressions outputted by sympt, but also for custom expression that may contain terms like the one provided in your example (i.e. having bosonic denominators not on the very left of positive powers of bosonic expressions).

Please let us know if you find some other issue, or if there was a misunderstanding from our part of what is that you were trying to point out.

[akhmerov]

Ah, thanks for the clarification, good to know. I've looked at the MulGroup, and indeed I see where the shift is applied.

I agree that the bug is insignificant for the sympt workflow.