UNIVERSIDADE FEDERAL DE SERGIPE Núcleo de Pós-Graduação … · Núcleo de Pós-Graduação em Física Método Kernel Polinomial aplicado a uma rede de spins em ambiente correlacionado

UNIVERSIDADE FEDERAL DE SERGIPENúcleo de Pós-Graduação em Física

DISSERTAÇÃO DE MESTRADO

Método Kernel Polinomial aplicado a umarede de spins em ambiente correlacionado

porGuilherme Martins Alves de Almeida

Cidade Universitária Prof. José Aloísio de CamposSão Cristóvão - SE - Brasil

2012

UNIVERSIDADE FEDERAL DE SERGIPENúcleo de Pós-Graduação em Física

Método Kernel Polinomial aplicado a umarede de spins em ambiente correlacionado

porGuilherme Martins Alves de Almeida

Dissertação apresentada ao Núcleo de Pós-Graduação em Física da Universidade Fe-deral de Sergipe como requisito parcialpara a obtenção do título de Mestre emFísica.Orientador: André Maurício Conceição deSouza

São Cristóvão2012

"I learned very early the difference between knowing the name of something andknowing something."

Richard Feynman

Agradecimentos

À minha família;

À minha namorada Patricia;

Ao meu orientador e grande amigo André Maurício;

Ao Professor Eduardo Mucciolo da University of Central Florida (UCF) por tornar

possível a minha visita aos Estados Unidos, pelo grande aprendizado e pelo apoio

que recebi após a volta ao Brasil. Sou muito grato por essa experiência ímpar;

Aos meus amigos Thiago Rocha do IFT-UNESP, Pejman Jouzdani e Javier Romero

da UCF;

À Sociedade Brasileira de Física;

À CAPES e ao CNPq pelo apoio financeiro;

Muito obrigado.

Resumo

Os bits quânticos, ou qubits, são altamente sensíveis a interações com o ambiente.

O estudo de protocolos visando proteger a informação quântica da descoerência é

essencial para a implementação da computação quântica em larga escala. Boa parte

dos modelos propostos para esta finalidade assume as correlações no ambiente como

inexistentes. Estas podem induzir uma dependência temporal na probabilidade de

erro, comprometendo efetivamente a confiabilidade da informação quântica ao longo

do tempo, mesmo na presença de um código de correção. Sendo assim, devemos le-

var em consideração possíveis limitações físicas na computação quântica tolerante a

falhas. Neste trabalho aplicamos o Método Kernel Polinomial (KPM) no cálculo da

densidade de estados e do decaimento da fidelidade para o código tórico 𝐿 = 3 sem

considerar a dinâmica entre os spins da rede. O modelo Hamiltoniano utilizado con-

siste em um ambiente bosônico livre e um acoplamento spin-bóson, com dois canais

de descoerência, 𝑋 e 𝑍. Uma interação efetiva de longo alcance, anisotrópica, entre

todos os pares de spins da rede é então proposta como um modelo correlacionado.

A correlação está diretamente associada à amplitude e ao alcance da interação entre

os spins. Mostramos que a escala de tempo do decaimento da fidelidade depende

destes fatores.

Palavras-chave: ambientes correlacionados, código tórico, método KPM

Abstract

Quantum bits, or qubits, are highly fragile due to interactions with the environment.

The search for good protocols for protecting quantum information from decoherence

is mandatory in order to make large-scale quantum computation possible. Most of

the models proposed for this assume that correlations in the environment do not

exist. Correlations can induce a time dependent error probability thus seriously

damaging the quantum information over the time even if a quantum correction

code is avaliable. In this way, we must taking into consideration possible physical

limitations to fault-tolerant quantum computing. In this work we apply the Kernel

Polynomial Method (KPM) to evaluate the density of states and fidelity decay of

a 𝐿 = 3 toric code without taking the lattice spin dynamics into account. The

Hamiltonian model is based in a free bosonic environment and a spin-boson coupling,

with two decoherence channels 𝑋 and 𝑍. A long-range, anisotropic interaction

between spin pairs is then proposed as a correlated model. This correlation is directly

related to the interaction strengh and range between spins. We show that the fidelity

decay time scale depends on these parameters.

Keywords: correlated environments, toric code, KPM method

Lista de Figuras

1.1 Rede de spins em um código tórico 𝐿 = 3. . . . . . . . . . . . . . . . 4

1.2 Sítio (s) e plaqueta (p) no código tórico. . . . . . . . . . . . . . . . . 5

1.3 Corrente de operações 𝑍 nos spins da rede. . . . . . . . . . . . . . . 5

3.1 Alguns polinômios de Chebyshev de primeira ordem. . . . . . . . . . 18

4.1 DOS para diferentes valores de 𝑁𝑇 . . . . . . . . . . . . . . . . . . . . 29

4.2 Curva do DOS comparando os kérneis de Dirichlet, Jackson e Lorentz. 30

4.3 DOS para os casos onde só há interações em 𝑋 ou em 𝑍. . . . . . . 31

4.4 Estimando o DOS de acordo com Δ𝐸. . . . . . . . . . . . . . . . . . 31

4.5 Divergência dos coeficientes de expansão 𝜇𝑛. . . . . . . . . . . . . . . 32

4.6 DOS considerando ambos os canais de interação para diferentes valo-

res de 𝛽 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

4.7 DOS para diferentes valores de 𝐽 . . . . . . . . . . . . . . . . . . . . . 33

4.8 Decaimento da fidelidade do estado inicial |0⟩ para diferentes valores

de 𝛽. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35

4.9 Decaimento da fidelidade do estado inicial (|0⟩ + |262143⟩)/√

2 para

diferentes valores de 𝛽. . . . . . . . . . . . . . . . . . . . . . . . . . . 36

4.10 Decaimento da fidelidade do estado inicial (|0⟩+|1⟩+...+|262143⟩)/√𝐷

para diferentes valores de 𝛽. . . . . . . . . . . . . . . . . . . . . . . . 36

4.11 Decaimento da fidelidade do estado inicial (|0⟩ + |262143⟩)/√

2 para

diferentes valores de 𝐽 . . . . . . . . . . . . . . . . . . . . . . . . . . . 37

4.12 Comparação entre o kérneis de Dirichlet, Jackson e Lorentz no decai-

mento da fidelidade. . . . . . . . . . . . . . . . . . . . . . . . . . . . 38

Lista de Tabelas

4.1 Parâmetros utilizados nas análises da densidade de estados e/ou da

fidelidade quântica. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

Sumário

1 Introdução 1

1.1 O código tórico . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

1.2 Objetivos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

2 O modelo de interação 8

2.1 Segunda quantização e o oscilador harmônico . . . . . . . . . . . . . 8

2.2 Acoplamento spin-bóson . . . . . . . . . . . . . . . . . . . . . . . . . 10

2.2.1 Função de correlação: caso bidimensional . . . . . . . . . . . 12

2.2.2 Função de correlação: caso tridimensional . . . . . . . . . . . 13

2.2.3 Função espectral . . . . . . . . . . . . . . . . . . . . . . . . . 13

2.2.4 A interação efetiva entre os spins . . . . . . . . . . . . . . . . 15

3 O Método Kernel Polinomial 16

3.1 Polinômios de Chebyshev . . . . . . . . . . . . . . . . . . . . . . . . 17

3.2 O método KPM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

3.3 Aplicações do método KPM . . . . . . . . . . . . . . . . . . . . . . . 20

3.3.1 Densidade de estados . . . . . . . . . . . . . . . . . . . . . . 20

3.3.2 Evolução temporal . . . . . . . . . . . . . . . . . . . . . . . . 23

4 Resultados e discussão 27

4.1 Cálculo da densidade de estados . . . . . . . . . . . . . . . . . . . . 28

4.2 Análise da fidelidade quântica . . . . . . . . . . . . . . . . . . . . . . 34

5 Conclusões 39

Referências Bibliográficas 40

A Programas em FORTRAN 90 44

A.1 Programas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44

A.2 Sub-rotinas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44

A.3 Códigos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46

A.3.1 mu_coeff . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46

A.3.2 dos_cheb . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50

A.3.3 fidelity_cheb . . . . . . . . . . . . . . . . . . . . . . . . . . 52

A.3.4 TnH . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56

Capítulo 1

Introdução

O estudo de sistemas quânticos visando o processamento de informações é um dos

mais importantes campos de pesquisa da atualidade. O computador quântico pode-

ria realizar certos protocolos significativamente mais rápido do que qualquer meio

clássico, tal como a fatoração de números inteiros através do algoritmo quântico

de Shor [1]. Este pode gerar sérias consequências aos campos da criptografia e da

comunicação1. Naturalmente, a simulação de sistemas quânticos de forma eficiente

é outra grande motivação para a implementação da computação quântica [2].

A unidade de informação de um computador quântico é o chamado bit quântico,

ou qubit, uma extensão do conceito de bit clássico. Os qubits podem assumir os

valores 0 ou 1, como também qualquer combinação linear entre eles [3],

|𝜓⟩ = 𝛼 |0⟩ + 𝛽 |1⟩ , (1.1)

já que estados quânticos podem ser sobrepostos. Esta superposição é a fonte de

todo o potencial da computação quântica. Em geral, um sistema com 𝑛 qubits pode

estar em 2𝑛 diferentes estados simultaneamente. Dessa forma, ao aplicar uma porta

lógica a um estado quântico, pode-se gerar um grande paralelismo computacional.

O qubit pode ser representado por sistemas quânticos de dois níveis bem definidos

1A fatoração de números inteiros grandes é computacionalmente inviável e nisso baseia-se acriptografia de chave pública RSA. Enquanto um algoritmo clássico levaria algo em torno de 104

anos para fatorar um número de 1024 bits, por exemplo, um computador quântico executando oalgoritmo de Shor levaria poucos minutos.

Capítulo 1. Introdução 2

tais como partículas de spin 1/2, o fóton2, entre outros. Contudo, a informação

quântica é bastante frágil. Processos dinâmicos envolvendo um grande número de

graus de liberdade, nos quais qubits interagem com o ambiente, podem resultar

em descoerência [4], fenômeno que gera o colapso da função de onda (o sistema

passa a adquirir comportamento clássico). Este é um dos principais obstáculos

na implementação da computação quântica em larga escala. Sistemas quânticos

reais sofrem interações com o meio externo [5] e para realizar uma computação

quântica confiável, devemos possuir controle sobre diversas formas de ruído. É nesse

sentido que a teoria da correção quântica de erro (QEC) [6] se destaca entre as

diversas áreas da informação quântica, buscando formas de proteger a informação

durante armazenamento, computação e transmissão. Os métodos de QEC consistem,

basicamente, em codificar a informação fazendo uso da redundância e então, após

passar por um canal ruidoso, ela pode ser decodificada recuperando a informação [3].

O código de três qubits para inversão de bit3, por exemplo, consiste em representar

um estado como o da Eq. 1.1 como

|𝜓⟩ = 𝛼 |000⟩ + 𝛽 |111⟩ ≡ 𝛼 |0⟩𝐿 + 𝛽 |1⟩𝐿 , (1.2)

após a inclusão de dois qubits auxiliares e aplicar, duas vezes, a porta CNOT4. O

termo | ⟩𝐿 se refere ao qubit lógico, ou codificado. Dessa forma podemos proteger o

estado contra inversões em apenas um qubit representando a síndrome de erro como

⟨𝜓|𝑃𝑖 |𝜓⟩, 𝑖 = 0,1,2,3, onde 𝑃0 = |000⟩ ⟨000| + |111⟩ ⟨111| é o operador projeção no

estado onde nenhuma inversão ocorreu, 𝑃1 = |100⟩ ⟨100|+ |011⟩ ⟨011| no estado onde

o 1º qubit foi invertido e assim por diante. Uma dada combinação dos valores da

síndrome de erro indicam qual procedimento deve ser realizado para a recuperação

do estado, isto é, informa qual qubit foi invertido.

Nas técnicas de correção quântica de erro, um qubit lógico é codificado em um

conjunto de qubits físicos, tal que os erros mais prováveis podem ser detectados e

2O fóton pode ser polarizado horizontal e verticalmente.3Inversão de bit: 𝜎𝑥 = 𝑋 =

(︂0 11 0

)︂; inversão de fase: 𝜎𝑧 = 𝑍 =

(︂1 00 −1

)︂.

4A porta quântica CNOT (não-controlado), inverte o segundo qubit (alvo) se o primeiro (con-trole) for 1


corrigidos periodicamente. Aplicando-se recursivamente estas codificações, proce-

dimento denominado concatenação, é possível atenuar os erros arbitrariamente ou,

equivalentemente, proteger a informação quântica de maneira robusta [3]. É impor-

tante destacar que as portas lógicas utilizadas para codificar e decodificar também

podem ser ruidosas. Entretanto, o teorema do limiar [7] estabelece que se o ruído

em um computador quântico estiver abaixo de um certo valor limite, digamos, se a

probabilidade de erro 𝑝 em qualquer componente deste for menor do que um 𝑝𝑙𝑖𝑚𝑖𝑎𝑟,

é possível implementar computação quântica em larga escala de forma eficiente.

Este teorema, no entanto, foi satisfatoriamente provado para modelos de erros es-

tocásticos ideais, sem considerar a origem física do fenômeno. Erros que surgem

de processos dinâmicos na interação de qubits com o ambiente possuem uma escala

de tempo característica, assim como o potencial que induz correlações temporais

e espaciais na evolução dos qubits no tempo. A diferença chave entre um modelo

de erro correlacionado e um estocástico é que no primeiro a probabilidade de erro

possui dependência temporal e/ou local [8]. Quando a taxa de erro não é constante,

torna-se mais difícil demonstrar se existe um limiar.

Na estimativa do valor do limiar em [7], o ruído foi assumido como Markoviano,

onde cada porta quântica no circuito é um mapeamento completamente positivo

que preserva o traço se aproximando então de uma porta quântica ideal. Também

foram assumidos modelos de erros nos quais as correlações espaciais e temporais

decaem exponencialmente. Estes modelos, portanto, não se relacionam diretamente

com modelos físicos descrevendo a descoerência de maneira mais detalhada. Neste

sentido, muito têm sido investigado a respeito de modelos realísticos de erros em

sistemas de qubits e suas implicações na computação quântica tolerante a falhas

[8, 9, 10, 11, 12].

1.1 O código tórico

Um estudo promissor se refere a como ambientes correlacionados afetam as chama-

das memórias quânticas que se auto-corrigem [13]. Nestes sistemas, a geometria e as


interações são dispostas de modo que um gap não-trivial separa o estado fundamen-

tal (que engloba a base lógica) dos estados excitados. Esta proteção é topológica,

portanto, leva um tempo exponencialmente grande ou uma perturbação altamente

não-local para induzir um erro lógico. O exemplo mais conhecido de memória quân-

tica topológica é o código tórico introduzido por Kitaev [14]. Neste código, dispomos

de partículas de spin 1/2 situadas nas arestas de uma rede bidimensional 𝐿x𝐿 com

periodicidade, formando um toróide. Para uma rede de tamanho 𝐿, há 2𝐿2 spins

(Fig. 1.1). Os operadores responsáveis pela síndrome de erro atuam em quatro spins

contidos em um sitio ou plaqueta na rede definida no toróide (Fig. 1.2). São eles os

