$\newcommand{\ket}[1]{\left|{#1}\right\rangle}$ $\newcommand{\bra}[1]{\left\langle{#1}\right|}$ $\newcommand{\braket}[2]{\left\langle{#1}\middle|{#2}\right\rangle}$ $\newcommand{\ketbra}[2]{\left|{#1}\right\rangle \left\langle{#2}\right|}$

Aplicações do operador densidade

Neste notebook, usaremos os conceitos aprendidos até agora, incluindo as ferramentas computacionais, em algumas aplicações físicas. Aproveitaremos, também, para extender algumas ideias, no sentido de descrever sistemas mais complexos, envolvendo múltiplas partículas.

import sympy as sy
import numpy as np
import sympy.physics as phys
from sympy import symbols, Matrix, sin, cos, exp, diff, integrate

sx = phys.matrices.msigma(1)  # sigma_x
sy = phys.matrices.msigma(2)  # sigma_y
sz = phys.matrices.msigma(3)  # sigma_z
I = Matrix([[1,0],[0,1]])     # unit matrix

Exemplo 1: sistemas compostos

Vamos considerar um sistema composto de duas partículas (ou dois sistemas físicos quaisquer) que, por simplicidade, vamos supor ser de dois níveis. Sejam A e B dois sistemas fechados de dois níveis representados pelas matrizes $\rho_A$ e $\rho_B$:

$$ \rho_A = \begin{pmatrix} 1/3 & 0 \\ 0 & 2/3 \end{pmatrix}; \quad \rho_B = \begin{pmatrix} 3/8 & 0 \\ 0 & 5/8 \end{pmatrix}, $$

:: ¡Desafio! :: _Para praticar, tente representar $\rho_A$ e $\rho_B$ como operadores..._

A matriz densidade do sistema composto é dada pelo produto tensorial

$$ \begin{aligned} \rho_{_{AB}} &=\rho_A \otimes \rho_B \\ \rho_A \otimes \rho_B &= \begin{pmatrix} 1/3 & 0 \\ 0 & 2/3 \end{pmatrix} \otimes \begin{pmatrix} 3/8 & 0 \\ 0 & 5/8 \end{pmatrix} \\ &= \begin{pmatrix} 1/8 & 0 & 0 & 0 \\ 0 & 5/24 & 0 & 0 \\ 0 & 0 & 1/4 & 0 \\ 0 & 0 & 0 & 5/12 \end{pmatrix} \end{aligned} $$

Como podemos facilmente verificar abaixo, usando sympy e numpy...

rho_A = Matrix([['1/3', '0'],['0', '2/3']])
rho_B = Matrix([['3/8', '0'],['0', '5/8']])

rho_A

$\displaystyle \left[\begin{matrix}\frac{1}{3} & 0\\0 & \frac{2}{3}\end{matrix}\right]$

rho_B

$\displaystyle \left[\begin{matrix}\frac{3}{8} & 0\\0 & \frac{5}{8}\end{matrix}\right]$

Oserve que o produto $\rho_A^2$ produz o resultado esperado..

rho_A*rho_A

$\displaystyle \left[\begin{matrix}\frac{1}{9} & 0\\0 & \frac{4}{9}\end{matrix}\right]$

(rho_A*rho_A).trace()

$\displaystyle \frac{5}{9}$

Assim como, usando os métodos definidos para o objeto Matrix, podemos verificar que $\rho_A \cdot \rho_A^{-1} = I$

rho_A*(rho_A.inv())

$\displaystyle \left[\begin{matrix}1 & 0\\0 & 1\end{matrix}\right]$

Usando a biblioteca sympy, calculamos o produto tensorial com TensorProduct( )

from sympy.physics.quantum import TensorProduct, Dagger
rho_AB = TensorProduct(rho_A, rho_B)
rho_AB

$\displaystyle \left[\begin{matrix}\frac{1}{8} & 0 & 0 & 0\\0 & \frac{5}{24} & 0 & 0\\0 & 0 & \frac{1}{4} & 0\\0 & 0 & 0 & \frac{5}{12}\end{matrix}\right]$

Neste exemplo, para facilitar a visualização e identificação das operações do produto tensorial, foram usadas matrizes diagonais (estados mistos, como mostra o $\text{Tr}(\rho_A^2)$, calculado abaixo)

rho_AB.trace()

$\displaystyle 1$

(rho_AB**2).trace()

$\displaystyle \frac{85}{288}$

mas poderíamos fazer o mesmo para qualquer matrix densidade, de qualquer dimensão. Mas adiante veremos outros exemplos, que exploram mais essa opções.

Isso também poderia ser feito, facilmente, usando a biblioteca numpy.

A = np.matrix([[1/3, 0],[0, 2/3]])
B = np.matrix([[3/8, 0],[0, 5/8]])
AB = np.kron(A,B)
AB

matrix([[0.125     , 0.        , 0.        , 0.        ],
        [0.        , 0.20833333, 0.        , 0.        ],
        [0.        , 0.        , 0.25      , 0.        ],
        [0.        , 0.        , 0.        , 0.41666667]])

AB.trace()

matrix([[1.]])

print(AB**2)
(AB**2).trace()

[[0.015625   0.         0.         0.        ]
 [0.         0.04340278 0.         0.        ]
 [0.         0.         0.0625     0.        ]
 [0.         0.         0.         0.17361111]]

matrix([[0.29513889]])

Em geral, nos casos que envolvem apenas cálculos numéricos e desempenho é importante, é sempre preferível usar a biblioteca numpy.

Exercício sugerido

Para praticar um pouco mais tudo que aprendemos até agora, expressse todos os operadores densidade acima usando a notação de operadores, ao invés da representação matricial.

Exemplo 2: operador reduzido e Traço parcial

Suponha, apenas para fixar ideia, que o sistema do Exemplo 1 representasse duas partículas de spin-1/2, e que inicialmente nos fosse dado apenas a matriz densidade do sistema composto AB.

Essa informação permite-nos, facilmente, calcular as probilidades dos estados produto. Por exemplo...

Faça você mesmo: Como poderíamos calcular a probabilidade de observar o estado $\ket{\downarrow\uparrow}$?

#Dica: a resposta é 1/4... 
#Tente calcular você mesmo, antes de resolvermos juntos na aula

Uma sitação diferente, mas ainda simples, é perguntar, para os estados de duas partículas, a probabilidade de observar uma medida onde a partícula B está no estado $\ket{\downarrow}$. Essa medida também pode ser respondida a pela observação direta da matriz densidade total do sistema.

DIY: Tente fazer sua previsão, antes de resolvermos...

Porém, se tivéssemos apenas a matriz total, como poderíamos responder, em geral, as perguntas abaixo?

No sistema do Exemplo 1, se fosse feita uma medida na partícula B para determinar a projeção do spin na direção $\hat z$, qual seria a probabilidade de medir o estado $\ket{\downarrow}$?

E se a medida fosse na direção $\hat x_B$, qual seria a probabilidade de observar o estado $\ket{\uparrow;x}$? Qual seria o valor esperado de $\langle S_x \rangle$ para cada partícula, se medidas separadamente?

DIY: Tente fazer sua previsão, antes de resolvermos.

Para perguntas desse tipo, podemos usar a matriz densidade reduzida do sistema composto (i.e., as densidades $\rho_A$ e $\rho_B$), que podem ser calculada a partir da matriz densidade total usando o conceito de Traço parcial

$$ \rho_A = \text{Tr}_B\,(\rho_{_{AB}}) $$

onde o símbolo $\text{Tr}_B(\cdot)$ significa, em palavras, fazer o traço apenas sobre o espaço de Hilbert do sistema B. Na prática, o traço parcial é um mapa linear definido abaixo

Traço parcial

O traço parcial pode ser definido por

$$\text{Tr}_B\,(\,\ketbra{a}{a} \otimes \ketbra{b}{b}\,) = \ketbra{a}{a}\,\text{Tr}(\,\ketbra{b}{b}\,),$$

portanto,

$$\text{Tr}_B\,(\,\rho_A \otimes \rho_B\,) = \rho_A \,\text{Tr}(\,\rho_B \,) = \rho_A $$

pois $\text{Tr}(\rho_A)=1$, como qualquer operador densidade.

---

?!🤔?:: Ok, mas como calcular na prática a partir da matriz total??

---

¡DIY!: Pense um pouco, e tente fazer sozinho(a) antes de resolvermos juntos.

sys_A = TensorProduct(rho_A,I)
sys_B = TensorProduct(I,rho_B)

sys_A

$\displaystyle \left[\begin{matrix}\frac{1}{3} & 0 & 0 & 0\\0 & \frac{1}{3} & 0 & 0\\0 & 0 & \frac{2}{3} & 0\\0 & 0 & 0 & \frac{2}{3}\end{matrix}\right]$

sys_B

$\displaystyle \left[\begin{matrix}\frac{3}{8} & 0 & 0 & 0\\0 & \frac{5}{8} & 0 & 0\\0 & 0 & \frac{3}{8} & 0\\0 & 0 & 0 & \frac{5}{8}\end{matrix}\right]$

sys_A*sys_B

$\displaystyle \left[\begin{matrix}\frac{1}{8} & 0 & 0 & 0\\0 & \frac{5}{24} & 0 & 0\\0 & 0 & \frac{1}{4} & 0\\0 & 0 & 0 & \frac{5}{12}\end{matrix}\right]$

Ou seja, temo o resultado (esperado):

$$ (\rho_A \otimes \hat I) \cdot (\hat I \otimes \rho_B) = \rho_{_AB} \\ S_A = \rho_A \otimes \hat I \\ S_B = \hat I \otimes \rho_B $$

Temos o produto das matrizes $S_A$ e $S_B$. Note que podermos fazer

$$ S_B = S_A^{-1} \cdot S_A \cdot S_B $$

sys_A.inv()*sys_A*sys_B

$\displaystyle \left[\begin{matrix}\frac{3}{8} & 0 & 0 & 0\\0 & \frac{5}{8} & 0 & 0\\0 & 0 & \frac{3}{8} & 0\\0 & 0 & 0 & \frac{5}{8}\end{matrix}\right]$

Aplicações nos problemas propostos

Vamos agora aplicar esse conceito nos problemas (perguntas) proposto acima.

• Qual a probabilidade de encontrar o sistema B (segunda partícula) no estado $\ket{\downarrow}$

  O projetor do estado $\ket{\downarrow}_B$ é $\hat P_d=\ketbra{\downarrow}{\downarrow}$

up = Matrix([[1],[0]])
down = Matrix([[0],[1]])
pd = TensorProduct(down.T,down)
op1 = TensorProduct(I,pd)

up

$\displaystyle \left[\begin{matrix}1\\0\end{matrix}\right]$

down

$\displaystyle \left[\begin{matrix}0\\1\end{matrix}\right]$

pd

$\displaystyle \left[\begin{matrix}0 & 0\\0 & 1\end{matrix}\right]$

op1

$\displaystyle \left[\begin{matrix}0 & 0 & 0 & 0\\0 & 1 & 0 & 0\\0 & 0 & 0 & 0\\0 & 0 & 0 & 1\end{matrix}\right]$

(rho_AB*op1)

$\displaystyle \left[\begin{matrix}0 & 0 & 0 & 0\\0 & \frac{5}{24} & 0 & 0\\0 & 0 & 0 & 0\\0 & 0 & 0 & \frac{5}{12}\end{matrix}\right]$

(rho_AB*op1).trace()

$\displaystyle \frac{5}{8}$

rho_AB[1,1]+rho_AB[3,3]

$\displaystyle \frac{5}{8}$

• Qual o valor esperado de $\langle S_x \rangle$ para a partícula B?

  O operador $\hat S_x $ é dado por $$S_x = \frac{\hbar}{2} \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix} $$

sx  # a menos de uma constante multiplicativa

$\displaystyle \left[\begin{matrix}0 & 1\\1 & 0\end{matrix}\right]$

op2 = TensorProduct(I,sx)
op2

$\displaystyle \left[\begin{matrix}0 & 1 & 0 & 0\\1 & 0 & 0 & 0\\0 & 0 & 0 & 1\\0 & 0 & 1 & 0\end{matrix}\right]$

Assim, podemos facilmente calcular o valor esperado $\langle S_x \rangle_B$

$$ \langle S_x \rangle_B = \text{Tr}\,(\rho\,(\hat I \otimes S_x)) $$

(rho_AB*op2)

$\displaystyle \left[\begin{matrix}0 & \frac{1}{8} & 0 & 0\\\frac{5}{24} & 0 & 0 & 0\\0 & 0 & 0 & \frac{1}{4}\\0 & 0 & \frac{5}{12} & 0\end{matrix}\right]$

(rho_AB*op2).trace()

$\displaystyle 0$

Portanto, o $\langle S_x \rangle_B=0$, como esperado.

• Qual a probabilidade de medir o estado $\ket{\uparrow;x}$ numa medida $\hat x_B$ ?

¡Faça você!: Pense um pouco, e tente fazer sozinho(a) essa parte...

#Dica: comece expressando o estado (up)_x na base original...

Exemplo 3: usando o QuTiP

Nos exemplos acima utilizamos o sympy, e mostramos como calcular o produto tensorial (Kronecker product) usando o numpy. Claro que poderíamos também, facilmente, fazer tudo numericamente, no numpy, dado que temos valores numéricos apenas, mas ao invés de repetir isso aqui, eu vou mostrar como fazer essas operações usando a outra ferramente que eu apresentei para vocês: o QuTiP (Quantum Toolbox in Python).

import qutip as qt
r_A = qt.Qobj(A)  # utilizando as matrizes A e B já definidas acima (numpy)
r_B = qt.Qobj(B)

A

matrix([[0.33333333, 0.        ],
        [0.        , 0.66666667]])

r_A

Quantum object: dims = [[2], [2]], shape = (2, 2), type = oper, isherm = True\begin{equation*}\left(\begin{array}{*{11}c}0.333 & 0.0\\0.0 & 0.667\\\end{array}\right)\end{equation*}

r_B

Quantum object: dims = [[2], [2]], shape = (2, 2), type = oper, isherm = True\begin{equation*}\left(\begin{array}{*{11}c}0.375 & 0.0\\0.0 & 0.625\\\end{array}\right)\end{equation*}

rho = qt.tensor(r_A,r_B)  # produto tensorial 
rho

Quantum object: dims = [[2, 2], [2, 2]], shape = (4, 4), type = oper, isherm = True\begin{equation*}\left(\begin{array}{*{11}c}0.125 & 0.0 & 0.0 & 0.0\\0.0 & 0.208 & 0.0 & 0.0\\0.0 & 0.0 & 0.250 & 0.0\\0.0 & 0.0 & 0.0 & 0.417\\\end{array}\right)\end{equation*}

rho.tr()

0.9999999999999999

Neste caso, para obter o traço parcial $\text{Tr}_B(\cdot)$, basta fazer

rho.ptrace(0)  # traço parcial Tr_B(), para obter operador da segunda partícula

Quantum object: dims = [[2], [2]], shape = (2, 2), type = oper, isherm = True\begin{equation*}\left(\begin{array}{*{11}c}0.333 & 0.0\\0.0 & 0.667\\\end{array}\right)\end{equation*}

Note que o mesmo não funcionaria partindo da matriz densidade completa (sem o produto tensorial).

Mostrarei isso usando a matriz AB, já definida acima.

AB  # matriz definida acima (usando numpy)

matrix([[0.125     , 0.        , 0.        , 0.        ],
        [0.        , 0.20833333, 0.        , 0.        ],
        [0.        , 0.        , 0.25      , 0.        ],
        [0.        , 0.        , 0.        , 0.41666667]])

r = qt.Qobj(AB)
r

Quantum object: dims = [[4], [4]], shape = (4, 4), type = oper, isherm = True\begin{equation*}\left(\begin{array}{*{11}c}0.125 & 0.0 & 0.0 & 0.0\\0.0 & 0.208 & 0.0 & 0.0\\0.0 & 0.0 & 0.250 & 0.0\\0.0 & 0.0 & 0.0 & 0.417\\\end{array}\right)\end{equation*}

Preste atenção nas características dos dois objetos Qobj(), em particular em dims.

Eles são, efetivamente, objetos diferentes, como podemos ver tentado tomar o traço parcial

r.ptrace(1)  # este erro é de propósito, para mostrar a dif. dos objetos envolvidos.

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-54-300a3eebe6c5> in <module>
----> 1 r.ptrace(1)  # este erro é de propósito, para mostrar a dif. dos objetos envolvidos.

~\Anaconda3\lib\site-packages\qutip\qobj.py in ptrace(self, sel, sparse)
   1353             return q.tidyup() if settings.auto_tidyup else q
   1354         else:
-> 1355             return _ptrace_dense(self, sel)
   1356 
   1357     def permute(self, order):

~\Anaconda3\lib\site-packages\qutip\qobj.py in _ptrace_dense(Q, sel)
   2192         sel = np.asarray(sel)
   2193     sel = list(np.sort(sel))
-> 2194     dkeep = (rd[sel]).tolist()
   2195     qtrace = list(set(np.arange(nd)) - set(sel))
   2196     dtrace = (rd[qtrace]).tolist()

IndexError: index 1 is out of bounds for axis 0 with size 1

Vamos calcular os mesmos resultados usando o QuTiP...

I2 = qt.qeye(2)
up = qt.Qobj(np.matrix([[1],[0]]))
down = qt.Qobj(np.matrix([[0],[1]]))
pd = qt.tensor(down.dag(),down)
op1 = qt.tensor(I2,pd)

down

Quantum object: dims = [[2], [1]], shape = (2, 1), type = ket\begin{equation*}\left(\begin{array}{*{11}c}0.0\\1.0\\\end{array}\right)\end{equation*}

down.dag()

Quantum object: dims = [[1], [2]], shape = (1, 2), type = bra\begin{equation*}\left(\begin{array}{*{11}c}0.0 & 1.0\\\end{array}\right)\end{equation*}

I2

Quantum object: dims = [[2], [2]], shape = (2, 2), type = oper, isherm = True\begin{equation*}\left(\begin{array}{*{11}c}1.0 & 0.0\\0.0 & 1.0\\\end{array}\right)\end{equation*}

a=I*up  # cuidado ao multiplicar "objetos" diferentes!
a       # embora o resultado possa parecer o mesmo...

$\displaystyle \left[\begin{matrix}1.0\\0\end{matrix}\right]$

a.dag()  # eles não são iguais!! Informação foi perdida no processo!

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-63-dff7409eea7d> in <module>
----> 1 a.dag()

AttributeError: 'MutableDenseMatrix' object has no attribute 'dag'

pd   # preste atenção nas informações (abaixo) sobre esse objeto

Quantum object: dims = [[1, 2], [2, 1]], shape = (2, 2), type = oper, isherm = True\begin{equation*}\left(\begin{array}{*{11}c}0.0 & 0.0\\0.0 & 1.0\\\end{array}\right)\end{equation*}

op1   # preste atenção nas informações (abaixo) sobre esse objeto

Quantum object: dims = [[2, 1, 2], [2, 2, 1]], shape = (4, 4), type = oper, isherm = True\begin{equation*}\left(\begin{array}{*{11}c}0.0 & 0.0 & 0.0 & 0.0\\0.0 & 1.0 & 0.0 & 0.0\\0.0 & 0.0 & 0.0 & 0.0\\0.0 & 0.0 & 0.0 & 1.0\\\end{array}\right)\end{equation*}

rho  # preste atenção nas informações (abaixo) sobre esse objeto

Quantum object: dims = [[2, 2], [2, 2]], shape = (4, 4), type = oper, isherm = True\begin{equation*}\left(\begin{array}{*{11}c}0.125 & 0.0 & 0.0 & 0.0\\0.0 & 0.208 & 0.0 & 0.0\\0.0 & 0.0 & 0.250 & 0.0\\0.0 & 0.0 & 0.0 & 0.417\\\end{array}\right)\end{equation*}

rho*op1  # este erro é de propósito, para mostrar a dif. dos objetos envolvidos.

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-67-3670142e2a67> in <module>
----> 1 rho*op1  # este erro é de propósito, para mostrar a dif. dos objetos envolvidos.

~\Anaconda3\lib\site-packages\qutip\qobj.py in __mul__(self, other)
    553 
    554             else:
--> 555                 raise TypeError("Incompatible Qobj shapes")
    556 
    557         elif isinstance(other, np.ndarray):

TypeError: Incompatible Qobj shapes

Ooops! O que aconteceu aqui?? 😲

(Dica: objserve os tipos de objetos que forma multiplicados...)

Ok, vamos fazer da maneira certa, desta vez. Para isso usaremos a função que cria a matriz densidade diretamente a partir do vetor de estado, como mostrado abaixo

pd = qt.ket2dm(down)  # forma segura de construir o op. dens. de um ket
pd          # preste atenção nas informações (abaixo) sobre esse objeto

Quantum object: dims = [[2], [2]], shape = (2, 2), type = oper, isherm = True\begin{equation*}\left(\begin{array}{*{11}c}0.0 & 0.0\\0.0 & 1.0\\\end{array}\right)\end{equation*}

op1 = qt.tensor(I2,pd)
op1        # preste atenção nas informações (abaixo) sobre esse objeto

Quantum object: dims = [[2, 2], [2, 2]], shape = (4, 4), type = oper, isherm = True\begin{equation*}\left(\begin{array}{*{11}c}0.0 & 0.0 & 0.0 & 0.0\\0.0 & 1.0 & 0.0 & 0.0\\0.0 & 0.0 & 0.0 & 0.0\\0.0 & 0.0 & 0.0 & 1.0\\\end{array}\right)\end{equation*}

rho*op1

Quantum object: dims = [[2, 2], [2, 2]], shape = (4, 4), type = oper, isherm = True\begin{equation*}\left(\begin{array}{*{11}c}0.0 & 0.0 & 0.0 & 0.0\\0.0 & 0.208 & 0.0 & 0.0\\0.0 & 0.0 & 0.0 & 0.0\\0.0 & 0.0 & 0.0 & 0.417\\\end{array}\right)\end{equation*}

Sucesso!! 😃

Agora basta tomar o traço, como fizemos antes.

(rho*op1).tr()

0.625

Na dúvida, vamos verificar se é o mesmo resultado...

5/8 == 0.625

True

Qual o valor esperado de $\langle S_z \rangle$ no estado $\rho$ (sistema total)

sx = qt.sigmax()
sy = qt.sigmay()
sz = qt.sigmaz()
sz

Quantum object: dims = [[2], [2]], shape = (2, 2), type = oper, isherm = True\begin{equation*}\left(\begin{array}{*{11}c}1.0 & 0.0\\0.0 & -1.0\\\end{array}\right)\end{equation*}

szA = qt.tensor(sz,I2)
szB = qt.tensor(I2,sz)
Sz = qt.tensor(sz,sz)
Sz

Quantum object: dims = [[2, 2], [2, 2]], shape = (4, 4), type = oper, isherm = True\begin{equation*}\left(\begin{array}{*{11}c}1.0 & 0.0 & 0.0 & 0.0\\0.0 & -1.0 & 0.0 & 0.0\\0.0 & 0.0 & -1.0 & 0.0\\0.0 & 0.0 & 0.0 & 1.0\\\end{array}\right)\end{equation*}

rho*Sz

Quantum object: dims = [[2, 2], [2, 2]], shape = (4, 4), type = oper, isherm = True\begin{equation*}\left(\begin{array}{*{11}c}0.125 & 0.0 & 0.0 & 0.0\\0.0 & -0.208 & 0.0 & 0.0\\0.0 & 0.0 & -0.250 & 0.0\\0.0 & 0.0 & 0.0 & 0.417\\\end{array}\right)\end{equation*}

(rho*Sz).tr()

0.08333333333333331

Outra forma de obter o mesmo resultado é usando a função qutip.expect( oper, state ).

qt.expect(Sz,rho)

0.08333333333333331

Exemplo 4: estado emaranhado

Vamos explorar um pouco mais o conceito de matriz densidade total e reduzida em um sistema composto de duas partículas de spin-1/2, mas agora considerando estados puros e superposições.

Vamos considerar os estados $\ket{\uparrow \uparrow}$ e $\ket{\downarrow \downarrow}$ e sua superposição

uu = qt.tensor(up,up)      # up, up
dd = qt.tensor(down,down)  # down, down 
uu  

Quantum object: dims = [[2, 2], [1, 1]], shape = (4, 1), type = ket\begin{equation*}\left(\begin{array}{*{11}c}1.0\\0.0\\0.0\\0.0\\\end{array}\right)\end{equation*}

dd

Quantum object: dims = [[2, 2], [1, 1]], shape = (4, 1), type = ket\begin{equation*}\left(\begin{array}{*{11}c}0.0\\0.0\\0.0\\1.0\\\end{array}\right)\end{equation*}

# superposiçao
b1 = (uu + dd).unit()
b1

Quantum object: dims = [[2, 2], [1, 1]], shape = (4, 1), type = ket\begin{equation*}\left(\begin{array}{*{11}c}0.707\\0.0\\0.0\\0.707\\\end{array}\right)\end{equation*}

O estado $\ket{b_1}$ é dado por

$$ \ket{b_1} = \frac{ \ket{\uparrow \uparrow} + \ket{\downarrow \downarrow} }{\sqrt2} $$

que é um estado puro.

Vejamos como fica a matriz densidade desse estado de duas partículas

dm_b1 = qt.ket2dm(b1)
dm_b1

Quantum object: dims = [[2, 2], [2, 2]], shape = (4, 4), type = oper, isherm = True\begin{equation*}\left(\begin{array}{*{11}c}0.500 & 0.0 & 0.0 & 0.500\\0.0 & 0.0 & 0.0 & 0.0\\0.0 & 0.0 & 0.0 & 0.0\\0.500 & 0.0 & 0.0 & 0.500\\\end{array}\right)\end{equation*}

(dm_b1**2).tr()

0.9999999999999996

Podemos determinar agora as matrizes densidades de cada partícula, usando o traço parcial...

dm_A = dm_b1.ptrace(0)
dm_A

Quantum object: dims = [[2], [2]], shape = (2, 2), type = oper, isherm = True\begin{equation*}\left(\begin{array}{*{11}c}0.500 & 0.0\\0.0 & 0.500\\\end{array}\right)\end{equation*}

dm_B = dm_b1.ptrace(1)
dm_B

Quantum object: dims = [[2], [2]], shape = (2, 2), type = oper, isherm = True\begin{equation*}\left(\begin{array}{*{11}c}0.500 & 0.0\\0.0 & 0.500\\\end{array}\right)\end{equation*}

Esse é um resultado surpreendente, pois as matrizes reduzidas corresponde a estados de máxima mistura, como podemos confirmar com o traço de $\rho_A^2$

dm_A**2

Quantum object: dims = [[2], [2]], shape = (2, 2), type = oper, isherm = True\begin{equation*}\left(\begin{array}{*{11}c}0.250 & 0.0\\0.0 & 0.250\\\end{array}\right)\end{equation*}

(dm_A**2).tr()

0.4999999999999998

Esse resultado ocorre devido ao fenômeno de emaranhamento, presente no estado $\ket{b_1}$, que corresponde a um dos chamados estados de Bell.

B1 = qt.bell_state()
B1

Quantum object: dims = [[2, 2], [1, 1]], shape = (4, 1), type = ket\begin{equation*}\left(\begin{array}{*{11}c}0.707\\0.0\\0.0\\0.707\\\end{array}\right)\end{equation*}

b1 == B1

True

Tente explorar, agora, usando o que aprendeu em um estado de superposição que não seja emaranhado e veja como são as matrizes densidades de cada partícula...

Farei um exemplo, usando um estado puro sem superposição, para você aprender o caminho.

dm_uu = qt.ket2dm(uu)
dm_uu

Quantum object: dims = [[2, 2], [2, 2]], shape = (4, 4), type = oper, isherm = True\begin{equation*}\left(\begin{array}{*{11}c}1.0 & 0.0 & 0.0 & 0.0\\0.0 & 0.0 & 0.0 & 0.0\\0.0 & 0.0 & 0.0 & 0.0\\0.0 & 0.0 & 0.0 & 0.0\\\end{array}\right)\end{equation*}

dm_uu.ptrace(0)

Quantum object: dims = [[2], [2]], shape = (2, 2), type = oper, isherm = True\begin{equation*}\left(\begin{array}{*{11}c}1.0 & 0.0\\0.0 & 0.0\\\end{array}\right)\end{equation*}

dm_uu.ptrace(1)

Quantum object: dims = [[2], [2]], shape = (2, 2), type = oper, isherm = True\begin{equation*}\left(\begin{array}{*{11}c}1.0 & 0.0\\0.0 & 0.0\\\end{array}\right)\end{equation*}

Note que agora a situação é bem diferente!

(dm_uu.ptrace(0)**2)

Quantum object: dims = [[2], [2]], shape = (2, 2), type = oper, isherm = True\begin{equation*}\left(\begin{array}{*{11}c}1.0 & 0.0\\0.0 & 0.0\\\end{array}\right)\end{equation*}

(dm_uu.ptrace(0)**2).tr()

1.0

Confirmando que trata-se de um estado puro!

Tente estender, agora, a ideia para outros estados de superposição e observe as matrizes densidades reduzidas de cada uma das partículas, verificando, ao final, se são estados puros ou mistos. Por exemplo, fazendo o exercício sugerido abaixo.

Exercício sugerido

Repita os passos acima para o estado $$ \ket{\psi} = \frac{ \ket{\uparrow \uparrow} + \ket{\downarrow \uparrow} }{\sqrt2} $$ Como seriam as matrizes densidades total e reduzidas em outra base, por exemplo, na base $\hat x$ ou $\hat y$?
Algumas delas tem temos fora da diagonal?

# Tente fazer em casa, usando este nobebook como modelo...
# explore também outras situações (por exemplo, p/ spins interagentes)
# "Enjoy & have fun!"