# Density Operator and Matrix 

## Imports


```python
from IPython.display import display
```


```python
# 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)
```


```python
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:


```python
psi = Ket('psi')
phi = Ket('phi')
```


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




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




```python
d.states()
```




    (|psi>, |phi>)




```python
d.probs()
```




    (0.5, 0.5)




```python
d.doit()
```




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




```python
Dagger(d)
```




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




```python
A = Operator('A')
```


```python
d.apply_op(A)
```




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



## Density operator for spin states

Now create a density operator using spin states:


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


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




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




```python
represent(d2)
```




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




```python
d2.apply_op(Jz)
```




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




```python
qapply(_)
```




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




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




    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


```python
entropy(d2)
```




    log(2)/2




```python
entropy(represent(d2))
```




    log(2)/2




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




    (0.69314718055994529-0j)




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




    (0.69314718055994529-0j)



## Density operators with tensor products


```python
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() 
```




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




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




    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


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




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




```python
t.doit()
```




    1.00000000000000



## Partial Trace on density operators with mixed state


```python
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
```




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




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




    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


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

Trace out the `0` index:


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




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




```python
t.doit()
```




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



Trace out the `1` index:


```python
t = Tr(d,[1])
t
```




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




```python
t.doit()
```




    |1,1><1,1|



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


```python
psi = Ket('psi')
phi = Ket('phi')

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

qapply(u*d)
```




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




```python
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)
```




    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


```python
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
```




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




```python
#after applying Not gate
qapply(uMat*d)
```




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



bd0b0a21c066c3979424793b744f88d6214a8f24
ipynb
Jupyter Notebook
notebooks/density.ipynb
gvvynplaine/quantum_notebooks
58783823596465fe2d6c494c2cc3a53ae69a9752
BSD-3-Clause
2017-10-17T22:44:27.000Z
2022-03-28T06:26:46.000Z
notebooks/density.ipynb
gvvynplaine/quantum_notebooks
58783823596465fe2d6c494c2cc3a53ae69a9752
BSD-3-Clause
2017-10-09T05:16:41.000Z
2018-09-22T03:08:29.000Z
notebooks/density.ipynb
gvvynplaine/quantum_notebooks
58783823596465fe2d6c494c2cc3a53ae69a9752
BSD-3-Clause
2017-10-09T04:22:19.000Z
2022-03-28T06:25:21.000Z
Qwen/Qwen-72B
1. YES
2. YES