operadores sítio,

𝐴𝑠 =∏︁𝑗∈𝑠

𝑋𝑗 , (1.3)

e plaqueta,

𝐵𝑝 =∏︁𝑗∈𝑝

𝑍𝑗 . (1.4)

Estes possuem autovalores 1 e −1. O subespaço protegido do código5 é formado

por todos os estados |𝜓⟩ de 2𝐿2 spins que satisfazem 𝐴𝑠 |𝜓⟩ = |𝜓⟩ e 𝐵𝑝 |𝜓⟩ = |𝜓⟩

para todo sítio e plaqueta. É fácil observar que 2𝐿2 − 2 geradores do código são

Figura 1.1: Rede de spins em um código tórico 𝐿 = 3. Os spins estão situados nas arestasdas plaquetas. Os marcadores sem preenchimento correspondem aos spins do respectivolado oposto (a rede é periódica).

5O código tórico é um código estabilizador, e por subespaço protegido nos referimos ao conjuntodos estados invariantes sob aplicação dos operadores que definem o código, denominados geradoresdo código. Portanto, este espaço é uma interseção dos subespaços definidos por cada gerador comautovalores iguais a 1. Para mais detalhes consultar a Ref. [3], Cap. 10.


p

s

Figura 1.2: Sítio (s) e plaqueta (p) no código tórico. Os operadores correspondentes, 𝐴𝑠 e𝐵𝑝, envolvem a interação entre quatro spins da rede. Um código tórico de tamanho 𝐿 possui𝐿2 sítios e plaquetas.

independentes. Desta maneira, este subespaço possui quatro dimensões, ou seja, o

código tórico protege dois qubits lógicos [13, 14].

Erros são causados por correntes de inversão de bit (𝑋) ou de fase (𝑍). Assu-

mindo que estas correntes não circulem ao redor do toróide (Fig. 1.3), nenhuma

operação lógica indesejada pode ocorrer, já que a corrente é um produto dos opera-

dores 𝐴𝑠 ou 𝐵𝑝 (o produto também está contido no subespaço protegido).

O código tórico nos permite realizar correções de erro a nível físico. Definimos a

Hamiltoniana [14]

𝐻tor = −∑︁

𝑠

𝐴𝑠 −∑︁

𝑝

𝐵𝑝, (1.5)

onde o estado fundamental é quatro vezes degenerado e os estados são os mesmos

C

C*

Figura 1.3: Corrente de operações 𝑍 nos spins da rede. O ciclo C é equivalente a um produtode operadores plaqueta e portanto não altera o estado codificado. Quando o ciclo contornao toróide (C*), o subespaço dos estados protegidos (o espaço estabilizador) é preservadomas ocorre uma operação lógica [13]. Para correntes de inversão bit, valem as mesmaspropriedades. Isto se torna evidente quando representamos o código pela rede dual, onde osoperadores sítio passam a ser plaqueta, mas envolvendo interações 𝑋 entre quatro spins.


dos contidos no subespaço protegido do código. Ao acoplar esta Hamiltoniana a um

banho térmico em baixa temperatura, erros podem ser removidos automaticamente

por processos de dissipação.

Neste trabalho não vamos incluir a proteção topológica. Não haverá, portanto,

dinâmica entre os spins da rede. Dito isso, é natural pensar sobre qual seria a

relevância em estudar o código tórico sem considerar a auto-correção. A vantagem

deste código está não só na maneira de codificar a informação, mas também no modo

como se extrai essa informação. Ao invés de ler spins individuais, as síndromes de

erro são feitas pelos operadores plaquetas e sítios, ou seja, em conjunto de spins.

Isso possibilita a correção de erros em ordens muito altas. Uma desvantagem é que

usa-se um grande numero de spins para codificar apenas dois qubits lógicos.

A geometria toroidal não é a única existente para formar este tipo de código

quântico de correção de erro. Na prática, ela se torna inconveniente se nos depa-

rarmos com situações onde é necessário que qubits residindo em diferentes toróides

interajam no curso da computação quântica [13]. Outras superfícies podem ser uti-

lizadas, com propriedades topológicas que determinam a degenerescência do estado

fundamental. Códigos de correção definidos em redes bidimensionais de spins são

conhecidos como códigos de superfície.

Em um protótipo experimental com doze qubits físicos feitos de junções de Jo-

sephson em nanoescala, foi observado que o qubit lógico é protegido contra variações

de fluxo magnético [15]. Também foi demonstrado experimentalmente, utilizando

um estado de oito fótons, que correlações permanecem invioladas ao ocorrer um

único erro em qualquer qubit [16]. Estes resultados sugerem que códigos topológicos

de correção quântica de erro são viáveis.

1.2 Objetivos

Considerando um código tórico 𝐿 = 3 em um ambiente bosônico, vamos calcular

a densidade de estados (DOS) e a evolução temporal do sistema para análise da

fidelidade quântica. Esta é uma das grandezas que consegue qualificar a informação


quântica após processos envolvendo descoerência e/ou erros unitários. O modelo

Hamiltoniano consistirá num conjunto infinito de osciladores harmônicos, ou banho

bosônico, e uma interação spin-bóson. A interação induz uma interação efetiva de

longo alcance entre todos os pares de spins. Há um parâmetro associado ao alcance e

outro associado à intensidade da interação. Estes parâmetros estão diretamente liga-

dos à correlação do sistema. Partindo desse modelo, vamos utilizar o Método Kernel

Polinomial (KPM), baseado na expansão de funções por polinômios de Chebyshev,

para calcular as grandezas mencionadas acima.

Esta dissertação esta organizada da seguinte forma: no Cap. 2 introduzimos o

modelo Hamiltoniano; no Cap. 3 discutimos aspectos gerais sobre do método KPM,

suficientes para a finalidade deste trabalho; no Cap. 4 apresentamos os resultados

obtidos para o DOS e para o decaimento da fidelidade no tempo; no Cap. 5 fazemos

as devidas considerações finais acerca dos resultados. O Apêndice A disponibiliza

alguns dos programas e sub-rotinas em FORTRAN 90 utilizados para obtenção dos

dados numéricos.

Capítulo 2

O modelo de interação

Neste capítulo vamos descrever a dinâmica de spins 1/2 sujeitos a interações com os

graus de liberdade externos. Uma maneira frequentemente empregada para descre-

ver o ambiente é considerá-lo como um um conjunto infinito de osciladores harmôni-

cos, ou um campo bosônico [17]. Nestes sistemas quânticos abertos, a Hamiltoniana

pode ser expressa da forma

𝐻(𝑡) = 𝐻𝑆 +𝐻𝐵 + 𝑉 (𝑡), (2.1)

onde 𝐻𝑆 descreve a dinâmica dos spins somente, 𝐻𝐵 o banho bosônico e 𝑉 (𝑡) é a

interação sistema-banho. Em geral, considera-se todo o sistema 𝑆+𝐵 como isolado.

Nas próximas seções, discutiremos o modelo Hamiltoniano levado em consideração

neste trabalho.

2.1 Segunda quantização e o oscilador harmônico

Sistemas envolvendo muitas partículas idênticas são frequentemente descritos pelo

formalismo de segunda quantização [18]. Seja um estado de 𝑁 bósons expresso por1

|𝑛1, 𝑛2, ...⟩ onde

∑︁𝑖

𝑛𝑖 = 𝑁, (2.2)

1Bósons não obedecem o princípio de exclusão de Pauli e portanto podem ocupar um mesmoestado quântico.

Capítulo 2. O modelo de interação 9

definimos o operador número �̂�𝑗 como

�̂�𝑗 |𝑛𝑗⟩ = 𝑛𝑗 |𝑛𝑗⟩ , (2.3)

onde 𝑛𝑗 é o número de bósons no estado 𝑗. Vamos agora introduzir os operadores

de criação e aniquilação,

𝑏†𝑗 |𝑛1, 𝑛2, ..., 𝑛𝑗 , ...⟩ =

√︁𝑛𝑗 + 1 |𝑛1, 𝑛2, ..., 𝑛𝑗 + 1, ...⟩ , (2.4)

𝑏𝑗 |𝑛1, 𝑛2, ..., 𝑛𝑗 , ...⟩ = √𝑛𝑗 |𝑛1, 𝑛2, ..., 𝑛𝑗 − 1, ...⟩ , (2.5)

respectivamente, que aumenta ou diminui o número de ocupação 𝑛𝑗 do estado |𝑛𝑗⟩.

Caso um estado possua 𝑛𝑗 = 0, obteremos 𝑏𝑗 |𝑛1, ..., 𝑛𝑗 , ...⟩ = 0.

O comutador [𝐴,𝐵] entre dois operadores 𝐴 e 𝐵 é definido por

[𝐴,𝐵] = 𝐴𝐵 −𝐵𝐴. (2.6)

Dessa forma, os operadores 𝑏†𝑗 e 𝑏𝑗 satisfazem as seguintes relações de comutação:

[𝑏†𝑗 , 𝑏

†𝑘] = 0; (2.7)

[𝑏𝑗 , 𝑏𝑘] = 0; (2.8)

[𝑏𝑗 , 𝑏†𝑘] = 𝛿𝑗𝑘. (2.9)

𝛿𝑗𝑘 é o Delta de Kronecker. Tendo introduzido as principais relações da segunda

quantização para bósons, podemos agora mostrar como o oscilador harmônico pode

ser formulado por este método. Em uma dimensão, sua Hamiltoniana é dada por

𝐻 = 𝑝2

2𝑚 + 12𝑚𝜔

2𝑥2. (2.10)

Definimos os operadores de criação e destruição como

𝑎 ≡ 𝑚𝜔

2~

(︂𝑥+ 𝑖𝑝

𝑚𝜔

)︂, (2.11)

𝑎† ≡ 𝑚𝜔

2~

(︂𝑥− 𝑖𝑝

𝑚𝜔

)︂, (2.12)

respectivamente. Sendo a relação de comutação entre 𝑥 e 𝑝 dada por

[𝑝, 𝑥] = ~/𝑖 (2.13)


é fácil verificar que 𝑎 e 𝑎† satisfazem as Eqs. (2.7) a (2.9). O operador número é

definido por

𝑁 = 𝑎†𝑎 = 𝐻

~𝜔− 1

2 , (2.14)

e assim obtemos o Hamiltoniano para o oscilador harmônico em segunda quantiza-

ção,

𝐻 = ~𝜔(︂𝑎†𝑎+ 1

2

)︂. (2.15)

Os auto-estados do operador 𝐻 são os mesmos de 𝑁 já que [𝑁,𝐻] = 0 e então eles

podem ser simultaneamente diagonalizáveis2 [19]. Os auto-estados de energia são

escritos como

|𝑛⟩ = (𝑎†)𝑛

√𝑛!

|0⟩ , (2.16)

e portanto

𝐻 |𝑛⟩ = ~𝜔(︂𝑛+ 1

2

)︂|𝑛⟩ = 𝐸𝑛 |𝑛⟩ . (2.17)

Podemos interpretar este resultado como a energia de 𝑛 bósons. Desta forma a

Hamiltoniana que descreve o ambiente bosônico considerado neste trabalho é dada

por (~ = 1)

𝐻𝐵 =∑︁k ̸=0

𝜔k

(︂𝑎†

k𝑎k + 12

)︂, (2.18)

onde 𝜔k = 𝑣|k|, ou seja, estamos considerando bósons acústicos sem massa.

2.2 Acoplamento spin-bóson

Diversos sistemas físicos ruidosos podem ser formalizados pelo modelo spin-bóson

[17]. Este modelo têm sido frequentemente utilizado no contexto de mecânica quân-

tica dissipativa. A interação entre bósons e os spins de uma rede pode ser descrita

2Neste caso, isto é óbvio já que 𝐻 é uma função linear de 𝑁 .


pela Hamiltoniana [12]

𝑉 = 𝜆

2∑︁

𝑗

∑︁k ̸=0

𝑐|k|𝑠(︁𝑒𝑖k·x𝑗𝑎†

k + 𝑒−𝑖k·x𝑗𝑎k)︁𝜎𝑧

𝑗 , (2.19)

onde 𝜆 é a intensidade da interação, 𝑠 um expoente dinâmico e 𝑐 um ajuste dimen-

sional. Sem considerar a dinâmica entre os spins (𝐻𝑆 = 0), vamos reescrever a Eq.

(2.19) na representação da interação como

𝑉 (𝑡) = 𝑒𝑖𝐻𝐵𝑡𝑉 𝑒−𝑖𝐻𝐵𝑡 = 𝜆

2∑︁

𝑗

∑︁k ̸=0

𝑐|k𝑠|(︁𝑒𝑖k·x𝑗𝑒𝑖𝜔k𝑡𝑎†

k + 𝑒−𝑖k·x𝑗𝑒−𝑖𝜔k𝑡𝑎k)︁𝜎𝑧

𝑗 . (2.20)

Definimos o campo bosônico como

𝜃(x, 𝑡) = −𝑖𝑐∑︁k ̸=0

(︁𝑒𝑖k·x𝑗𝑒𝑖𝜔k𝑡𝑎†

k − 𝑒−𝑖k·x𝑗𝑒−𝑖𝜔k𝑡𝑎k)︁ |k|𝑠

𝜔k, (2.21)

de modo que

𝜕𝑡𝜃(x, 𝑡) = 𝑐∑︁k ̸=0

(︁𝑒𝑖k·x𝑗𝑒𝑖𝜔k𝑡𝑎†

k + 𝑒−𝑖k·x𝑗𝑒−𝑖𝜔k𝑡𝑎k)︁

|k|𝑠 (2.22)

e reescrevemos a Eq. (2.23) como

𝑉 (𝑡) = 𝜆

2∑︁

𝑗

𝜕𝑡𝜃(x𝑗 , 𝑡)𝜎𝑧𝑗 . (2.23)

O termo de interação spin-bóson induz uma interação efetiva entre os pares de spins.

Para encontrar sua forma, vamos calcular a função de correlação,

⟨𝑉 (𝑡)𝑉 (𝑡′)⟩ = 𝜆2

4∑︁𝑗,𝑙

⟨𝜕𝑡𝜃(x𝑗 , 𝑡)𝜕𝑡𝜃(x𝑙, 𝑡′)⟩𝜎𝑧

𝑗𝜎𝑧𝑙 . (2.24)

Como os campos bosônicos são gaussianos, a Hamiltoniana efetiva quadrática é a

única que devemos levar em conta já que todos os termos de acoplamento de ordem

alta se anulam. Calculamos então o correlator bosônico como

𝐹 (x𝑗 − x𝑙, 𝑡− 𝑡′) = ⟨𝜕𝑡𝜃(x𝑗 , 𝑡)𝜕𝑡𝜃(x𝑙, 𝑡′)⟩

= 𝑐2 ∑︁k ̸=0

|k|2𝑠𝑒−𝑖k·(x𝑗−x𝑙)𝑒−𝑖𝜔k(𝑡−𝑡′), (2.25)


onde os valores esperados são calculados no vácuo dos campos bosônicos, isto é, nos

quais

⟨𝑎†k𝑎k⟩ = 0, (2.26)

⟨𝑎k𝑎†k⟩ = 0. (2.27)

Reescrevendo a soma da Eq. (2.25) como contínua temos

