por
Franklin Marquezino
Universidade Federal do Rio de Janeiro
URL: http://www.cos.ufrj.br/~franklin
Em desenvolvimento!
Nesse arquivo vou descrever as matrizes que serão necessárias para nossa simulação. Para tornar o estudo mais fácil de acompanhar, vou começar pelo caso 1D e depois passar para o caso 2D antes de generalizar para outros grafos regulares. É possível visualizar este arquivo usando o nbviewer. No entanto, é altamente aconselhável instalar o IPython Notebook para rodar os exemplos. Aliás, o IPython Notebook pode ser uma ferramenta muito útil para pesquisadores!
Primeiramente, temos que importar algumas bibliotecas do Python. Estou usando o Python 3, mas creio que tudo (ou quase tudo) deve funcionar também no Python 2. Talvez você precise instalar as bibliotecas Matplotlib, Numpy, entre outras. Geralmente, isso é feito no Ubuntu por meio do comando 'sudo apt-get install python-biblioteca', ou 'sudo apt-get install python3- biblioteca', ou algo parecido.
In [1]:
import numpy as np
In [2]:
import matplotlib.pyplot as plt
In [3]:
%matplotlib inline
In [4]:
size = 201 # Tamanho da linha
verbose = False # Define se alguns prints devem ser ativados
if size>10 and verbose:
print('Cuidado! Algumas saídas podem ficar ilegíveis para esse tamanho de malha.')
print('Sugiro definir verbose=False ou diminuir a malha.')
De acordo com o primeiro postulado, o estado do caminhante quântico deve ser um vetor unitário num espaço de Hilbert. Em nosso exemplo, consideramos que o caminhante ocupa uma posição no grafo linha com número de vértices (nós) igual a size. Portanto, o espaço de Hilbert associado a posição é $\mathcal{H}_\mbox{P}$ de dimensão igual a size. A base de $\mathcal{H}_\mbox{size}$ pode ser denotada $\{ | x \rangle : 0 \leq x < \mbox{size} \}$.
Para funcionar como esperamos, o caminhante quântico discreto no tempo ainda precisa de um grau de liberdade adicional, além da posição. Esse grau de liberdade nós chamamos de moeda, pois desempenha um papel semelhante ao lançamento de uma moeda quando comparado ao random walk clássico. Para entender melhor a necessidade de termos essa "moeda quântica", consulte minha tese de doutorado. Por enquanto, basta consideramos que, de acordo com o primeiro postulado, essa moeda quântica também precisa ser representada por um vetor unitário num espaço de Hilbert. O espaço de Hilbert associado a moeda é $\mathcal{H}_\mbox{M}$ de dimensão igual a 2. (Note que só temos dois resultados possíveis para o "lançamento da moeda" nesse caso: o caminhante pode ir para a esquerda ou para a direita.) A base de $\mathcal{H}_\mbox{M}$ pode ser denotada $\{ | m \rangle : 0 \leq m \leq 1 \}$.
De acordo com o terceiro postulado da mecânica quântica, o estado do caminhante é um vetor unitário no espaço $\mathcal{H}_\mbox{M} \otimes \mathcal{H}_\mbox{P}$. Ou seja, é um vetor que num instante de tempo $t$ qualquer pode ser escrito como $$| \Psi (t) \rangle = \sum_{m=0}^{1} \sum_{x=0}^{\mbox{size}-1} \psi_{m;x}(t) | m \rangle \otimes | x \rangle,$$ com $\psi_{m;x} \in \mathbb{C}$. Para que esse vetor seja unitário, precisamos que $\sum_{m}\sum_{x} |\psi_{m;x}(t)|^2 = 1.$
Vejamos como fica isso no Python! Nosso estado inicial será $$| \Psi (0) \rangle = \left(\dfrac{ |0\rangle + i|1\rangle }{\sqrt{2}} \right) \otimes | \frac{\mbox{size}}{2} \rangle,$$ ou seja, a moeda é uma superposição dos dois estados possíveis; e a posição é o "meio" da linha. Mais tarde, você pode trocar essa condição inicial e re-executar o notebook para ver como ela influencia o resultado final.
In [5]:
coinState = np.array([1, 1j], dtype=complex)/np.sqrt(2)
print(coinState)
In [6]:
positionState = np.zeros(size, dtype=complex)
positionState[int(size/2)] = 1.0
""" Para ativar os prints abaixo, a varíavel 'verbose' deve ter sido definida como True. Isso só deve ser feito
juntamente com um valor pequeno para a variável size; digamos, size=10. Caso contrário, a saída ficará ilegível.
Recomendo que você execute essa seção do notebook pelo menos uma vez com os prints ativados, para ajudar no
entendimento. Alguns gráficos a partir da próxima seção, entretanto, ficarão legíveis somente com valores altos
em size. """
if verbose:
print(positionState)
In [7]:
Psi = np.kron(coinState, positionState)
if verbose:
print(Psi)
Vamos confirmar que a norma de $| \Psi \rangle$ é igual a 1. Ficamos satisfeitos se for um número muito próximo de 1, levando em consideração que pode haver erros de arredondamento típicos dos cálculos numéricos.
In [8]:
np.linalg.norm(Psi)
Out[8]:
Para saber a distribuição de probabilidades na posição do caminhante, precisamos aplicar o quarto postulado da mecânica quântica. Você pode encontrar detalhes na minha tese de doutorado. Por enquanto, podemos resumir no seguinte. A probabilidade de, no instante $t$, realizar uma medição e encontrar o caminhante na posição $x$, é dada por $$p(x, t) = \sum_{m=0}^{1} |\psi_{m;x}(t)|^2.$$
Como vamos implementar isso no Python? Precisamos fazer um somatório fixando a posição e variando a moeda. No entanto, depois que aplicamos o produto de Kronecker, nosso vetor de estado ficou aparentemente misturado. Ou seja, só tem um índice para percorrer. Complicado? Nem tanto. Se analisarmos com cuidado como funciona o produto de Kronecker, veremos que podemos identificar precisamente quais entradas do vetor correspondem ao estado da moeda $m=0$, e quais correspondem ao estado da moeda $m=1$. Fica mais fácil ainda se reorganizarmos o nosso vetor, usando o comando reshape do Numpy. Se você tiver definido a variável verbose como True, poderá ver o comando reshape em funcionamento logo abaixo:
In [9]:
if verbose:
print('Assim é o vetor: \n', Psi)
print('\nReorganizando as entradas (sem alterar os valores) ele vai ficar assim: ')
Psi = Psi.reshape((2,size))
print(Psi)
print('Esse formado facilita o cálculo das probabilidades na malha 1D.')
print('\nA qualquer momento posso voltar ao que tínhamos antes. Veja: ')
Psi = Psi.reshape((2*size,))
print(Psi)
Agora sim, usando essa ideia fica fácil calcular as probabilidades no Python! Vejamos:
In [10]:
def get_prob(stateVector, position):
stateVector = stateVector.reshape((2,size))
prob = np.abs(stateVector[0, position])**2 + np.abs(stateVector[1, position])**2
stateVector = stateVector.reshape((2*size,))
return prob
O caminhante inicialmente está localizado no meio da linha, pois assim o definimos. Portanto, a probabilidade de encontrá-lo nessa posição deve ser 1. Vejamos:
In [11]:
get_prob(Psi, size/2)
Out[11]:
E a probabilidade de encontrá-lo em qualquer outra posição deve ser 0. Vejamos:
In [12]:
get_prob(Psi, 0)
Out[12]:
Vamos agora plotar o gráfico da distribuição de probabilidades.
In [13]:
def get_pdf(stateVector):
""" Returns the probability distribution function (pdf) of
a given state vector. """
stateVector = stateVector.reshape((2,size))
pdf = np.sum(np.abs(stateVector)**2, axis=0)
stateVector = stateVector.reshape((2*size,))
return pdf
In [14]:
Prob = get_pdf(Psi)
In [15]:
plt.xlabel('posição')
plt.ylabel('probabilidade')
plt.grid(True)
plt.plot(range(size), Prob)
Out[15]:
Agora que definimos o estado do caminhante, podemos finalmente colocá-lo para caminhar. Para isso, devemos continuar respeitando os postulados da mecânica quântica. Veremos isso na próxima seção.
A evolução da caminhada quântica se dá por um operador unitário da forma $U = S(C \otimes I)$. A matriz $C$ é o operador moeda, e atua somente no subespaço $\mathcal{H}_M$. Em princípio, $C$ pode ser qualquer operador unitário de dimensão compatível com $\mathcal{H}_M$. Nesse exemplo, vamos tomar $C$ como a matriz de Hadamard. Veja a definição na minha tese de doutorado. A matriz $I$ é a identidade que atua no espaço $\mathcal{H}_P$. A matriz $S$ é o operador deslocamento, que atua no espaço $\mathcal{H}_M \otimes \mathcal{H}_P$.
Vamos definir esses operadores, um de cada vez.
In [16]:
hadamardCoin = np.array([[1.0, 1.0],
[1.0, -1.0]]) / np.sqrt(2)
print(hadamardCoin)
A matriz anterior atua somente no espaço moeda. Para termos o operador moeda que atua no espaço composto, ainda precisamos fazer o produto tensorial com o operador identidade que atua no espaço posição. Faremos isso a seguir:
In [17]:
Coin = np.kron(hadamardCoin, np.eye(size))
In [18]:
# checando se é unitária
np.testing.assert_allclose( np.dot(Coin, Coin.T.conj()), np.eye(2*size) )
Agora vamos definir o operador deslocamento. Esse operador modifica a posição de forma dependente da moeda. A estrutura dele, portanto, é um pouco mais complicada. Para entender como ele funciona, vamos por partes. Primeiro: quando a moeda é, digamos, $|0\rangle$, queremos que o caminhante se mova para a direita; e quando a moeda é $|1\rangle$, queremos que o caminhante se mova para a esquerda. Portanto, ao aplicar o operador $S$ queremos o seguinte efeito:
$$S|m\rangle |x\rangle = |m\rangle |x + (-1)^m\rangle,$$
para $m\in\{0,1\}$, e para $0< x < \mbox{size-1}$. Note que aqui usei desigualdade estrita, pois os contornos devem ser tratado como casos à parte -- afinal, quando o caminhante já está em $|0\rangle$, ele não consegue se mover para a esquerda; e quando já está em $|1\rangle$, ele não consegue se mover para a direita. Os casos especiais podem ser definidos de diversas formas, dependendo se a condição de contorno é reflexiva, periódica, etc. Por exemplo, por simplicidade, vamos considerar que seja periódica. Assim, temos uma notação bastante compacta para a matriz associada ao operador $S$, a saber,
$$S = \sum_{m=0}^{1}\sum_{x=0}^{\mbox{size}-1} |m\rangle \langle m| \otimes | x+(-1)^m \rangle \langle x |,$$
em que a operação de soma deve ser realizada modulo size.
In [19]:
def ket(index, dimension=2):
vec = np.zeros((dimension,1), dtype=complex)
vec[index,0] = 1.0
return vec
In [20]:
def bra(index, dimension=2):
vec = np.zeros((1,dimension), dtype=complex)
vec[0,index] = 1.0
return vec
In [21]:
Shift = np.zeros((2*size, 2*size), dtype=complex)
for x in range(0,size):
for m in [0,1]:
Shift += np.kron( np.dot(ket(m, 2),bra(m, 2)),
np.dot( ket((x+(-1)**m) % size, size), bra(x, size)))
In [22]:
# checando se é unitária
np.testing.assert_allclose( np.dot(Shift, Shift.T.conj()), np.eye(2*size, dtype=complex) )
Agora sim, temos tudo que precisamos para definir o operador de evolução $U$.
In [23]:
U = np.dot(Coin, Shift)
In [24]:
# checando se é unitária
np.testing.assert_allclose( np.dot(U, U.T.conj()), np.eye(2*size) )
Para simular a caminhada quântica na reta, defina o número de passos por meio da variável tempo abaixo. Se quiser, você pode voltar ao comando em que definimos o estado inicial do caminhante, e pode modificá-lo para testar diferentes cenários. Dado que o caminhante inicialmente encontrava-se no estado $|\Psi (0) \rangle$, depois de $t$ passos da caminhada ele estará no estado $|\Psi (t)\rangle = U^t |0\rangle.$
In [25]:
tempo = 100
In [26]:
Psit = [Psi]
for t in range(tempo):
newPsi = np.dot(U, Psit[-1])
np.testing.assert_allclose(np.linalg.norm(Psi), 1.0)
Psit.append( newPsi )
A solução acima não é a mais eficiente, principalmente em relação ao consumo de memória. Porém, será bastante útil para fins didáticos. Você agora pode plotar o gráfico da função de distribuição de probabilidades para qualquer instante de tempo simulado. Primeiro, veja como fica para $t=100$. Depois modifique esse valor. Observe como a partícula vai se "espalhando" pela reta.
In [27]:
Prob = get_pdf(Psit[100])
In [28]:
plt.xlabel('posição')
plt.ylabel('probabilidade')
plt.title('Função de Distribuição de Probabilidades')
plt.grid(True)
plt.plot(range(size), Prob)
Out[28]:
Espero que tenha sido útil até aqui. Você pode agora clicar em Kernel / Restart, voltar ao início do notebook e repetir os comandos, fazendo modificações de acordo com sua curiosidade.
Além da função de distribuição de probabilidades, há muitas coisas que podemos querer calcular em nossa simulação de caminhadas quânticas. Por exemplo, nas caminhadas aleatórias clássicas, normalmente queremos conhecer a distribuição estacionária do caminhante, para então calcular o tempo que leva para a caminhada convergir para essa distribuição, etc.
Já mostramos como calcular a distribuição de probabilidades da posição do caminhante quântico no instante $t$. Uma primeira tentativa de definir a distribuição estacionária, portanto, poderia ser algo como $$\bar{P}(x) = \lim_{t\rightarrow \infty} p(x,t).$$ No entanto, essa definição não funciona! Observe que, como $U$ é um operador unitário, $p(x,t)$ não converge com $t\rightarrow \infty$.
Se definirmos a função $$\bar{P}(x) = \lim_{T\rightarrow \infty} \frac{1}{T} \sum_{t=0}^{T} p(x,t),$$ podemos provar que ela converge. Portanto, podemos usar essa função como uma boa definição para a distribuição estacionária de um caminhante quântico.
Vamos criar uma função em Python para calcular $\bar{P}(x)$.
In [29]:
def get_stationary(Psit):
stac = np.zeros(size)
for t in range(len(Psit)):
stac += get_pdf(Psit[t])
stac /= len(Psit)
return stac
Agora, devemos continuar simulando a caminhada quântica por tempos suficientemente longos. Quando maior o valor usado em maisTempo, mais lenta será a simulação, porém melhor será a aproximação de $\bar{P}(x)$.
In [30]:
maisTempo = 10000
for t in range(maisTempo):
newPsi = np.dot(U, Psit[-1])
np.testing.assert_allclose(np.linalg.norm(Psi), 1.0)
Psit.append( newPsi )
In [31]:
print(len(Psit))
Podemos ver agora, com uma boa aproximação, como fica a distribuição estacionária da caminhada quântica na reta.
In [32]:
Stac = get_stationary(Psit)
In [33]:
plt.xlabel('posição')
plt.ylabel('probabilidade')
plt.title('Distribuição Estacionária')
plt.grid(True)
plt.plot(range(size), Stac)
Out[33]:
In [ ]: