In sympy, simplify an expression using a canonical commutation relation

Question

I have a ladder operator â, which satisfies this commutator relation with its own adjoint:

[â, â⁺] = 1

In sympy I have written this code:

import sympy
from sympy import *
from sympy.physics.quantum import *

a = Operator('a')
ad = Dagger(a)

ccr = Eq( Commutator(a, ad),  1 )

Now I need to expand and simplify an expression like this:

(â⁺ + â)⁴

If I just use ((ad + a)**4).expand(), sympy doesn't use the commutator relation. How do I simplify the expression while using the canonical commutator relation?

Answer 1

Inspired by m93a's answer, I wrote an implementation that simultaneously handles an arbitrary number of commutation relations. I needed this because I was using multiple bosonic operators.

from sympy import *

def expand_powers(args):
    
    args = list(args)
    if len(args) == 0:
        return []

    rest = [] if len(args) == 1 else args[1:]

    if isinstance(power := args[0], Pow) and not isinstance(power.base, Number):
        return [power.base for i in range(power.exp)] + expand_powers(rest)

    expanded = [args[0]] + expand_powers(rest)
    return expanded
            
        

def qsimplify(expr, commutator_relations=None):
    
    commutator_relations = [] if commutator_relations is None else commutator_relations

    
    
    def qsimplify_walktree(expr):

        if isinstance(expr, Number):
            return expr

        if not isinstance(expr, AssocOp) and not isinstance(expr, Function):
            return expr.copy()    

        elif not isinstance(expr, Mul):
            return expr.func(*(qsimplify_walktree(node) for node in expr.args)) 

        args = expand_powers(list(expr.args))
        for i in range(len(args)-1):


            A = args[i]
            B = args[i+1]

            for pair, value in ordered_commutator_relations:
                if (B, A) == pair:
                    args = args[0:i] + [B*A - value] + args[i+2:]
                    result = Mul(*args).expand()
                    return qsimplify_walktree(result)

        return expr.copy()
    
    ordered_commutator_relations = []
    for relation in commutator_relations:
        lhs, rhs = relation
        for idx, node in enumerate(nodes := list(preorder_traversal(lhs))):
            if isinstance(node, Commutator):
                comm = node
                if idx == 0:
                    ordered_commutator_relations.append((comm.args, rhs))
                    break

                if isinstance(sign := nodes[idx-1], Number):
                    ordered_commutator_relations.append((comm.args, -rhs))
                    break
                    
    return qsimplify_walktree(expr)

Example usage (writing the commutators in this order leads to normal-ordered expressions):

ai = Operator('{a}_{i}')
adi = Dagger(ai)


aj = Operator('{a}_{j}')
adj = Dagger(aj)


commutator_relations = [(Commutator(ai, Dagger(ai)), 1),
                        (Commutator(aj, Dagger(aj)), 1),
                        (Commutator(ai, aj), 0),
                        (Commutator(Dagger(ai), aj), 0),
                        (Commutator(ai, Dagger(aj)), 0),
                        (Commutator(Dagger(ai), Dagger(aj)), 0)
                       ]

((Dagger(ai) + ai + Dagger(aj))**4).expand().qsimplify(commutator_relations)

Answer 2

I couldn't find any built-in way to do it, so I wrote a very basic algorithm for it. It's used like this:

((ad + a)**4).expand().apply_ccr(ccr)

Result

3 + 12 a⁺ a + 4 a⁺ a³ + 6 a⁺² + 6 a⁺² a² + 4 a⁺³ a + a⁺⁴ + 6a² + a⁴

.

There is an optional argument called reverse which would rearange the expression to be a first and then a⁺. This is necessary to overcome the limitations of sympy which doesn't let you to specify the Commutator in a different order [source].

---

This is the implementation of apply_ccr:

from sympy.core.operations import AssocOp

def apply_ccr(expr, ccr, reverse=False):
    if not isinstance(expr, Basic):
        raise TypeError("The expression to simplify is not a sympy expression.")

    if not isinstance(ccr, Eq):
        if isinstance(ccr, Basic):
            ccr = Eq(ccr, 0)
        else:
            raise TypeError("The canonical commutation relation is not a sympy expression.")

    comm = None

    for node in preorder_traversal(ccr):
        if isinstance(node, Commutator):
            comm = node
            break

    if comm is None:
        raise ValueError("The cannonical commutation relation doesn not include a commutator.")

    solutions = solve(ccr, comm)

    if len(solutions) != 1:
        raise ValueError("There are more solutions to the cannonical commutation relation.")

    value = solutions[0]

    A = comm.args[0]
    B = comm.args[1]

    if reverse:
        (A, B) = (B, A)
        value = -value

    def is_expandable_pow_of(base, expr):
        return isinstance(expr, Pow) \
            and base == expr.args[0] \
            and isinstance(expr.args[1], Number) \
            and expr.args[1] >= 1


    def walk_tree(expr):
        if isinstance(expr, Number):
            return expr

        if not isinstance(expr, AssocOp) and not isinstance(expr, Function):
            return expr.copy()

        elif not isinstance(expr, Mul):
            return expr.func(*(walk_tree(node) for node in expr.args))

        else:
            args = [arg for arg in expr.args]

            for i in range(len(args)-1):
                x = args[i]
                y = args[i+1]

                if B == x and A == y:
                    args = args[0:i] + [A*B - value] + args[i+2:]
                    return walk_tree( Mul(*args).expand() )

                if B == x and is_expandable_pow_of(A, y):
                    ypow = Pow(A, y.args[1] - 1)
                    args = args[0:i] + [A*B - value, ypow] + args[i+2:]
                    return walk_tree( Mul(*args).expand() )

                if is_expandable_pow_of(B, x) and A == y:
                    xpow = Pow(B, x.args[1] - 1)
                    args = args[0:i] + [xpow, A*B - value] + args[i+2:]
                    return walk_tree( Mul(*args).expand() )

                if is_expandable_pow_of(B, x) and is_expandable_pow_of(A, y):
                    xpow = Pow(B, x.args[1] - 1)
                    ypow = Pow(A, y.args[1] - 1)
                    args = args[0:i] + [xpow, A*B - value, ypow] + args[i+2:]
                    return walk_tree( Mul(*args).expand() )

            return expr.copy()


    return walk_tree(expr)


Basic.apply_ccr = lambda self, ccr, reverse=False: apply_ccr(self, ccr, reverse)

(No rights reserved.)