𝐹 (x, 𝑡) = 𝑐2𝐿𝑑∫︁ d𝑑𝑘

(2𝜋)𝑑|k|2𝑠𝑒−𝑖k·x𝑒−𝑖𝜔k𝑡. (2.28)

A partir daqui devemos especificar a dimensão espacial 𝑑.

2.2.1 Função de correlação: caso bidimensional

Para 𝑑 = 2 temos,

𝐹 (x, 𝑡) = 𝑐2𝐿2∫︁ Λ

0

d𝑘 𝑘2𝑠+1

(2𝜋)2 𝑒−𝑖𝜔𝑘𝑡∫︁ 2𝜋

0d𝜑 𝑒𝑖𝑘𝑥cos𝜑

= 𝑐2𝐿2∫︁ Λ

0

d𝑘 𝑘2𝑠+1

(2𝜋)2 𝑒−𝑖𝜔𝑘𝑡𝐽0(𝑘𝑥), (2.29)

onde Λ é um corte ultravioleta3 e 𝐽0(𝑘𝑥) é a função de Bessel de primeira espécie de

ordem 0. Esta integral pode ser calculada exatamente se considerarmos interação

instantânea, isto é, 𝑡 = 0. Obtemos então,

𝐹 (x, 0) = 𝑐2𝐿2

2𝜋𝑥2𝑠+2

∫︁ Λ𝑥

0d𝑞 𝑞2𝑠+1𝐽0(𝑞). (2.30)

Consultando a Ref. [21], Eq. 6.561.14,

∫︁ ∞

0d𝑥𝑥𝜇𝐽𝜈(𝑎𝑥) = 2𝜇𝑎−𝜇−1

Γ(︂1

2 + 12𝜈 + 1

2𝜇)︂

Γ(︂1

2 + 12𝜈 − 1

2𝜇)︂ (2.31)

e fazendo Λ𝑥 → ∞, obtemos

∫︁ ∞

0d𝑞 𝑞2𝑠+1𝐽0(𝑞) = 22𝑠+1 Γ(1 + 𝑠)

Γ(−𝑠) , (2.32)

3Em Teoria Quântica de Campos, este termo se refere a um parâmetro introduzido para evitardivergências nas funções de correlação. O corte ultravioleta é então um corte nas altas energias.Mais detalhes podem ser encontrados na Ref. [20].


desde que −1 < 𝑠 < −1/4. Para 𝑠 < −1, a integral possui uma divergência infraver-

melho, que leva em consideração o volume do sistema bosônico, e para 𝑠 > −1/4,

o limite Λ𝑥 → ∞ não pode ser estabelecido, e a integral passa a depender do corte

ultravioleta Λ.

2.2.2 Função de correlação: caso tridimensional

Para 𝑑 = 3 temos,

𝐹 (x, 𝑡) = 𝑐2𝐿3∫︁ Λ

0

d𝑘 𝑘2𝑠+2

(2𝜋)3 𝑒−𝑖𝜔𝑘𝑡∫︁ 2𝜋

0d𝜑∫︁ 𝜋

0d𝜃 sin𝜃𝑒𝑖𝑘𝑥cos𝜃

= 2𝑐2𝐿3

𝑥

∫︁ Λ

0

d𝑘 𝑘2𝑠+1

(2𝜋)2 𝑒−𝑖𝜔𝑘𝑡sin(𝑘𝑥). (2.33)

Novamente, considerando 𝑡 = 0, obtemos

𝐹 (x, 0) = 2𝑐2𝐿3

(2𝜋)2𝑥2𝑠+3

∫︁ Λ𝑥

0d𝑞 𝑞2𝑠+1sin(𝑞). (2.34)

Consultando a Ref. [21], Eq. 3.761.4,

∫︁ ∞

0d𝑥𝑥𝜇−1sin(𝑎𝑥) = Γ(𝜇)

𝑎𝜇sin𝜇𝜋2 =

𝜋sec𝜇𝜋22𝑎𝜇Γ(1 − 𝜇) , (2.35)

obtemos, com Λ𝑥 → ∞,

∫︁ ∞

0d𝑞 𝑞2𝑠+1sin(𝑞) = Γ(2𝑠+ 2)sin[(𝑠+ 1)𝜋], (2.36)

para o intervalo −1 < 𝑠 < −1/2. Observe que esta integral possui divergência

infravermelha quando 𝑠 < −1 e uma dependência do corte ultravioleta quando

𝑠 > −1/2.

2.2.3 Função espectral

Para fornecer uma boa descrição do mecanismo de dissipação, é conveniente calcular

a função espectral do banho bosônico. Para um sistema interagente descrito por

𝐻𝑖𝑛𝑡 =∑︁

𝑖

∑︁𝑘

𝑔𝑘𝜎𝑧𝑖 𝐹𝑘(𝑖), (2.37)


e

𝐻𝐵 =∑︁

𝑖

∑︁𝑘

𝜔𝑘𝐹†𝑘 (𝑖)𝐹𝑘(𝑖), (2.38)

a função espectral é definida como

𝜌(𝜔) =∑︁

𝑘

|𝑔𝑘|2𝛿(𝜔 − 𝜔𝑘), (2.39)

onde 𝛿(𝜔−𝜔𝑘) é o delta de Dirac. Para o nosso modelo, portanto, a função espectral

é dada por

𝜌(𝜔) = 𝜆2

4∑︁k ̸=0

𝑐2|k|2𝑠𝛿(𝜔 − 𝜔k). (2.40)

Considerando o limite contínuo, temos

𝜌(𝜔) = 𝜆2

4 𝑐2𝐿𝑑

∫︁ d𝑑𝑘

(2𝜋)𝑑𝛿(𝜔 − 𝑣𝑘)𝑘2𝑠. (2.41)

Se 𝑑 = 2,

𝜌(𝜔) = 𝜆2

8𝜋𝑐2𝐿2

∫︁ Λ

0d𝑘 𝑘2𝑠+1𝛿(𝜔 − 𝑣𝑘)

= 𝜆2

8𝜋𝑐2𝐿2𝜔

2𝑠+1

𝑣2𝑠+2 Θ(Λ𝑣 − 𝜔). (2.42)

Dessa maneira, para baixas frequências, 𝜌(𝜔) ∼ 𝜔2𝑠+1, de acordo com a classificação

padrão [17], podemos distinguir os regimes: a) 𝑠 > 0, super-ôhmico; b) 𝑠 = 0,

ôhmico e c) 𝑠 < 0, sub-ôhmico. Para o regime ôhmico, a função de correlação, Eq.

(2.30), decai como 1/𝑥2. Ela decai mais rapidamente no caso super-ôhmico e mais

devagar no caso sub-ôhmico. Para 𝑑 = 3, a função espectral é dada por

𝜌(𝜔) = 𝜆2

8𝜋2 𝑐2𝐿3

∫︁ Λ

0d𝑘 𝑘2𝑠+2𝛿(𝜔 − 𝑣𝑘)

= 𝜆2

8𝜋2 𝑐2𝐿3𝜔

2𝑠+2

𝑣2𝑠+3 Θ(Λ𝑣 − 𝜔). (2.43)

Para baixas frequências, notamos que 𝜌(𝜔) ∼ 𝜔2𝑠+2. Da mesma maneira, identi-

ficamos os três regimes: a) 𝑠 > −1/2, super-ôhmico; b) 𝑠 = −1/2, ôhmico e c)

𝑠 < −1/2, sub-ôhmico. Assim como no caso bidimensional para o caso ôhmico, a

função de correlação, Eq. (2.34) também decai como 1/𝑥2. O decaimento quadrático


é então a principal característica do regime ôhmico em qualquer dimensão espacial

se os bósons são acústicos e não possuem massa.

2.2.4 A interação efetiva entre os spins

Sendo ⟨𝑉 (𝑡)𝑉 (𝑡′)⟩ o único termo de correlação não nulo e 𝐹 (𝑥) ∼ 1/𝑥2𝑠+𝑑 com

𝐹 (𝑥) ∼ 1/𝑥2 quando o regime é ôhmico, redefinimos 2𝑠 + 𝑑 ≡ 2 + 𝛽 onde 𝛽 torna-

se o novo parâmetro que caracteriza o tipo de dissipação. Dessa forma, podemos

modelar uma interação efetiva entre os pares de spins induzida pelo banho bosônico

a 𝑇 = 0 de acordo com

𝑉𝑒𝑓𝑓 = 12∑︁𝑗 ̸=𝑙

∑︁𝛼=𝑥,𝑦,𝑧

𝐽𝛼

|x𝑗 − x𝑙|2+𝛽𝜎𝛼

𝑗 𝜎𝛼𝑙 , (2.44)

onde os vetores posição são dados em unidades da constante de rede e 𝐽𝛼 > 0 para

qualquer 𝛼. Foi acrescentado o canal de interação 𝑋. A Hamiltoniana descreve um

acoplamento anisotrópico, anti-ferromagnético de longo alcance entre os pares de

spins. O parâmetro 𝛽 está associado ao alcance da interação e identificamos os três

regimes: a) 𝛽 > 0, super-ôhmico; b) 𝛽 = 0, ôhmico e c) 𝛽 < 0, sub-ôhmico. Vamos

assumir essa interação efetiva como instantânea, desconsiderando, portanto, efeitos

de retardação relacionados à velocidade finita dos bósons. Isto é válido no caso de

fraco acoplamento e bósons com altas velocidades.

Capítulo 3

O Método Kernel Polinomial

A análise de sistemas cujas propriedades dependem da ação conjunta de inúme-

ros graus de liberdade ou várias interações atuando numa mesma escala energética,

requer métodos numéricos eficientes e compatíveis com o poder computacional aces-

sível. O cálculo das propriedades de diversos sistemas microscópicos está limitado

aos recursos disponíveis para a diagonalização da Hamiltoniana, ou seja, à obtenção

dos autovalores e autovetores. Este procedimento demanda uma memória compu-

tacional de 𝑂(𝐷2) e um tempo de processamento 𝑂(𝐷3), onde 𝐷 é a dimensão

da matriz. Neste capítulo vamos brevemente descrever o Método Kernel Polino-

mial (KPM) [22, 23] baseado na expansão de funções através dos polinômios de

Chebyshev de primeira ordem. Este método demanda recursos em 𝑂(𝐷) para ma-

trizes esparsas possibilitando, então, a análise de sistemas com espaço de Hilbert de

dimensões maiores [22]. Na Seção 3.1 introduziremos os polinômios de Chebyshev e

em seguida, os fundamentos do método KPM na Seção 3.2. Por fim, na Seção 3.3

mostramos como aplicar o método no cálculo da densidade de estados e da evolução

temporal de um sistema.

Capítulo 3. O Método Kernel Polinomial 17

3.1 Polinômios de Chebyshev

Sejam 𝑓(𝑥) e 𝑔(𝑥) funções contínuas 𝑓, 𝑔 : [𝑎, 𝑏] → R, onde, dado um conjunto

completo de polinômios 𝑝𝑛(𝑥), podem ser expandidas da seguinte maneira:

𝑓(𝑥) =∞∑︁

𝑛=0𝛼𝑛𝑝𝑛(𝑥), (3.1)

𝑔(𝑥) =∞∑︁

𝑚=0𝛽𝑚𝑝𝑚(𝑥), (3.2)

onde {𝛼𝑛} e {𝛽𝑚} são os coeficientes da expansão. Dada uma função peso 𝑊 (𝑥),

𝑊 : [𝑎, 𝑏] → R, podemos escrever o produto interno como

⟨𝑓 |𝑔⟩ =∫︁ 𝑏

𝑎d𝑥𝑊 (𝑥)𝑓(𝑥)𝑔(𝑥). (3.3)

A escolha do conjunto ortogonal {𝑝𝑛} vai depender da finalidade em questão. Para

diversas, uma boa opção são os polinômios de Chebyshev [22, 23, 24, 25, 26, 27],

principalmente devido à rápida convergência de sua expansão, como veremos no

Cap. 4. Vamos utilizar apenas os polinômios de primeira ordem 𝑇𝑛(𝑥). Estes são

soluções da equação diferencial de Chebyshev1 ,

(1 − 𝑥2)d2𝑦

d𝑥2 − 𝑥d𝑦d𝑥 + 𝑝2𝑦 = 0, (3.4)

onde 𝑝 é número real constante. Os polinômios de Chebyshev de primeira ordem

obedecem a uma simples relação de recorrência,

𝑇𝑛+1(𝑥) = 2𝑥𝑇𝑛(𝑥) − 𝑇𝑛−1(𝑥), (3.5)

𝑛 > 0, com 𝑇0(𝑥) = 1 e 𝑇1(𝑥) = 𝑥. A Fig. 3.1 mostra os primeiros polinômios. Eles

são definidos no intervalo [−1, 1], com função peso 𝑊 (𝑥) = (𝜋√

1 − 𝑥2)−1 de modo

que reescrevemos a Eq. (3.3) como

⟨𝑓 |𝑔⟩ =∫︁ 1

−1d𝑥 𝑓(𝑥)𝑔(𝑥)

𝜋√

1 − 𝑥2. (3.6)

A relação de ortogonalidade é então expressa por

1As soluções são obtidas por série de potências 𝑦 =∑︀∞

𝑗=0 𝑎𝑗𝑥𝑗 [28].


-1,0 -0,5 0,0 0,5 1,0

-1,0

-0,5

0,0

0,5

1,0

n = 2

n = 3n = 4

n = 1

Tn(x

)

x

Figura 3.1: Primeiros polinômios de Chebyshev 𝑇𝑛(𝑥) para 𝑥 ∈ [−1, 1].

⟨𝑇𝑛|𝑇𝑚⟩ =∫︁ 1

−1d𝑥 𝑇𝑛(𝑥)𝑇𝑚(𝑥)

𝜋√

1 − 𝑥2. (3.7)

Podemos resolver esta integral fazendo a transformação 𝑥 = cos(𝜃), 𝑇𝑛(cos(𝜃)) =

cos(𝑛𝜃), obtendo

⟨𝑇𝑛|𝑇𝑚⟩ =∫︁ 𝜋

0d𝜃 𝜋−1cos(𝑛𝜃)cos(𝑚𝜃)

= 1 + 𝛿𝑛,02 𝛿𝑛, 𝑚. (3.8)

Os aspectos gerais sobre os polinômios de Chebyshev descritos anteriormente são

suficientes para a finalidade deste trabalho. Mais detalhes podem ser encontrados

na Ref. [28].

3.2 O método KPM

O método KPM aproxima funções definidas em um intervalo fechado por uma soma

polinomial truncada. Seja 𝑓(𝑥) uma função real definida no intervalo [−1, 1], pode-

mos aproximá-la como [24]

𝑓𝐾𝑃 𝑀 (𝑥) =∫︁ 1

−1d𝑦 𝑓(𝑦)𝜋

√︁1 − 𝑦2𝐾𝑁𝑇

(𝑥, 𝑦), (3.9)


onde o kernel2 é definido por, para uma expansão de ordem 𝑁𝑇 ,

𝐾𝑁𝑇(𝑥, 𝑦) = 𝑔0𝜑0(𝑥)𝜑0(𝑦) + 2

𝑁𝑇 −1∑︁𝑛=1

