Search code examples
operatorssympyphysicssymbolic-math

Substitute Pauli Operators in SymPy

I want to calculate the trace of a tensor product of Pauli matrices but I'm not sure if it is possible.

import sympy as sp

from sympy.physics.quantum import TensorProduct as ts

from sympy.physics.matrices import msigma
from sympy.physics.quantum import qapply

from sympy.physics.paulialgebra import Pauli
from sympy.physics.quantum import tensor_product_simp as tpsimp

sI , sx, sy ,sz = 1, Pauli(1), Pauli(2),Pauli(3)
Id, x, y, z = sp.eye(2), msigma(1), msigma(2), msigma(3)

convrt = {sI:Id, sx:x, sy:y, sz:z}

r = ts(sx,sy,sz,sI)
print(tpsimp(qapply(r.subs(convrt))).doit())

The result is:

enter image description here

But I expected a 16x16 Matrix:

enter image description here

With the result I got I can't calculate the trace of the matrix. Is it possible to do it? This is just an small example, my r isn't easy like this one, it involve functions and some constants and has more terms but my main problem is to convert properly tensor products of Pauli operators to a matrix

Solution

  • Looking at the documentation of the TensorProduct, it states:

    For matrices, this uses matrix_tensor_product to compute the Kronecker or tensor product matrix. For other objects a symbolic TensorProduct instance is returned.

    In particular, r is a tensor product of symbolic objects. When you substitute the matrices:

    expr = r.subs(convrt)

    This is a tensor product of matrices, but because it was originally a tensor product of symbolic object, the matrix_tensor_product was never called.

    I would consider opening an issue for this, and ask to implement an as_explicit method, just like what's found on matrix expressions.

    Usually, I would reconstruct the object (the tensor product in this case) with this syntax, in order to force the matrix_tensor_product to be executed:

    expr.func(*expr.args)

    However, this still doesn't evaluate the product. Note that your Pauli matrices are instances of Matrix (alias of MutableDenseMatrix). But when we substituted them into the symbolic tensor product, somewhere they were changed to ImmutableDenseMatrix. Sadly, the constructor of TensorProduct looks for instances of Matrix, hence ImmutableDenseMatrix are considered as symbolic objects. In my opinion, there are two bugs:

    1. The conversion from MutableDenseMatrix to ImmutableDenseMatrix shouldn't have happened.
    2. The constructor of TensorProduct should also look for instances of ImmutableDenseMatrix.

    All of this to say that your are forced to execute the undocumented matrix_tensor_product:

    from sympy.physics.quantum.matrixutils import matrix_tensor_product
matrix_tensor_product(*expr.args)