Density Operator and Matrix¶

Imports¶

In [1]:

from IPython.display import display

In [2]:

# TODO: there is a bug in density.py that is preventing this from working, uncomment to reproduce
# from sympy import init_printing
# init_printing(use_latex=True)

In [3]:

from sympy import *
from sympy.core.trace import Tr
from sympy.physics.quantum import *
from sympy.physics.quantum.density import *
from sympy.physics.quantum.spin import (
    Jx, Jy, Jz, Jplus, Jminus, J2,
    JxBra, JyBra, JzBra,
    JxKet, JyKet, JzKet,
)

Basic density operator¶

Create a density matrix using symbolic states:

In [4]:

psi = Ket('psi')
phi = Ket('phi')

In [5]:

d = Density((psi,0.5),(phi,0.5));
d

Out[5]:

Density((|psi>, 0.5),(|phi>, 0.5))

In [6]:

d.states()

Out[6]:

(|psi>, |phi>)

In [7]:

d.probs()

Out[7]:

(0.5, 0.5)

In [8]:

d.doit()

Out[8]:

0.5*|phi><phi| + 0.5*|psi><psi|

In [9]:

Dagger(d)

Out[9]:

Density((|psi>, 0.5),(|phi>, 0.5))

In [10]:

A = Operator('A')

In [11]:

d.apply_op(A)

Out[11]:

Density((A*|psi>, 0.5),(A*|phi>, 0.5))

Density operator for spin states¶

Now create a density operator using spin states:

In [12]:

up = JzKet(S(1)/2,S(1)/2)
down = JzKet(S(1)/2,-S(1)/2)

In [13]:

d2 = Density((up,0.5),(down,0.5)); d2

Out[13]:

Density((|1/2,1/2>, 0.5),(|1/2,-1/2>, 0.5))

In [14]:

represent(d2)

Out[14]:

Matrix([
[0.5,   0],
[  0, 0.5]])

In [15]:

d2.apply_op(Jz)

Out[15]:

Density((Jz*|1/2,1/2>, 0.5),(Jz*|1/2,-1/2>, 0.5))

In [16]:

qapply(_)

Out[16]:

Density((hbar*|1/2,1/2>/2, 0.5),(-hbar*|1/2,-1/2>/2, 0.5))

In [17]:

qapply((Jy*d2).doit())

Out[17]:

0.5*Jy*|1/2,-1/2><1/2,-1/2| + 0.5*Jy*|1/2,1/2><1/2,1/2|

Evaluate entropy of the density matrices¶

In [18]:

entropy(d2)

Out[18]:

log(2)/2

In [19]:

entropy(represent(d2))

Out[19]:

log(2)/2

In [20]:

entropy(represent(d2,format="numpy"))

Out[20]:

(0.69314718055994529-0j)

In [21]:

entropy(represent(d2,format="scipy.sparse"))

Out[21]:

(0.69314718055994529-0j)

Density operators with tensor products¶

In [22]:

A, B, C, D = symbols('A B C D',commutative=False)

t1 = TensorProduct(A,B,C)

d = Density([t1, 1.0])
d.doit()

t2 = TensorProduct(A,B)
t3 = TensorProduct(C,D)

d = Density([t2, 0.5], [t3, 0.5])
d.doit() 

Out[22]:

0.5*(A*Dagger(A))x(B*Dagger(B)) + 0.5*(C*Dagger(C))x(D*Dagger(D))

In [23]:

d = Density([t2+t3, 1.0])
d.doit() 

Out[23]:

1.0*(A*Dagger(A))x(B*Dagger(B)) + 1.0*(A*Dagger(C))x(B*Dagger(D)) + 1.0*(C*Dagger(A))x(D*Dagger(B)) + 1.0*(C*Dagger(C))x(D*Dagger(D))

Trace operators on density operators with spin states¶

In [24]:

d = Density([JzKet(1,1),0.5],[JzKet(1,-1),0.5]);
t = Tr(d);
t

Out[24]:

Tr(Density((|1,1>, 0.5),(|1,-1>, 0.5)))

In [25]:

t.doit()

Out[25]:

1.00000000000000

Partial Trace on density operators with mixed state¶

In [26]:

A, B, C, D = symbols('A B C D',commutative=False)

t1 = TensorProduct(A,B,C)

d = Density([t1, 1.0])
d.doit()

t2 = TensorProduct(A,B)
t3 = TensorProduct(C,D)

d = Density([t2, 0.5], [t3, 0.5])
d

Out[26]:

Density((AxB, 0.5),(CxD, 0.5))

In [27]:

tr = Tr(d,[1])
tr.doit()

Out[27]:

0.5*A*Dagger(A)*Tr(B*Dagger(B)) + 0.5*C*Dagger(C)*Tr(D*Dagger(D))

Partial trace on density operators with spin states¶

In [28]:

tp1 = TensorProduct(JzKet(1,1), JzKet(1,-1))

Trace out the 0 index:

In [29]:

d = Density([tp1,1]);
t = Tr(d,[0])
t

Out[29]:

Tr((|1,1>x|1,-1>, 1))

In [30]:

t.doit()

Out[30]:

|1,-1><1,-1|

Trace out the 1 index:

In [31]:

t = Tr(d,[1])
t

Out[31]:

Tr((|1,1>x|1,-1>, 1))

In [32]:

t.doit()

Out[32]:

|1,1><1,1|

Examples of `qapply()` on density matrices with spin states¶

In [33]:

psi = Ket('psi')
phi = Ket('phi')

u = UnitaryOperator()
d = Density((psi,0.5),(phi,0.5)); d

qapply(u*d)

Out[33]:

O*Density((|psi>, 0.5),(|phi>, 0.5))

In [34]:

up = JzKet(S(1)/2, S(1)/2)
down = JzKet(S(1)/2, -S(1)/2)
d = Density((up,0.5),(down,0.5))

uMat = Matrix([[0,1],[1,0]])
qapply(uMat*d)

Out[34]:

Matrix([
[                                          0, Density((|1/2,1/2>, 0.5),(|1/2,-1/2>, 0.5))],
[Density((|1/2,1/2>, 0.5),(|1/2,-1/2>, 0.5)),                                           0]])

Example of `qapply()` on density matrices with qubits¶

In [35]:

from sympy.physics.quantum.gate import UGate
from sympy.physics.quantum.qubit import Qubit

uMat = UGate((0,), Matrix([[0,1],[1,0]]))
d = Density([Qubit('0'),0.5],[Qubit('1'), 0.5])
d

Out[35]:

Density((|0>, 0.5),(|1>, 0.5))

In [36]:

#after applying Not gate
qapply(uMat*d)

Out[36]:

U((0,),Matrix([
[0, 1],
[1, 0]]))*Density((|0>, 0.5),(|1>, 0.5))