𝑔𝑛𝜑𝑛(𝑥)𝜑𝑛(𝑦) (3.10)

e

𝜑𝑛(𝑥) = 𝑇𝑛(𝑥)𝜋

√1 − 𝑥2

. (3.11)

Diferentes coeficientes {𝑔𝑛} definem um determinado kernel. Entre os mais utiliza-

dos, podemos citar [22]:

1. O kernel de Dirichlet,

𝑔𝑛 = 1. (3.12)

Esta é escolha mais simples. É geralmente utilizado em situações onde há

rápida convergência da expansão de Chebyshev como veremos adiante para o

operador evolução temporal (Sec. 3.3.2). Caso contrário, este fornece pouca

precisão e flutuações relacionadas à descontinuidades em funções (caso haja).

2. O kernel de Jackson,

𝑔𝑛 = 1𝑁𝑇 + 1

[︂(𝑁𝑇 − 𝑛+ 1) cos

(︂𝜋𝑛

𝑁𝑇 + 1

)︂+ sin

(︂𝜋𝑛

𝑁𝑇 + 1

)︂cot

(︂𝜋

𝑁𝑇 + 1

)︂]︂. (3.13)

Este é o mais conveniente para a maioria das aplicações, fornecendo resolução

proporcional a 1/𝑁𝑇 para ordens mais altas.

3. O kernel de Lorentz,

𝑔𝑛 = sinh[𝜆(1 − 𝑛/𝑁𝑇 )]sinh(𝜆) , (3.14)

onde 𝜆 é um parâmetro real, é mais apropriado em situações envolvendo cálculo

das funções de Green [25].

2Na matemática, a palavra kernel possui diversos significados. Nesse contexto nos referimos àfunção kernel de uma transformação integral 𝑓(𝑥) → (𝑇 𝑓)(𝑥) =

∫︀ 𝑦2𝑦1

d𝑦 𝑘(𝑥, 𝑦)𝑓(𝑦).


Substituindo a Eq. (3.10) na Eq. (3.9) e utilizando a ortogonalidade dos polinô-

mios de Chebyshev de primeira ordem, Eq. (3.8), obtemos

𝑓𝐾𝑃 𝑀 (𝑥) = 1𝜋

√1 − 𝑥2

⎡⎣𝑔0𝜇0 + 2𝑁𝑇 −1∑︁𝑛=1

𝑔𝑛𝜇𝑛𝑇𝑛(𝑥)

⎤⎦ , (3.15)

onde os coeficientes da expansão são dados por

𝜇𝑛 =∫︁ 1

−1d𝑥𝑇𝑛(𝑥)𝑓(𝑥). (3.16)

Na análise de sistemas físicos, ao expandir uma função qualquer da energia em

polinômios de Chebyshev, primeiramente devemos reescalonar a Hamiltoniana de

modo a ajustar o espectro para o intervalo [−1, 1]. Caso contrário, a expansão não

converge. Sejam {|𝜑𝑖⟩} os auto-estados de um Hamiltoniano, satisfazendo 𝐻 |𝜑𝑖⟩ =

𝐸𝑖 |𝜑𝑖⟩ com 𝐸𝑚𝑖𝑛 < 𝐸𝑖 < 𝐸𝑚𝑎𝑥. Aplicando a transformação linear [22]

�̃� = 2𝐻𝑎

− 𝑏

𝑎𝐼, (3.17)

𝐸𝑖 = 2𝐸𝑖

𝑎− 𝑏

𝑎, (3.18)

onde 𝑎 = 𝐸𝑚𝑎𝑥 − 𝐸𝑚𝑖𝑛 e 𝑏 = 𝐸𝑚𝑎𝑥 + 𝐸𝑚𝑖𝑛, temos uma nova distribuição de au-

tovalores {𝐸𝑖} onde −1 < �̃�𝑖 < 1. O deslocamento −𝑏/𝑎 serve apenas para que o

meio da banda de autovalores passe a ser zero. Os limites 𝐸𝑚𝑎𝑥 e 𝐸𝑚𝑖𝑛 são definidos

como 𝐸𝑚𝑎𝑥 = 𝑚𝑎𝑥{𝐸𝑖} + 𝜖 e 𝐸𝑚𝑖𝑛 = 𝑚𝑖𝑛{𝐸𝑖} − 𝜖, com 𝜖 > 0 de modo a garantir

que o espectro não exceda o intervalo [−1, 1]. O parâmetro 𝜖 escala com a resolução

da expansão, 𝜖 ∝ 1/𝑁𝑇 , quando conhecemos os extremos do espectro, 𝑚𝑎𝑥{𝐸𝑖} e

𝑚𝑖𝑛{𝐸𝑖}.

3.3 Aplicações do método KPM

3.3.1 Densidade de estados

A aplicação mais fundamental do método KPM é o cálculo da densidade de estados

(DOS) [23, 29]. Esta é, inclusive, a forma mais conveniente de testar a eficiência do


método. A densidade de estados

𝑁(𝐸) = Tr[𝛿(𝐸 − �̂�)] (3.19)

pode então ser aproximada da seguinte forma:

𝑁𝐾𝑃 𝑀 (𝐸) = 1𝜋

√1 − 𝐸2

⎡⎣𝑔0𝜇0 + 2𝑁𝑇 −1∑︁𝑛=1

𝑔𝑛𝜇𝑛𝑇𝑛(𝐸)

⎤⎦ , (3.20)

após o reescalonamento dos autovalores (a partir daqui vamos omitir o sinal "∼").

Os coeficientes da expansão são dados por

𝜇𝑛 =∫︁ 1

−1d𝐸𝑁(𝐸)𝑇𝑛(𝐸) =

∫︁ 1

−1d𝐸 Tr[𝛿(𝐸 − �̂�)]𝑇𝑛(𝐸) = Tr[𝑇𝑛(�̂�)]. (3.21)

O cálculo do traço de matrizes de dimensões altas é um procedimento computaci-

onalmente exigente. Existe, porém, uma forma bem simples e precisa de realizar

tal tarefa. A aproximação estocástica do traço [22, 23] consiste, primeiramente, em

definir uma base truncada {|𝑟⟩}, com 𝑟 = 0, 1, ..., 𝑅 − 1. Os vetores dessa base

são expressos por combinações lineares aleatórias dos vetores de base do espaço de

Hilbert do sistema,

|𝑟⟩ =𝐷−1∑︁𝑘=0

𝜉𝑟𝑘 |𝑘⟩ , (3.22)

onde os coeficientes {𝜉𝑟𝑘} são números aleatórios independentes e, em geral, 𝜉𝑟𝑘 ∈ C.

Uma maneira conveniente de definir estes números aleatórios é considerar 𝜉𝑟𝑘 = 𝑒𝑖𝜑𝑟𝑘

onde as fases {𝜑𝑟𝑘} são uma distribuição uniforme no intervalo [0, 2𝜋]. Outra opção

é gerar distribuições gaussianas para as partes reais e imaginárias de 𝜉𝑟𝑘. Podemos

até considerar 𝜉𝑟𝑘 ∈ R em situações envolvendo matrizes simétricas reais. De uma

forma ou de outra, o fundamental é que {𝜉𝑟𝑘} satisfaça:

𝜉𝑟𝑘 = 0; (3.23)

𝜉*𝑟𝑘𝜉𝑟′𝑘′ = 𝛿𝑟,𝑟′𝛿𝑘,𝑘′ , (3.24)


para coeficientes complexos ou

𝜉𝑟𝑘 = 0; (3.25)

𝜉𝑟𝑘𝜉𝑟′𝑘′ = 𝛿𝑟,𝑟′𝛿𝑘,𝑘′ , (3.26)

para coeficientes reais. O traço da Eq. (3.21) é então aproximado por

Tr[𝑇𝑛(�̂�)] ≈ 1𝑅

𝑅−1∑︁𝑟=0

⟨𝑟|𝑇𝑛(�̂�) |𝑟⟩ . (3.27)

Definindo |𝑟⟩𝑛 = 𝑇𝑛(�̂�) |𝑟⟩, reescrevemos os coeficientes da expansão de Chebyshev

como

𝜇𝑛 ≈ 1𝑅

𝑅−1∑︁𝑟=0

⟨𝑟|𝑟⟩𝑛, (3.28)

onde, de acordo com a relação de recorrência, Eq. (3.5),

|𝑟⟩𝑛+1 = 2�̂� |𝑟⟩𝑛 − |𝑟⟩𝑛−1 , (3.29)

𝑛 > 0, com |𝑟⟩0 = |𝑟⟩ e |𝑟⟩1 = �̂� |𝑟⟩. Para determinar um número de estados

aleatórios, 𝑅, que forneça um resultado com boa precisão, vamos calcular a flutuação

da estimativa do traço,

(𝛿𝜇𝑛)2 = 𝜇2𝑛 − 𝜇𝑛

2. (3.30)

Podemos considerar 𝜇1 ≈ (1/𝑅)∑︀𝑅−1𝑟=0 ⟨𝑟| �̂� |𝑟⟩ sem perda de generalidade. Seja

𝐻𝑖,𝑗 = ⟨𝑖| �̂� |𝑗⟩ os elementos de matriz da Hamiltoniana, temos [22]

𝜇21 ≈ 1

𝑅2

𝑅−1∑︁𝑟,𝑟′=0

⟨𝑟| �̂� |𝑟⟩ ⟨𝑟′| �̂� |𝑟′⟩

= 1𝑅2

𝑅−1∑︁𝑟,𝑟′=0

𝐷−1∑︁𝑖,𝑗,𝑖′,𝑗′=0

𝜉*𝑟𝑖𝜉𝑟𝑗𝜉*

𝑟′𝑖′𝜉𝑟′𝑗′𝐻𝑖,𝑗𝐻𝑖′,𝑗′

= (Tr �̂�)2 + 1𝑅

⎡⎣Tr(�̂�2) +(︁|𝜉𝑟𝑖|4 − 2

)︁𝐷−1∑︁𝑗=0

𝐻2𝑗𝑗

⎤⎦ , (3.31)

e

𝜇12 ≈

⎛⎝ 1𝑅

𝑅−1∑︁𝑟=0

𝐷−1∑︁𝑖,𝑗=0

𝜉*𝑟𝑖𝜉𝑟𝑗𝐻𝑖,𝑗

⎞⎠2

=(︃

𝐷−1∑︁𝑖=0

𝐻𝑖𝑖

)︃2

= (Tr �̂�)2, (3.32)


onde utilizamos a Eqs. (3.23) e (3.24), e ∑︀𝑖 |𝑖⟩ ⟨𝑖| = ∑︀𝑗 |𝑗⟩ ⟨𝑗| = 𝐼. A Eq. (3.30) é

então reescrita como

(𝛿𝜇1)2 ≈ 1𝑅

⎡⎣Tr(�̂�2) +(︁|𝜉𝑟𝑖|4 − 2

)︁𝐷−1∑︁𝑗=0

𝐻2𝑗𝑗

⎤⎦ . (3.33)

Sendo Tr(�̂�2) de ordem 𝑂(𝐷), o erro relativo 𝛿𝜇1/𝜇1 é da ordem de 𝑂(1/√𝑅𝐷).

Consequentemente, manter 𝑅 ≪ 𝐷 garante uma boa aproximação desde que, evi-

dentemente, estejamos tratando de sistemas com muitas dimensões 𝐷. Inclusive,

este deve ser o caso para a distribuição {𝜉𝑟𝑘} satisfazer as Eqs. (3.23) à (3.26) da

melhor forma possível. Veremos no Cap. 4 que 𝑅 < 10 já fornece bons resultados

para um sistema da ordem de 105 dimensões.

3.3.2 Evolução temporal

Nesta seção vamos mostrar como os polinômios de Chebyshev podem ser utilizados

para tratar da evolução dinâmica de sistemas envolvendo inúmeros graus de liber-

dade [27, 30, 31, 32, 33]. Seja 𝑓(𝐸) uma função já reescalonada para o intervalo

[−1, 1]. Como anteriormente, podemos a expandir como

𝑓𝐾𝑃 𝑀 (𝐸) = 1𝜋

√1 − 𝐸2

⎡⎣𝑔0𝜇0 + 2𝑁𝑇 −1∑︁𝑛=1


⎤⎦ , (3.34)

onde os coeficientes {𝜇𝑛} são fornecidos por

𝜇𝑛 =∫︁ 1

−1d𝐸 𝑓(𝐸)𝑇𝑛(𝐸). (3.35)

Vamos reescrever a expansão de outra forma fazendo a transformação [24]

𝑓(𝐸) → 𝑓(𝐸)/(𝜋√︀

1 − 𝐸2), (3.36)

e então a Eq. (3.34) se torna

𝑓𝐾𝑃 𝑀 (𝐸) = 𝑔0𝛾0 + 2𝑁𝑇 −1∑︁𝑛=1

𝑔𝑛𝛾𝑛𝑇𝑛(𝐸), (3.37)

onde os novos coeficientes são dados por

𝛾𝑛 =∫︁ 1

−1d𝐸 𝑓(𝐸)𝑇𝑛(𝐸)

𝜋√

1 − 𝐸2. (3.38)


Esta é apenas uma forma equivalente de escrever a expansão, mais conveniente para

tratar da evolução temporal como veremos adiante. Antes, vejamos um pouco de

mecânica quântica.

A equação de Schrödinger para o operador evolução temporal 𝑈(𝑡, 𝑡0) é dada por

(} = 1)

𝑖𝜕

𝜕𝑡𝑈(𝑡, 𝑡0) = 𝐻𝑈(𝑡, 𝑡0). (3.39)

Seja |𝜓(𝑡0)⟩ o estado num instante inicial 𝑡0, se soubermos a maneira pela qual o

operador 𝑈(𝑡, 𝑡0) satisfaz

𝑈(𝑡, 𝑡0) |𝜓(𝑡0)⟩ = |𝜓(𝑡)⟩ , (3.40)

podemos obter o auto-estado em qualquer instante 𝑡. Devemos então encontrar as

soluções para a Eq. (3.39). Há três casos a serem considerados [19]. No primeiro,

quando o operador Hamiltoniano é independente do tempo, a solução é simplesmente

dada por

𝑈(𝑡, 𝑡0) = 𝑒−𝑖𝐻(𝑡−𝑡0). (3.41)

No segundo caso, o Hamiltoniano é dependente do tempo e [𝐻(𝑡1), 𝐻(𝑡2)] = 0 para

qualquer 𝑡1 ̸= 𝑡2. A solução é então escrita como

𝑈(𝑡, 𝑡0) = 𝑒−𝑖∫︀ 𝑡

𝑡0d𝑡′ 𝐻(𝑡′)

. (3.42)

Quando [𝐻(𝑡1), 𝐻(𝑡2)] ̸= 0, no terceiro caso, a solução é similar a Eq. (3.42),

mas como nem todos os Hamiltonianos comutam, não podemos mais escrever a

solução como uma exponencial apenas. Introduzimos então o chamado operador

ordenamento temporal 𝑇 e a solução se torna

𝑈(𝑡, 𝑡0) = 𝑇𝑒−𝑖∫︀ 𝑡

𝑡0d𝑡′ 𝐻(𝑡′)

, (3.43)

que é uma expansão em séries de Dyson. Maiores detalhes não serão necessários

já que o método KPM para a evolução temporal que será formulado a seguir acaba

por utilizar um operador no formato da Eq. (3.41) mesmo quando a Hamiltoniana


possui dependência temporal.

Vamos considerar que 𝑓(𝐸) = 𝑒−𝑖𝐸𝜏 . Os coeficientes 𝛾𝑛 se tornam

𝛾𝑛 =∫︁ 1

−1d𝐸 𝑒−𝑖𝐸𝜏𝑇𝑛(𝐸)

𝜋√

1 − 𝐸2= (−𝑖)𝑛𝐽𝑛(𝜏), (3.44)

onde

𝐽𝑛(𝜏) = 12𝜋

∫︁ 𝜋

−𝜋d𝜃 𝑒−𝑖(𝑛𝜃−𝜏cos𝜃) (3.45)

é a n-ésima função de Bessel de primeira espécie. Estas satisfazem a relação de

recorrência

𝐽𝑛+1(𝑥) = 2𝑛𝑥𝐽𝑛(𝑥) − 𝐽𝑛−1(𝑥). (3.46)

Podemos agora calcular a evolução temporal de um estado inicial |𝜓(𝑡0)⟩ para um

instante 𝑡, 𝑡− 𝑡0 = 𝜏 , de acordo com (após o reescalonamento da Hamiltoniana)

|𝜓(𝑡)⟩ = 𝑈(𝜏) |𝜓(𝑡0)⟩ = 𝑒−𝑖𝐻𝜏 |𝜓(𝑡0)⟩

= 𝑔0𝛾0(𝜏) |𝜓(𝑡0)⟩ + 2𝑁𝑇 −1∑︁𝑛=1

𝑔𝑛𝛾𝑛(𝜏)𝑇𝑛(�̂�) |𝜓(𝑡0)⟩

= 𝑔0𝐽0(𝜏) |𝜓(𝑡0)⟩0 + 2𝑁𝑇 −1∑︁𝑛=1

𝑔𝑛𝐽𝑛(𝜏) |𝜓(𝑡0)⟩𝑛 , (3.47)

onde |𝜓(𝑡0)⟩𝑛 = 𝑇𝑛(�̂�) |𝜓(𝑡0)⟩ e

|𝜓(𝑡0)⟩𝑛+1 = 2𝐻 |𝜓(𝑡0)⟩𝑛 − |𝜓(𝑡0)⟩𝑛−1 (3.48)

para 𝑛 > 1 com |𝜓(𝑡0)⟩0 = |𝜓(𝑡0)⟩ e |𝜓(𝑡0)⟩1 = 𝐻 |𝜓(𝑡0)⟩, de acordo com a relação

de recorrência dos polinômios de Chebyshev, Eq. (3.5).

A aproximação descrita anteriormente é válida para Hamiltonianas independen-

tes do tempo. Se este não for o caso, podemos dividir o tempo em 𝐾 pequenas

partes Δ𝑡 = 𝑡/𝐾 de maneira que

𝑈(𝑡) =𝐾−1∏︁𝜆=0

𝑆((𝜆+ 1)Δ𝑡, 𝜆Δ𝑡), (3.49)


onde

𝑆((𝜆+ 1)Δ𝑡, 𝜆Δ𝑡) = 𝑇𝑒−𝑖∫︀ 𝑡+Δ𝑡

𝑡d𝑡′ 𝐻(𝑡′)

≈ 𝑒−𝑖𝐻(𝑡)Δ𝑡, (3.50)

caso Δ𝑡 → 0. Evidentemente, a Hamiltoniana deve ser praticamente constante

em escalas temporais da ordem Δ𝑡. Se deve, então, escolher o 𝐾 de maneira que

𝐾 ≫ 𝑡/𝜏 onde 𝜏 é a escala de tempo induzida por 𝐻(𝑡). Sendo assim, expandimos

o operador evolução temporal em polinômios de Chebyshev para cada Δ𝑡 da forma

[30, 32, 33]

𝑒−𝑖𝐻(𝑡)Δ𝑡 = 𝑔0𝐽0(Δ𝑡) + 2𝑁𝑇 −1∑︁𝑛=1

(−𝑖)𝑛𝑔𝑛𝐽𝑛(Δ𝑡)𝑇𝑛(�̂�(𝑡)). (3.51)

Antes de prosseguir para o próximo capítulo, seria conveniente mostrar o efeito

do reescalonamento da Hamiltoniana no operador evolução temporal. Seja

𝑈(𝑡) = 𝑒−𝑖𝐻𝑡/~. (3.52)

De acordo com a Eq. 3.17 reescrevemos a equação acima como

𝑈(𝑡) = 𝑒−𝑖𝑡𝑎/2~ (�̃�+𝑏/𝑎)

= 𝑒−𝑖�̃�𝑎𝑡/2~ 𝑒−𝑖𝑏𝑡/2~

= 𝑒−𝑖�̃�𝑡′ (3.53)

onde consideramos 𝑒−𝑖𝑏𝑡/2~ = 1 e 𝑡′ = 𝑎𝑡/2~, lembrando que 𝑎 = 𝐸𝑚𝑎𝑥 −𝐸𝑚𝑖𝑛. Esta

é, portanto, a forma do operador no qual aplicamos o método KPM.

Capítulo 4

Resultados e discussão

Neste capítulo mostramos os resultados obtidos aplicando o método KPM para o

cálculo da densidade de estados (DOS) e do decaimento da fidelidade considerando

o código tórico 𝐿 = 3, ou seja, uma rede bidimensional de 2𝐿2 = 18 spins com

condições de contorno periódicas.

Como já foi mencionado no Cap. 1, as propriedades topológicas de auto-correção

do código não serão levadas em conta. Sendo assim, não há dinâmica entre os spins

(𝐻𝑆 = 0), apenas uma interação efetiva entre os pares destes,

𝑉𝑒𝑓𝑓 = 12∑︁𝑗 ̸=𝑙

(︃𝐽𝑥

|x𝑗 − x𝑙|2+𝛽𝑋𝑗𝑋𝑙 + 𝐽𝑧

|x𝑗 − x𝑙|2+𝛽𝑍𝑗𝑍𝑙

)︃, (4.1)

induzida pelo acoplamento spin-bóson. O operador 𝑍 está representado na base

|0⟩ , |1⟩ de um spin 1/2 de modo que 𝑍 |0⟩ = |0⟩ e 𝑍 |1⟩ = − |1⟩ e então o opera-

dor 𝑋 satisfaz 𝑋 |0⟩ = |1⟩ e 𝑋 |1⟩ = |0⟩. O espaço de Hilbert do sistema possui

218 dimensões1. Os estados da base {|𝑎1⟩ ⊗ |𝑎2⟩ ⊗ ... ⊗ |𝑎18⟩}, 𝑎𝑖 ∈ {0, 1}, serão

representados por |𝑘⟩ onde

𝑘 = 𝑎120 + 𝑎221 + ..., + 𝑎18217, (4.2)

e portanto os estados com todos os spins up ou down seriam |0⟩ ≡ |000...0⟩ e

|262143⟩ ≡ |111...1⟩, respectivamente. Os spins da rede estão sujeitos a inversões de

1Os operadores serão produtos tensoriais dos operadores de um spin. Supondo que o spin 3 foiinvertido. Este operador é dado por 𝐼1 ⊗ 𝐼2 ⊗ 𝑋3 ⊗ 𝐼4 ⊗ ... ⊗ 𝐼18 onde 𝐼𝑖 é o operador identidade noi-ésimo spin.

Capítulo 4. Resultados e discussão 28

fase e de bit. Não há necessidade de considerar interações 𝑌𝑗𝑌𝑙 já que 𝑌 = 𝑖𝑋𝑍. De

acordo com a Eq. (4.1), é fácil verificar que o canal 𝑋 gera apenas elementos fora da

diagonal da Hamiltoniana e 𝑍 gera os elementos da diagonal. Estaremos assumindo

a constante de rede 𝑎 igual a 1.

Os dados foram obtidos através de programas na linguagem FORTRAN 90

(Apêndice A) e os gráficos plotados utilizando o Origin 6.0 (OriginLab, Northamp-

ton, MA). Na tabela 4.1, estão descritos todos os parâmetros levados em consideração

nas análises que serão apresentadas a seguir.

Tabela 4.1: Parâmetros utilizados nas análises da densidade de estados e/ou dafidelidade quântica.

Parâmetro Descrição𝑁𝑇 ordem da expansão de Chebyshev𝑅 número de bases truncadas {|𝑟⟩} na estimativa de 𝜇𝑛

𝛽 determina o alcance da interação𝐽𝑥 intensidade da interação do canal 𝑋𝐽𝑧 intensidade da interação do canal 𝑍Δ𝐸 largura (geralmente estimada) da banda energética|𝜓(0)⟩ estado inicialKernel Dirichlet, Jackson ou Lorentz𝑡 tempo (adimensional)𝑡* = 2~𝑡/Δ𝐸𝐾 número de passos no tempoΔ𝑡 = 𝑡/𝐾

4.1 Cálculo da densidade de estados

A densidade de estados é fornecida por

𝑁𝐾𝑃 𝑀 (𝐸) = 1𝜋

√1 − 𝐸2

⎡⎣𝑔0𝜇0 + 2𝑁𝑇 −1∑︁𝑛=1


⎤⎦ (4.3)

onde 𝜇𝑛 ≈ (1/𝑅)∑︀𝑅−1𝑟=0 ⟨𝑟|𝑟⟩𝑛 (aproximação estocástica do traço) sendo |𝑟⟩ uma

combinação linear dos vetores de base |𝑘⟩ com coeficientes {𝜉𝑟𝑘}.

Inicialmente vamos supor que 𝐽𝑥 = 0, ou seja, a Hamiltoniana é diagonal. Dessa

forma, já conhecemos o espectro e torna-se mais fácil selecionar os parâmetros𝑁𝑇 e𝑅


mais convenientes. A Fig. 4.1 mostra a densidade de estados (DOS) para diferentes

valores de 𝑁𝑇 comparados com a curva exata2. A curva está normalizada para

área unitária (idem para todas as figuras envolvendo o DOS nesta seção). Nota-se

que quanto maior a ordem 𝑁𝑇 , ocorrem pequenas oscilações. Este efeito é comum

e ocorre porque para um sistema discreto e finito a densidade de estados é uma

sequência de funções delta, Eq. (3.19). Basicamente, os polinômios de Chebyshev

estão tentando resolver a estrutura de níveis do sistema. O fato de observarmos uma

curva contínua para ordens mais baixas tem a ver com uma suavização inerente ao

método KPM. O ganho varia inversamente com a maior ordem polinomial usada.

Quando esta ordem é muito alta, portanto, a expansão passa a resolver a estrutura de

níveis discretos do sistema e então surgem as oscilações. Podemos observá-las duma

forma mais nítida na Fig. 4.2 onde empregamos o kernel de Dirichlet (𝑔𝑛 = 1). Em

geral, esta é a razão pela qual a aplicação deste kernel é recomendável (as vezes

necessária) apenas em expansões fortemente convergentes, como veremos na Sec.

4.2 para o operador evolução temporal.

Note que até agora foi utilizado 𝑅 = 3 somente. Este valor é suficiente, tendo em

-10 0 10 200.00

0.02

0.04

0.06

0.08

0.10

0.12

0.14

0.16

DO

S

E/Jz

curva exata

NT = 50

NT = 100

NT = 200

Figura 4.1: Densidade de estados para diferentes valores de 𝑁𝑇 . A curva está normalizadapara área unitária. Parâmetros: 𝑅 = 3, Δ𝐸/𝐽𝑧 = 44,21, 𝛽 = 0 e kernel de Jackson.

2A curva exata é, na verdade, uma curva ajustada ao histograma dos autovalores da Hamilto-niana.


-10 0 10 200.00

0.02

0.04

0.06

0.08

0.10

0.12

0.14

0.16

DO

S

E/Jz

Dirichlet

Jackson

Lorentz

Figura 4.2: Curva da densidade de estados comparando os kérneis de Dirichlet, Jackson eLorentz. Este último não apresenta boa convergência nos extremos do espectro e daqui emdiante vamos utilizar o kernel de Jackson. Parâmetros: 𝑅 = 3, Δ𝐸/𝐽𝑧 = 44,21, 𝛽 = 0 e𝑁𝑇 = 100.

vista que o erro relativo, 𝛿𝜇𝑛/𝜇𝑛 é da ordem de 𝑂(1/√𝑅𝐷) (Sec. 3.3.1). Poderíamos

mostrar uma comparação de curvas com diferentes valores de 𝑅 mas as variações

seriam imperceptíveis. Provavelmente, mesmo com 𝑅 ≫ 3 as mudanças seriam

pouco significativas não compensando o custo computacional.

Uma boa maneira de verificar a confiabilidade deste procedimento seria consi-

derar apenas os elementos fora da diagonal da Hamiltoniana gerados pela interação

𝑋𝑗𝑋𝑙 fazendo 𝐽𝑧 = 0. Os resultados devem corresponder aos anteriores com 𝐽𝑥 = 0,

se mantermos a mesma escala de energia, evidentemente (Fig. 4.3). Isto ocorre

devido à invariância rotacional [19] das propriedades deste sistema. Para melhor vi-

sualização, vamos considerar o estado de apenas um qubit. Transformar o operador

𝑍 para 𝑋 é equivalente a gerar uma rotação no sistema,

|0⟩ → 1√2

(|0⟩ + |1⟩) ≡ |+⟩ , (4.4)

|1⟩ → 1√2

(|0⟩ − |1⟩) ≡ |−⟩ , (4.5)

onde, portanto, 𝑍 |+⟩ = |−⟩ e 𝑍 |−⟩ = |+⟩. É claro que, no caso da Fig. 4.3,

tivemos a vantagem de já conhecer previamente o espectro. Se esta não fosse a


-10 0 10 200,00

0,02

0,04

0,06

0,08

0,10

0,12

0,14

0,16

DO

S

E/Jα

curva exata

Jx=0

Jz=0

Figura 4.3: Densidade de estados para os casos onde só há interações em 𝑋 (𝐽𝑧 = 0) ouem 𝑍 (𝐽𝑥 = 0). A energia está expressa em unidades de 𝐽𝛼 onde 𝐽𝛼 = 𝐽𝑥 se 𝐽𝑧 = 0 evice-e-versa. Parâmetros: 𝑅 = 3, Δ𝐸/𝐽𝛼 = 44,21, 𝛽 = 0, 𝑁𝑇 = 100 e kernel de Jackson.

situação poderíamos estimá-lo até o ponto onde a expansão de Chebyshev passa a

divergir. A Fig. 4.4 mostra estas estimativas partindo de larguras de banda maiores

até valores próximos de Δ𝐸 real. Na Fig. 4.5 podemos notar a divergência dos

-10 0 10 200.00

0.02

0.04

0.06

0.08

0.10

0.12

0.14

0.16

0.18

DO

S

E/Jx

∆E/Jx = 44,19 (valor exato)

∆E/Jx = 46,19

∆E/Jx = 54,19

∆E/Jx = 64,19

Figura 4.4: Estimando a densidade de estados de acordo com Δ𝐸. Parâmetros: 𝑅 = 3,𝛽 = 0, 𝑁𝑇 = 100 e kernel de Jackson.


0 10 20 30 40-2x10

5

-1x105

0

1x105

2x105

3x105

µ(n

)

n

∆E/Jx = 44,19 (valor exato)

∆E/Jx = 44,21

∆E/Jx = 42,19

Figura 4.5: Divergência dos coeficientes de expansão 𝜇𝑛. Parâmetros: 𝑅 = 3, 𝐽𝑧 = 0, 𝛽 = 0e 𝑁𝑇 = 40.

coeficientes da expansão 𝜇𝑛 para Δ𝐸 pouco menor do que valor exato. A divergência

ocorre, evidentemente, quando alguns autovalores (reescalonados) são encontrados

fora do intervalo [−1, 1].

Por fim, vamos considerar os dois canais de interação com 𝐽𝑥 = 𝐽𝑧 = 𝐽 . Na

Fig. 4.6 podemos ver que, ao aumentar o alcance da interação, ou seja, diminuindo

𝛽 (caso sub-ôhmico), o espectro tende para maiores valores de energia, como era

de se esperar. A Fig. 4.7 mostra o comportamento da densidade de estados com a

variação da intensidade da interação 𝐽 . Este parâmetro define linearmente a escala

energética.

As curvas da densidade de estados do sistema obtidas nesta seção mostra que

o método KPM é bastante conveniente para tratar deste aspecto. Em relação a

métodos que envolvem a diagonalização direta da Hamiltoniana, o método KPM

exige, significativamente, menos recursos computacionais.


-10 0 10 20 300,00

0,02

0,04

0,06

0,08

0,10

DO

S

E/J

β = -1

β = 0

β = 1

Figura 4.6: Densidade de estados, considerando ambos os canais de interação para diferentesvalores de 𝛽. Parâmetros: 𝑅 = 3, 𝑁𝑇 = 100 e kernel de Jackson.

-20 -10 0 10 20 300.00

0.02

0.04

0.06

0.08

0.10

DO

S

E/J*

J* = J

J* = J / 1,25

J* = J / 1,50

J* = J / 1,75

Figura 4.7: Densidade de estados para diferentes valores de 𝐽 . Parâmetros: 𝑅 = 3, 𝛽 = 1,𝑁𝑇 = 100 e kernel de Jackson.


4.2 Análise da fidelidade quântica

A aplicação do método KPM que vem a seguir é mais voltada para o contexto da

Informação Quântica. Vamos calcular o decaimento da fidelidade do sistema. A

fidelidade é uma medida que informa o quão sensível é a dinâmica de um sistema

quântico sujeito a interações com o ambiente externo. Seria então uma medida de

distância entre um estado quântico inicial e o estado após a evolução temporal. Ela

é definida como [3]

𝐹 (𝜌1, 𝜌2) ≡ Tr√︁𝜌

1/21 𝜌2𝜌

1/21 , (4.6)

onde 𝜌1 e 𝜌2 são dois estados arbitrários. No caso de dois estados puros, digamos, o

estado inicial |𝜓(0)⟩ e o estado |𝜓(𝑡)⟩ após um certo instante 𝑡, teremos

𝐹 (|𝜓(0)⟩ , |𝜓(𝑡)⟩) = Tr[︂√︁

⟨𝜓(0)|𝜓(𝑡)⟩ ⟨𝜓(𝑡)|𝜓(0)⟩ |𝜓(0)⟩ ⟨𝜓(0)|]︂

=√︁

⟨𝜓(0)|𝜓(𝑡)⟩ ⟨𝜓(𝑡)|𝜓(0)⟩

= |⟨𝜓(0)|𝜓(𝑡)⟩|. (4.7)

Portanto, a fidelidade nada mais é do que o valor absoluto do produto interno entre

os dois estados. Seja |𝜓(𝑡)⟩ = 𝑈(𝑡) |𝜓(0)⟩, o operador evolução temporal é expresso

por, de acordo com o método KPM,

𝑈(𝑡) = 𝑔0𝐽0(𝑡) + 2𝑁𝑇 −1∑︁𝑛=1

(−𝑖)𝑛𝑔𝑛𝐽𝑛(𝑡)𝑇𝑛(�̂�). (4.8)

Para análise do decaimento, a fidelidade deve ser calculada em diversos instantes

e portanto será mais conveniente dividir o tempo 𝑡 em 𝐾 partes de Δ𝑡. Note que

este procedimento é similar ao o descrito para Hamiltonianas dependentes do tempo

na Sec. 3.3.2. Mas como 𝑉𝑒𝑓𝑓 é independente do tempo, torna-se desnecessário

requerer que Δ𝑡 → 0. Considerando diversos intervalos de tempo 𝑡𝜆 = 𝜆Δ𝑡, onde

𝜆 = 0, 1, ...,𝐾, temos que

⟨𝜓(0)|𝜓(𝑡𝜆+1)⟩ = ⟨𝜓(0)| 𝑒−𝑖�̂�Δ𝑡 |𝜓𝑡𝜆⟩

= 𝑔0𝐽0(Δ𝑡)𝛾0𝜆 + 2

𝑁𝑇 −1∑︁𝑛=1

(−𝑖)𝑛𝑔𝑛𝐽𝑛(Δ𝑡)𝛾𝑛𝜆 , (4.9)


onde 𝛾𝑛𝜆 = ⟨𝜓(0)|𝑇𝑛(�̂�) |𝜓(𝑡𝜆)⟩. Sendo assim, a fidelidade a cada passo no tempo é

obtida por

𝐹𝜆+1 = 𝐹 (|𝜓(0)⟩ , |𝜓(𝑡𝜆+1)⟩) = | ⟨𝜓(0)|𝜓(𝑡𝜆+1)⟩|. (4.10)

Considerando o estado inicial |𝜓(0)⟩ = |0⟩ e 𝐽𝑧 = 𝐽𝑥 = 𝐽 a Fig. 4.8 mostra o

comportamento da fidelidade ao longo do tempo para diferentes valores de 𝛽. O

tempo 𝑡 = 𝑡*Δ𝐸/2~ é adimensional. Torna-se claro, portanto, que a escala de tempo

do decaimento da fidelidade depende da largura da banda Δ𝐸. Quanto menor for

𝛽, maior Δ𝐸 e mais rapidamente a fidelidade decai. Isso está associado à correlação

do sistema, ou seja, ao alcance das interações entre os pares de spins. As Figs. 4.9 e

4.10 fornecem a mesma análise considerando os estados iniciais (|0⟩+|262143⟩)/√

2 e

(|0⟩+|1⟩+...+|262143⟩)/√𝐷, respectivamente. Outro parâmetro que tem influência

na escala de tempo do decaimento, é a intensidade da interação 𝐽 como podemos

observar na Fig. 4.11. Isso é evidente tendo em vista a dependência em Δ𝐸.

A única razão pela qual há oscilações na fidelidade ao longo do tempo é que no

modelo de interação o sistema mantem coerência quântica. Se levássemos em conta a

0 10 20 30 40 50 600,0

0,2

0,4

0,6

0,8

1,0

F(t

)

t

β = 0,5

β = 0,25

β = 0

β = -0,25

β = -0,5

Figura 4.8: Decaimento da fidelidade para diferentes valores de 𝛽. Parâmetros: |𝜓(0)⟩ = |0⟩,Δ𝑡 = 3, 𝑁𝑇 = 15 e kernel de Dirichlet.


0 10 20 30 40 50 60

0,0

0,2

0,4

0,6

0,8

1,0

β = 1

β = 0,5

β = 0

β = -0,5

β = -1

F(t

)

t

Figura 4.9: Decaimento da fidelidade para diferentes valores de 𝛽. Parâmetros: |𝜓(0)⟩ =(|0⟩ + |262143⟩)/

√2, Δ𝑡 = 3, 𝑁𝑇 = 15 e kernel de Dirichlet.

0 10 20 30 40 50 600,0

0,2

0,4

0,6

0,8

1,0

β = 1

β = 0,5

β = 0

β = -0,5

β = -1

F(t

)

t

Figura 4.10: Decaimento da fidelidade para diferentes valores de 𝛽. Parâmetros: |𝜓(0)⟩ =(|0⟩ + |1⟩ + ...+ |262143⟩)/

√𝐷, Δ𝑡 = 3, 𝑁𝑇 = 15 e kernel de Dirichlet.


0,0 0,2 0,4 0,6 0,8 1,0 1,20,0

0,2

0,4

0,6

0,8

1,0

F(t

*)

t*

J = J'

J = 1,25 J'

J = 1,50 J'

J = 1,75 J'

Figura 4.11: Decaimento da fidelidade para diferente valores de 𝐽 . O tempo 𝑡* está expressoem unidades de ~/𝐽 . Observe que na curva onde a fidelidade decai mais rápido (maior 𝐽)percebe-se um maior número de pontos (maior resolução). Isso ocorreu porque calculamostodos estes casos para um mesmo 𝑡 e 𝐾, e ajustamos a escala temporal (eixo 𝑥) posterior-mente. Parâmetros: |𝜓(0)⟩ = (|0⟩ + |262143⟩)/

√2, 𝐽𝑧 = 𝐽𝑥 = 𝐽 , 𝛽 = 1, 𝑁𝑇 = 15 e kernel

de Dirichlet.

dinâmica do banho, ou seja, o fato de que as interações efetivas não são instantâneas

e flutuam com o tempo, o sistema não evoluiria de forma unitária. Assim, o estado

deixaria de ser puro, haveria descoerência da fase e as oscilações desapareceriam

rapidamente.

Nas Figs. 4.8 a 4.11, a curva não deve se estabilizar em um valor finito da fide-

lidade, mas sim oscilar com uma amplitude que vai decrescendo com o tempo. Para

escalas de tempo muito longas, é possível que haja uma recuperação da fidelidade,

mas como o numero de estados do sistema é muito grande, essa escala deve ser

exponencialmente longa. Neste sentido, a única parte da curva do decaimento da

fidelidade em que se pode confiar é a inicial.

O operador evolução temporal é unitário e não deve alterar a norma do vetor de

estado. Verificar o valor da norma em função do tempo é, portanto, uma maneira de

se certificar que a aproximação está funcionando. O kernel de Dirichlet satisfaz este

requerimento da melhor forma possível. O kernel de Jackson, porém, não é indicado


em algumas expansões rapidamente convergentes (𝑁𝑇 = 15) como é caso. Note que

a expansão também envolve as funções de Bessel e para garantir uma resolução,

no mínimo razoável, Δ𝑡 deve valer algo entre 3 e 5. Então para baixas precisões

numéricas, 𝑔𝑛𝐽𝑛(Δ𝑡) = 0 logo nos primeiros termos da expansão. O mesmo é válido

para o kernel de Lorentz. A Fig. 4.12 mostra o comportamento da norma do estado

após cada passo Δ𝑡.

0 10 20 30 40 50 600,0

0,2

0,4

0,6

0,8

1,0

F(t

)

t

Dirichlet

Jackson

Lorentz

Figura 4.12: Comparação entre o kérneis de Dirichlet, Jackson e Lorentz no que se refereà conservação da norma do estado inicial. A linha contínua é o decaimento da fidelidade etracejada representa a norma. Parâmetros: |𝜓(0)⟩ = (|0⟩ + |262143⟩)/

√2, 𝑁𝑇 = 15.

Em suma, nesta seção mostramos explicitamente a dependência da escala de

tempo do decaimento da fidelidade em relação à correlação espacial do sistema.

Capítulo 5

Conclusões

O Método Kernel Polinomial (KPM) mostrou ser uma forma numericamente con-

veniente na análise de códigos de superfície. Os polinômios de Chebyshev possuem

boas propriedades de convergência e uma simples relação de recorrência.

Consideramos o código tórico 𝐿 = 3 com uma interação efetiva entre todos os pa-

res de spins. Contudo, o formalismo adotado neste trabalho pode ser estendido para

outros códigos quânticos de superfície, inclusive considerando a proteção topológica.

É importante deixar claro que ao aumentar o tamanho do sistema, necessita-se de

muito mais recursos computacionais.

A principal razão de termos calculado a densidade de estados do sistema foi para

mostrar a funcionalidade do método KPM, analisando o comportamento de seus

principais parâmetros. Em seguida, mostramos que a escala de tempo associada

ao decaimento da fidelidade depende da correlação espacial da interação entre os

pares de spins da rede, sendo esta associada ao alcance da interação e à intensidade

do acoplamento. Quanto maior for a correlação do sistema, mais rapidamente a

fidelidade decai.

A correlação espacial da interação entre os spins da rede faz com que a acumula-

ção de erros nos qubits do código tórico passem a não seguir um processo Markovi-

ano. Sendo assim, a probabilidade de erro não é constante e torna-se mais complexa

a estimativa de um limiar, se houver, para a computação quântica tolerante a falhas.

Referências Bibliográficas

[1] P. W. Shor, "Algorithms for quantum computation: discrete logarithms and fac-

toring", Proceedings of the 35th Annual IEEE Symposium on the Foundations

of Computer Science (1994).

[2] C. Zalka, "Efficient Simulation of Quantum Systems by Quantum Computers",

Proc. Roy. Soc. Lond. A 454, 313 (1998).

[3] M. A. Nielsen and I. L. Chuang, Quantum Computation and Quantum Infor-

mation. Cambridge University Press, Cambridge, UK, 2000.

[4] W. H. Zurek, "Decoherence and the transition from quantum to classical", Phys.

Today 44, 36 (1991).

[5] H.-P. Breuer and F. Petruccione, The Theory of Open Quantum Systems. Ox-

ford University Press, Oxford, UK, 2002.

[6] A. R. Calderbank and P. W. Shor, "Good quantum error-correcting codes exist",

Phys. Rev. A 54, 1098 (1996); D. Gottesman, "Theory of fault-tolerant quantum

computation", Phys. Rev. A 57, 127 (1998).

[7] E. Knill, R. Laflamme, and W. H. Zurek, "Resilient Quantum Computation",

Science 279, 342 (1998); D. Aharonov and M. Ben-Or, "Fault-Tolerant Quan-

tum Computation With Constant Error Rate", e-print arXiv:quant-ph/9906129.

[8] E. Novais, E. R. Mucciolo, and H. U. Baranger, "Resilient Quantum Computa-

tion in Correlated Environments: A Quantum Phase Transition Perspective",

Phys. Rev. Lett. 98, 040501 (2007).


[9] B. M. Terhal and G. Burkard, "Fault-tolerant quantum computation for local

non-Markovian noise", Phys. Rev. A 71, 012336 (2005).

[10] H. K. Ng and J. Preskill, "Fault-tolerant quantum computation versus Gaussian

noise", Phys. Rev. A 79, 032318 (2009).

[11] D. Aharonov, A. Kitaev, and J. Preskill, "Fault-tolerant quantum computation

with long-range correlated noise", Phys. Rev. Lett. 96, 050504 (2006).

[12] E. Novais, E. R. Mucciolo, and H. U. Baranger, "Bound on quantum computa-

tion time: Quantum error correction in a critical environment", Phys. Rev. A

82, 020303(R) (2010).

[13] E. Dennis, A. Kitaev, A. Landahl, and J. Preskill, "Topological quantum me-

mory", J. Math. Phys. 43, 4452 (2002).

[14] A. Kitaev, "Fault-tolerant quantum computation by anyons", Ann. Phys. 303,

2 (2003).

[15] S. Gladchenko, D. Olaya, E. Dupont-Ferrier, B. Douçot, L. B. Ioffe, and M. E.

Gershenson, "Superconducting nanocircuits for topologically protected qubits",

Nat. Phys. 5, 48 (2009).

[16] X. Yao et al., "Experimental demonstration of topological error correction",

Nature 482, 489 (2012).

[17] A. J. Leggett et al., "Dynamics of the dissipative two-state system", Rev. Mod.

Phys. 59, 1 (1987).

[18] H. Bruus and K. Flensberg, Many-Body Quantum Theory in Condensed Matter

Physics. Oxford University Press, Oxford, UK, 2004.

[19] J. J. Sakurai, Modern Quantum Mechanics. Addison Wesley, 1993.

[20] M. Maggiore, A Modern Introduction to Quantum Field Theory. Oxford Uni-

versity Press, Oxford, UK, 2004.


[21] I. S. Gradshteyn and I. M. Ryzhik, Table of Integrals, Series, and Products.

Academic Press, 5th edition, 1994.

[22] A. Weisse, G. Wellein, A Alvermann, and H. Fehske, "The kernel polynomial

method", Rev. Mod. Phys. 78, 275 (2006).

[23] R. N. Silver and H. Röder, "Calculation of densities of states and spectral

functions by Chebyshev recursion and maximum entropy", Phys. Rev. E 56,

4822 (1997).

[24] E. Mucciolo, "Time Evolution Calculations using the Kernel Polynomial

Method"(não publicado).

[25] L. Covaci, F. M. Peeters, and M. Berciu, "Efficient Numerical Approach

to Inhomogeneous Superconductivity: The Chebyshev-Bogoliubov–de Gennes

Method", Phys. Rev. Lett. 105, 167006 (2010).

[26] L. W. Wang, "Calculating the density of states and optical-absorption spectra

of large quantum systems by the plane-wave moments method", Phys. Rev. B

49, 10154 (1994).

[27] Y. L. Loh, S. N. Taraskin, and S. R. Elliott, "Fast Time-Evolution Method for

Dynamical Systems", Phys. Rev. Lett. 84, 2290 (2000).

[28] G. B. Arfken, Mathematical Methods for Physicists. Academic Press, 6th edi-

tion, 2005.

[29] H. Röder, R. N. Silver, D. A. Drabold, and Jian Jun Dong, "Kernel polyno-

mial method for a nonorthogonal electronic-structure calculation of amorphous

diamond", Phys. Rev. B 55, 15382 (1997).

[30] H. Tal-Ezer and R. Kosloff, "An accurate and efficient scheme for propagating

the time dependent Schrodinger equation", J. Chem. Phys. 81, 3967 (1984).

[31] H. Fehske, J. Schleede, G. Schubert, G. Wellein, V. S. Filinov, and A. R. Bishop,

"Numerical approaches to time evolution of complex quantum systems", Phys.

Lett. A 373, 2182 (2009).


[32] H. Ishii, N. Kobayashi, and K. Hirose, "Quantum transport properties of carbon

nanotube field-effect transistors with electron-phonon coupling", Phys. Rev. B

76, 205432 (2007).

[33] S. Roche, J. Jiang, F. Triozon, and R. Saito, "Quantum Dephasing in Carbon

Nanotubes due to Electron-Phonon Coupling", Phys. Rev. Lett. 95, 076803

(2005).

Apêndice A

Programas em FORTRAN 90

Fornecemos neste apêndice uma lista dos principais programas e sub-rotinas em

FORTRAN 90 utilizados no cálculo da densidade de estados e evolução temporal.

Em seguida, disponibilizamos os códigos mais relevantes.

A.1 Programas

• mu_coeff. Calcula os coeficientes {𝜇𝑛} da expansão de Chebyshev para a

densidade de estados considerando os canais 𝑋 e 𝑍.

• dos_cheb. Fornece a densidade de estados.

• fidelity_cheb. Calcula a evolução temporal do estado e a fidelidade, para

cada passo no tempo Δ𝑡 até um instante 𝑡.

A.2 Sub-rotinas

• neighb_hilbert. Para um determinado estado, fornece os "vizinhos" no es-

paço de Hilbert considerando a interação efetiva entre todos os pares de q-bits

da rede.

• hamiltonian. Sub-rotina que coordena todas as outras relacionadas ao cálculo

dos elementos de matriz da Hamiltoniana.

• matrix_element. Calcula qualquer elemento de matriz da Hamiltoniana.

Apêndice A. Programas em FORTRAN 90 45

• xy_index. Atribui as coordenadas 𝑥 e 𝑦 de cada q-bit da rede.

• distance_periodic. Fornece a distância entre cada par de q-bit na rede

bidimensional com condições de contorno periódicas.

• Tnx. Calcula os polinômios de Chebyshev recursivamente dado um número

real qualquer.

• TnH. Calcula os polinômios de Chebyshev recursivamente para a Hamiltoni-

ana efetiva.

• gen_rand_vec. Gera os coeficientes 𝜉𝑟𝑘 da base truncada |𝑟⟩.

• dirichlet_kernel. Gera os coeficientes {𝑔𝑛} do kernel de Dirichlet.

• jackson_kernel. Gera os coeficientes {𝑔𝑛} do kernel de Jackson.

• lorentz_kernel. Gera os coeficientes {𝑔𝑛} do kernel de Lorentz.

• besselcoeff. Fornece as funções de Bessel de primeira espécie até determinada

ordem, dado um certo valor.


A.3 Códigos

A.3.1 mu_coeff

1 ! This program compute the mu(n) c o e f f i c i e n t s o f Chebyshev expansion2 ! We c o n s i d e r here a b id imens iona l l a t t i c e with p e r i o d i c boundary3 ! c o n d i t i o n s ( t o r i c code ) and long−range a n i s o t r o p i c i n t e r a c t i o n s4 ! between sp in p a i r s .5 !6 ! ND = H i l b e r t space dimension7 ! NR = number o f t runcated v e c t o r s8 ! NT = Chebyshev expansion order9 ! t2 = number o f ne ighbours in h i l b e r t space f o r each s t a t e vec to r

10 ! t = number o f qub i t s11 ! L = t o r i c code s i z e12 ! beta = i n t e r a c t i o n range13 ! Jx = i n t e r a c t i o n s t r enght f o r X operator14 ! Jz = i n t e r a c t i o n s t r enght f o r Z operator15 ! H = matrix element16 ! Hres = r e s c a l e d matrix element17 ! Emax = maximum e igenva lue ( est imated )18 ! Emin = minimum e igenva lue ( est imated )19 ! mu = Chebyshev expansion c o e f f i c i e n t s20 !21 ! ================22 ! Status : working OK23 ! Last update : 12/29/1124 ! Guilherme M. A. de Almeida25 ! ================26 !27 program mu_coeff28 !29 i n t ege r , parameter : : ND=2∗∗18 , NR=3, NT=10030 double p r e c i s i o n , parameter : : beta =1.d0 , Jz =1.5d0 , Jx=1.d031 double p r e c i s i o n , parameter : : p i = 3.1415926535898 d032 i n t ege r , parameter : : L=3, t =18, t2=t ∗( t −1)/233 i n t e g e r : : i , j , r , n , bra , temp34 double complex , dimension ( 0 :NR−1,ND) : : r v e c s e t35 double p r e c i s i o n : : E36 double p r e c i s i o n , dimension ( 0 :NT−1) : : mu37 double p r e c i s i o n , dimension ( 1 :NT−1, 0 :NR−1, 0 :ND−1) : : rvecaux38 i n t ege r , dimension ( t ) : : V, b39 double p r e c i s i o n : : H, Hres , Emax, Emin , Emax0 , Emin0 , f a t o r40 i n t ege r , dimension ( t2 +1) : : K41 i n t ege r , dimension ( t2 , t ) : : neighb , d iag42 i n t ege r , dimension ( t2 ) : : a43 double p r e c i s i o n , dimension ( t2 ) : : Haux44 double p r e c i s i o n , dimension ( 0 :ND−1) : : Hauxdiag45 !46 ! don ’ t worry about the se because they ’ re gonna be rep laced soon47 Emax=10.d048 Emin=−4.d049 Emax0=0.d050 Emin0=0.d051 !52 ! s e t t i n g Emax and Emin53 wr i t e (∗ ,∗ ) ’ S e t t i ng Emax and Emin ’54 !55 do i =0,ND−1


56 K(1)=i57 c a l l ne ighb_hi lber t (V, neighb , K, t , t2 , diag , a , b )58 bra=K(1)59 c a l l hami l tonian ( bra , K, H, Hres , beta , Jx , Jz , Emax, Emin , &60 L , t , t2 , neighb , diag , a )61 i f (H. gt . Emax0) Emax0=H62 i f (H. l t . Emin0) Emin0=H63 end do64 wr i t e (∗ ,∗ ) Emax0 , ’ and ’ ,Emin065 !66 ! s e t t i n g energy bounds ( gues s ing )67 open ( un i t =47, f i l e=’ bandwidth . txt ’ , s t a t u s=’unknown ’ )68 open ( un i t =48, f i l e=’ Emedio . txt ’ , s t a t u s=’unknown ’ )69 f a t o r =0.1d0∗Emax070 Emax=Emax0+f a t o r71 Emin=Emin0−f a t o r72 wr i t e (47 ,∗ ) Emax−Emin73 wr i t e (48 ,∗ ) (Emax+Emin) /2 . d074 c l o s e ( un i t =47)75 c l o s e ( un i t =48)76 !77 ! Generate | r> expansion c o e f f i c i e n t s78 c a l l gen_rand_vec ( rvec se t , ND, NR)79 !80 ! Creat ing o f f −d iagona l matrix e lements sequence81 K(1)=082 c a l l ne ighb_hi lber t (V, neighb , K, t , t2 , diag , a , b )83 do i =2, t2+184 bra=K( i )85 c a l l hami l tonian ( bra , K, H, Hres , beta , Jx , Jz , Emax, Emin , &86 L , t , t2 , neighb , diag , a )87 Haux( i −1)=Hres88 end do89 !90 ! Creat ing d iagona l matrix e lements sequence91 do i =0,ND−192 K(1)=i93 c a l l ne ighb_hi lber t (V, neighb , K, t , t2 , diag , a , b )94 bra=K(1)95 c a l l hami l tonian ( bra , K, H, Hres , beta , Jx , Jz , Emax, Emin , &96 L , t , t2 , neighb , diag , a )97 !98 Hauxdiag ( i )=Hres99 end do

100 !101 ! S e t t i ng | r>_n = 0102 do i =0,ND−1103 do r =0,NR−1104 do n=1,NT−1105 rvecaux (n , r , i ) =0.d0106 end do107 end do108 end do109 !110 ! the ac t i on s t a r t s r i g h t here111 wr i t e (∗ ,∗ ) ’ 1 ’112 do i =1,ND113 K(1)=i −1114 c a l l ne ighb_hi lber t (V, neighb , K, t , t2 , diag , a , b )


115 do r =0,NR−1116 do j =1, t2+1117 i f (K( j ) . ne .K(1) ) then118 rvecaux (1 , r ,K( j ) )=rvecaux (1 , r ,K( j ) )+ &119 s q r t ( 2 . d0 ) ∗ r e a l ( r v e c s e t ( r , i ) ) ∗Haux( j −1)120 e l s e121 rvecaux (1 , r ,K( j ) )=rvecaux (1 , r ,K( j ) )+ &122 s q r t ( 2 . d0 ) ∗ r e a l ( r v e c s e t ( r , i ) ) ∗Hauxdiag (K(1) )123 end i f124 end do125 end do126 end do127 !128 ! second step129 wr i t e (∗ ,∗ ) ’ 2 ’130 do i =1,ND131 K(1)=i −1132 c a l l ne ighb_hi lber t (V, neighb , K, t , t2 , diag , a , b )133 do r =0,NR−1134 do j =1, t2+1135 i f (K( j ) . ne .K(1) ) then136 rvecaux (2 , r ,K( j ) )=rvecaux (2 , r ,K( j ) )+ &137 2 . d0∗ rvecaux (1 , r , i −1)∗Haux( j −1)138 e l s e139 rvecaux (2 , r ,K( j ) )=rvecaux (2 , r ,K( j ) )+ &140 2 . d0∗ rvecaux (1 , r , i −1)∗Hauxdiag ( i −1)141 end i f142 end do143 end do144 end do145 !146 do i =1,ND147 do r =0,NR−1148 rvecaux (2 , r , i −1)=rvecaux (2 , r , i −1)− &149 r e a l ( r v e c s e t ( r , i ) ) ∗ s q r t ( 2 . d0 )150 end do151 end do152 !153 ! u n t i l the end154 do n=3,NT−1155 wr i t e (∗ ,∗ ) n156 do i =1,ND157 K(1)=i −1158 c a l l ne ighb_hi lber t (V, neighb , K, t , t2 , diag , a , b )159 do r =0,NR−1160 do j =1, t2+1161 i f (K( j ) . ne .K(1) ) then162 rvecaux (n , r ,K( j ) )=rvecaux (n , r ,K( j ) )+ &163 2 . d0∗ rvecaux (n−1, r , i −1)∗Haux( j −1)164 e l s e165 rvecaux (n , r ,K( j ) )=rvecaux (n , r ,K( j ) )+ &166 2 . d0∗ rvecaux (n−1, r , i −1)∗Hauxdiag ( i −1)167 end i f168

169 end do170 end do171 end do172 !173 do i =1,ND


174 do r =0,NR−1175 rvecaux (n , r , i −1)=rvecaux (n , r , i −1)−rvecaux (n−2, r , i −1)176 end do177 end do178 end do179 !180 ! computing the mu(n) c o e f f i c i e n t s181 open ( un i t =11, f i l e=’mu. txt ’ , s t a t u s=’unknown ’ )182 do n=0,NT−1183 mu(n) =0.d0184 end do185 do r =0,NR−1186 do i =1,ND187 mu(0) =2.d0 ∗ ( 1 . 0 d0/NR) ∗ r e a l ( r v e c s e t ( r , i ) ) ∗ r e a l ( r v e c s e t ( r , i ) )+mu(0)188 mu(1) =(1.0 d0/NR) ∗ r e a l ( r v e c s e t ( r , i ) ) ∗ rvecaux (1 , r , i −1)∗ s q r t ( 2 . d0 )+mu(1)189 do n=2,NT−1190 mu(n) =(1.0 d0/NR) ∗ r e a l ( r v e c s e t ( r , i ) ) ∗ rvecaux (n , r , i −1)∗ s q r t ( 2 . d0 )+mu(n)191 end do192 end do193 end do194 !195 ! sav ing them196 do n=0,NT−1197 wr i t e (11 ,∗ ) mu(n)198 end do199 !200 c l o s e ( un i t =11)201 !202 end program203 !204 ! − ∗ Inc lud ing subrout ine s205 i n c lude " ne ighb_hi lber t . f 90 "206 i n c lude " hami l tonian . f90 "207 i n c lude " gen_ran_vec . f90 "


A.3.2 dos_cheb

1 ! Program to eva luate the dens i ty o f s t a t e s (DOS) .2 ! I t depends on a mu data f i l e c r ea ted by the3 ! program mu_n.4 !5 ! ================6 ! Status : working OK7 ! Last update : 12/29/118 ! Guilherme M. A. de Almeida9 ! ================

10 !11 program dos_cheb12 i n t ege r , parameter : : ND=2∗∗18 , NT=10013 i n t e g e r : : i , j , n14 double p r e c i s i o n , dimension ( 0 :NT−1) : : t t15 double p r e c i s i o n : : dens i ty , E, deltaE , Emed16 double p r e c i s i o n , parameter : : p i = 3.1415926535898 d017 double p r e c i s i o n , dimension ( 0 :NT−1) : : mu, g18 double p r e c i s i o n , parameter : : lambda=3.d019 !20 open ( un i t =11, f i l e=’mu. txt ’ , s t a t u s=’unknown ’ , a c t i on=’ read ’ )21 open ( un i t =12, f i l e=’ dosdata . txt ’ , s t a t u s=’unknown ’ )22 open ( un i t =54, f i l e=’ bandwidth . txt ’ , s t a t u s=’unknown ’ , a c t i on=’ read ’ )23 open ( un i t =55, f i l e=’ Emedio . txt ’ , s t a t u s=’unknown ’ , a c t i on=’ read ’ )24 !25 do n=0, NT−126 read (11 ,∗ ) mu(n)27 end do28 !29 read (54 ,∗ ) deltaE30 read (55 ,∗ ) Emed31 c l o s e ( un i t =54)32 c l o s e ( un i t =55)33 !34 ! s e l e c t the ke r ne l you want35 c a l l d i r i c h l e t _ k e r n e l ( NT, g )36 ! c a l l jackson_kerne l ( NT, g , p i )37 ! c a l l l o r en t z_ke rne l ( NT, lambda , g )38 !39 ! Loop over energy40 do i =−99 ,99 ,141 E=0.01d0 ∗( i )42 dens i ty = 0 . d043 do n=1,NT−144 c a l l Tnx(E, NT, t t )45 dens i ty = g (n) ∗2 . d0∗mu(n) ∗ t t (n) + dens i ty46 end do47 dens i ty= dens i ty /( p i ∗ s q r t ( 1 . d0−E∗E) )+(mu(0) ∗g (0 ) ) /( p i ∗ s q r t ( 1 . d0−E∗E) )48 wr i t e (12 ,∗ ) dens i ty49 end do50 !51 ! S e t t i ng X−a x i s52 open ( un i t =55, f i l e=’ ene rgyax i s . txt ’ , s t a t u s=’unknown ’ )53 do i =−99 ,99 ,154 E=0.01d0∗ i ∗ deltaE /2 . d0+Emed55 wr i t e (55 ,∗ ) E56 end do57 c l o s e ( un i t =55)


58 !59 c l o s e ( un i t =12)60 c l o s e ( un i t =11)61 !62 end program63 !64 ! − ∗ Inc lud ing subrout ine s65 i n c lude "Tnx . f90 "66 i n c lude " d i r i c h l e t _ k e r n e l . f 90 "67 i n c lude " jackson_kerne l . f 90 "68 i n c lude " l o r en t z_kerne l . f 90 "


A.3.3 fidelity_cheb

1 ! Program to compute the time e v o l u t i o n o f a s t a t e vec to r and so the f i d e l i t y2 ! a f t e r each time step with r e s p e c t to the chosen i n i t i a l s t a t e .3 ! Here we c o n s i d e r a L = 3 t o r i c code .4 !5 ! ================6 ! Status : working OK7 ! Last update : 1/10/128 ! Guilherme M. A. de Almeida9 ! ================

10 !11 program f i d e l i t y _ c h e b12 !13 i n t ege r , parameter : : ND=2∗∗18 , NT=1514 r ea l , parameter : : Jx =1.0 , Jz =0.5 , beta =1.015

16 r ea l , parameter : : p i = 3.141592617 i n t ege r , parameter : : L=3, t =18, t2=t ∗( t −1)/218 !19 r ea l , parameter : : t o t a l t i m e =60.020 i n t ege r , parameter : : s t e p s =2021 !22 i n t e g e r : : i , j , r , n , bra , temp , nothing23 r ea l , dimension ( 0 :NT−1) : : g24 i n t ege r , dimension ( t ) : : V, b25 r e a l : : H, Hres , Emax, Emin , Emax0 , Emin0 , f a t o r26 i n t ege r , dimension ( t2 +1) : : K27 i n t ege r , dimension ( t2 , t ) : : neighb , d iag28 i n t ege r , dimension ( t2 ) : : a29 r ea l , dimension ( t2 ) : : Haux30 r ea l , dimension ( 0 :ND−1) : : Hauxdiag31 r e a l : : time32 r ea l , dimension ( 0 :NT−1) : : J b e s s e l33 complex , dimension ( 0 :NT−1) : : gama34 complex , dimension ( 0 :NT−1, 0 :ND−1) : : p s i35 complex , dimension ( 0 :ND−1) : : p s i 036 complex : : Fid , imfac37 r e a l : : f i d e l i t y38 r e a l : : norm , normcheck39 r ea l , parameter : : lambda=3.040 !41 ! s e t t i n g dt=time and c a l l i n g the k e r n e l and Bes s e l f u n c t i o n s42 time=t o t a l t i m e / s t e p s43 c a l l b e s s e l c o e f f ( NT, time , J b e s s e l )44 !45 ! s e l e c t the ke r ne l you want46 c a l l d i r i c h l e t _ k e r n e l ( NT, g )47 ! c a l l jackson_kerne l ( NT, g , p i )48 ! c a l l l o r en t z_ke rne l ( NT, lambda , g )49 !50 ! t h i s i s only to check which order o f NT may be used51 wr i t e (∗ ,∗ ) ( J b e s s e l (n) ∗g (n) , n=0,NT−1)52 wr i t e (∗ ,∗ ) ’___’53 wr i t e (∗ ,∗ ) ( g (n) , n=0,NT−1)54 !55 ! don ’ t worry about the se because they ’ re gonna be rep laced56 Emax=16.057 Emin=−4.0


58 Emax0=0.059 Emin0=0.060 !61 ! s e t t i n g Emax and Emin62 wr i t e (∗ ,∗ ) ’ S e t t i ng Emax and Emin ’63 !64 do i =0,ND−165 K(1)=i66 c a l l ne ighb_hi lber t (V, neighb , K, t , t2 , diag , a , b )67 bra=K(1)68 c a l l hami l tonian ( bra , K, H, Hres , beta , Jx , Jz , Emax, Emin , &69 L , t , t2 , neighb , diag , a )70 i f (H. gt . Emax0) Emax0=H71 i f (H. l t . Emin0) Emin0=H72 end do73 wr i t e (∗ ,∗ ) Emax0 , ’ and ’ ,Emin074 !75 ! s e t t i n g energy bounds ( gues s ing )76 open ( un i t =47, f i l e=’ bandwidth . txt ’ , s t a t u s=’unknown ’ )77 f a t o r=Emax0∗0 .178 Emax=Emax0+f a t o r79 Emin=Emin0−f a t o r80 wr i t e (47 ,∗ ) Emax−Emin81 c l o s e ( un i t =47)82 !83 ! c r e a t i n g o f f −d iagona l matrix e lements sequence84 wr i t e (∗ ,∗ ) ’ Creat ing o f f −d iagona l matrix e l emetents sequence ’85 K(1)=086 !87 ! f i n d i n g ne ighbours in H i l b e r t space88 c a l l ne ighb_hi lber t (V, neighb , K, t , t2 , diag , a , b )89 do i =2, t2+190 bra=K( i )91 !92 ! computing matrix e lements93 c a l l hami l tonian ( bra , K, H, Hres , beta , Jx , Jz , Emax, Emin , &94 L , t , t2 , neighb , diag , a )95 !96 ! s e t t i n g the hami l tonian matrix e lements sequence to the a n c i l l a r y vec to r97 Haux( i −1)=Hres98 wr i t e (∗ ,∗ ) H99 end do

100 !101 ! c r e a t i n g d iagona l matrix e lements sequence102 open ( un i t =13, f i l e=’ checkdiag . txt ’ , s t a t u s=’unknown ’ )103 wr i t e (∗ ,∗ ) ’ Creat ing d iagona l matrix e l emetents sequence ’104 !105 do i =0,ND−1106 K(1)=i107 c a l l ne ighb_hi lber t (V, neighb , K, t , t2 , diag , a , b )108 bra=K(1)109 c a l l hami l tonian ( bra , K, H, Hres , beta , Jx , Jz , Emax, Emin , &110 L , t , t2 , neighb , diag , a )111 Hauxdiag ( i )=Hres112 wr i t e (13 ,∗ ) H113 end do114 c l o s e ( un i t =13)115 !116 ! s e t t i n g the f i x e d i n i t i a l s t a t e do |0>


117 do i =0,ND−1118 ps i 0 ( i ) =(0 . 0 , 0 . 0 )119 end do120 ps i 0 (0 ) =(1 . 0 , 0 . 0 )121 !122 ! normal i z ing123 norm=0.0124 do i =0,ND−1125 norm=norm+conjg ( p s i 0 ( i ) ) ∗ ps i 0 ( i )126 end do127 norm=1.0/norm128 norm=s q r t (norm)129 do i =0,ND−1130 ps i 0 ( i )=ps i 0 ( i ) ∗norm131 end do132 !133 ! s e t t i n g the f i r s t r e c u r s i v e i n i t i a l s t a t e134 do i =0,ND−1135 p s i (0 , i )=ps i 0 ( i )136 end do137 !138 ! loop over time s t e p s139 open ( un i t =12, f i l e=’ f i d e l i t y . txt ’ , s t a t u s=’unknown ’ )140 open ( un i t =14, f i l e=’ norma . txt ’ , s t a t u s=’unknown ’ )141 !142 do j =1, s t e p s143 wr i t e (∗ ,∗ ) ’ running step ’ , j144 !145 ! s e t t i n g a l l vec to r c o e f f i c i e n t s equal to zero146 do n=1,NT−1147 do i =0,ND−1148 p s i (n , i ) =(0 . 0 , 0 . 0 )149 end do150 end do151 !152 ! computing the Chebyshev p o l i n o m i a l s153 c a l l TnH( ps i0 , ps i , Haux , &154 Hauxdiag , Emax, Emin , gama , NT, ND, t , t2 )155 !156 ! s e t t i n g a new i n i t i a l s t a t e ( t h i s i s the time e v o l u t i o n )157 do i =0,ND−1158 p s i (0 , i )=g (0 ) ∗ J b e s s e l (0 ) ∗ p s i (0 , i )159 imfac =(1 . 0 , 0 . 0 )160 do n=1,NT−1161 imfac = imfac ∗ (0 .0 , −1 .0)162 p s i (0 , i )=p s i (0 , i ) +2.0∗g (n) ∗ J b e s s e l (n) ∗ p s i (n , i ) ∗ imfac163 end do164 ! wr i t e (50 ,∗ ) p s i (0 , i )165 end do166 !167 ! Checking i f the s t a t e i s normal ized168 normcheck =0.0169 do i =0,ND−1170 normcheck=normcheck+conjg ( p s i (0 , i ) ) ∗ p s i (0 , i )171 end do172 wr i t e (∗ ,∗ ) ’ norma=’ , normcheck173 wr i t e (14 ,∗ ) normcheck174 !175 ! computing the f i d e l i t y


176 f i d =(0 . 0 , 0 . 0 )177 do i =0,ND−1178 f i d=f i d+conjg ( p s i 0 ( i ) ) ∗ p s i (0 , i )179 end do180 f i d e l i t y=abs ( f i d )181 wr i t e (∗ ,∗ ) ’ F ide l i dade ’ , j , ’= ’ , f i d e l i t y182 wr i t e (12 ,∗ ) f i d e l i t y183 end do ! j184 !185 c l o s e ( un i t =12) ! f i d e l i t y186 c l o s e ( un i t =14) ! norm187 !188 stop189 end program190 !191 ! − ∗ Inc lud ing subrout ine s192 i n c lude " b e s s e l c o e f f . f 90 "193 i n c lude " jackson_kerne l . f 90 "194 i n c lude " l o r en t z_kerne l . f 90 "195 i n c lude " d i r i c h l e t _ k e r n e l . f 90 "196 i n c lude " ne ighb_hi lber t . f 90 "197 i n c lude " hami l tonian . f90 "198 i n c lude "TnH. f90 "


A.3.4 TnH

1 ! Subrout ine to compute the Chebyshev polynomia l s T_n(H) r e c u r s i v e l y .2 !3 ! ================4 ! Status : working OK5 ! Last update : 1/10/126 ! Guilherme M. A. de Almeida7 ! ================8 !9 subrout ine TnH( ps i0 , ps i , Haux , &

10 Hauxdiag , Emax, Emin , gama , NT, ND, t , t2 )11 !12 i n t e g e r : : i , j , m, n , t , t2 , ND, NT13 i n t ege r , dimension ( t ) : : V, b14 i n t ege r , dimension ( t2 , t ) : : neighb , d iag15 i n t ege r , dimension ( t2 +1) : : K16 i n t ege r , dimension ( t −1) : : f17 i n t ege r , dimension ( t2 ) : : a18 r e a l : : Emax, Emin19 r ea l , dimension ( t2 ) : : Haux20 r ea l , dimension ( 0 :ND−1) : : Hauxdiag21 complex , dimension ( 0 :NT−1) : : gama22 complex , dimension ( 0 :NT−1, 0 :ND−1) : : p s i23 complex , dimension ( 0 :ND−1) : : p s i 024 r e a l : : norm , normtest25 !26 ! s t a r t i n g with n=127 wr i t e (∗ ,∗ ) ’ 1 ’28 !29 do i =0,ND−130 K(1)=i31 c a l l ne ighb_hi lber t (V, neighb , K, t , t2 , diag , a , b )32 do j =2, t2+133 p s i (1 ,K( j ) )=p s i (1 ,K( j ) )+Haux( j −1)∗ p s i (0 , i )34 end do35 p s i (1 ,K(1) )=p s i (1 ,K(1) )+Hauxdiag (K(1) ) ∗ p s i (0 , i )36 end do37 wr i t e (∗ ,∗ ) p s i ( 1 , 0 )38 !39 do n=2,NT−140 wr i t e (∗ ,∗ ) n41 do i =0,ND−142 K(1)=i43 c a l l ne ighb_hi lber t (V, neighb , K, t , t2 , diag , a , b )44 do j =2, t2+145 p s i (n ,K( j ) )=p s i (n ,K( j ) ) +2.0∗ p s i (n−1, i ) ∗Haux( j −1)46 end do47 p s i (n ,K(1) )=p s i (n ,K(1) ) +(2.0∗ p s i (n−1, i ) ) ∗Hauxdiag (K(1) )48 end do49 do i =0,ND−150 p s i (n , i )=p s i (n , i )−p s i (n−2, i )51 end do52 !53 end do54 !55 re turn56 end subrout ine

Documents

UNIVERSIDADE FEDERAL DE SERGIPE Núcleo de Pós-Graduação … · Núcleo de Pós-Graduação em Física Método Kernel Polinomial aplicado a uma rede de spins em ambiente correlacionado