UNIVERSIDADE ESTADUAL PAULISTA “JÚLIO DE MESQUITA FILHO”
FACULDADE DE CIÊNCIAS/CAMPUS DE BAURU

PROGRAMA DE PÓS-GRADUAÇÃO EM CIÊNCIA E TECNOLOGIA DE
MATERIAIS

Estudo de Materiais Supercondutores em Forma de
SQUID

com uma Constrição Usando Métodos de Paralelização
Computacional

André Luiz Severino

Orientador:

Edson Sardella

BAURU
2015



André Luiz Severino

Estudo de Materiais Supercondutores em Forma de SQUID
com uma Constrição Usando Métodos de Paralelização Computacional

Dissertação apresentada ao programa de
Pós-Graduação em Ciência e Tecnologia

de Materiais da Universidade Estadual
Paulista como parte dos pré-requisitos

para a obtenção do título de Mestre.

BAURU
2015



 

 

 

 

 

 

 

 

 

 

 

   

 

  
 

 

Severino, André Luiz. 

   Estudo de materiais supercondutores em forma de 

SQUID com uma constrição usando métodos de 

paralelização computacional / André Luiz Severino, 

2015 

   100 f.   

 

   Orientador: Edson Sardella 

    

   Dissertação (Mestrado)– Universidade Estadual 

Paulista. Faculdade de Ciências, Bauru, 2015 

 

1. Supercondutividade. 2. Vórtices. 3. Equações de 
London. 4. Equações de Ginzburg-Landau. 5. Método psi-

U. 6. Squid. 7. Paralelização. 8. CUDA. 9. Linguagem 

C. 10. Mesoscópico. I. Universidade Estadual Paulista. 

Faculdade de Ciências.  

 

    

 

  

 

 

 



 



Dedicatória

Dedico este trabalho à pessoa que encheu minha vida de amor: Sandra Alves
Pereira. Obrigado por apoiar-me.



Agradecimentos
Agradeço à toda minha família, e em especial, ao meu avô Agenor, que sem-

pre compartilhou muitos momentos de sabedoria comigo.

Agradeço ao meu orientador, Professor Edson Sardella, que ensinou-me, além
de Física e tantas outras coisas, também a ser sempre alegre, mesmo nas
dificuldades. Tive confiança nele, e ele confiou em mim. Confiança recí-
proca.

Agradeço aos meus amigos de classe que, além de enfrentarem as dificul-
dades e desafios das disciplinas comigo, também alegraram minha vida.
Tempos memoráveis.

Agradeço aos meus amigos de Jaú, minha cidade natal, e de Bauru. Cada um
sabe o quanto os considero.

Agradeço aos meus professores da pós-graduação, Alexandre, Alexys, Pa-
blo, Paulo e Vicente, pelo conhecimento que adquiri graças à excelente
didática de todos.

Agradeço aos meus colegas de trabalho da EPC de Jaú pelas risadas e ami-
zade.

Agradeço à Natureza por formar e desenvolver minha consciência - mesmo
que de maneira possivelmente aleatória.



Resumo

É notável o desenvolvimento atual da Ciência e Tecnologia relacionada aos materiais

supercondutores. Exemplos de utilização deste tipo de material, são o armazenamento de

energia em forma de supercorrentes, trens do tipo Maglev, cujo uso das supercorrentes

também proporciona a geração de altos valores de campo magnético que podem propor-

cionar levitação e propulsão. Aplicações desta natureza, dentre outras, requerem uma

profunda compreensão do comportamento das propriedades físicas fundamentais dos su-

percondutores, tanto do ponto de vista macroscópico quanto microscópico. Como exem-

plo, podemos citar a formação da rede de vórtices em materiais supercondutores na pre-

sença de campos magnéticos aplicados e/ou de correntes de transporte, cujas magnitu-

des, quando ultrapassam determinados valores críticos, permitem a entrada quantizada de

fluxo magnético.

Com relação aos métodos utilizados neste trabalho, estudamos o comportamento e

as características de um material de dimensões mesoscópicas em formato de dispositivo

supercondutor de interferência quântica - "SQUID", um tipo de magnetômetro muito sen-

sível, capaz de realizar medidas de campos magnéticos muito tênues - usando um algo-

ritmo de simulação computacional. A simulação basicamente atua no sentido de se aplicar

campos magnéticos em passos de tempo determinados - a temperatura fixa - ao material

modelado, o que permite o estudo das características deste tipo de material por meio dos

dados de saída gerados pelo simulador, os quais permitem a construção de gráficos para

estudo do comportamento do parâmetro de ordem, densidades de corrente, formação de

vórtices, entre outros parâmetros relacionados. Assim, fazemos uma revisão das teorias

fenomenológicas e desenvolvemos um algorítimo de solução numérica das equações en-

volvidas, as quais não possuem solução exata. Em seguida usamos a tecnologia CUDA



com a linguagem C para implementar este algoritmo e aplicamos à solução do problema

de um SQUID com uma pequena constrição. Por fim, os gráficos são construídos utili-

zando os dados gerados pelo simulador para análise e interpretação do comportamento do

material em estudo.

Palavras chave: Supercondutividade, vórtices, equações de London, equações de

Ginzburg-Landau, método ψ − U , SQUID, paralelização, CUDA, linguagem C, mesos-

cópico.



Abstract

It is remarkable the current development of science and technology related to super-

conducting materials. Examples of using this type of material, such as storage of energy

as supercurrents, the Maglev-type trains, which also provides use of supercurrents ge-

nerating high values of magnetic field that can provide levitation and propulsion. Ap-

plications of this nature, among others, require a deep understanding of the behavior of

the fundamental physical properties of superconductors, both from the macroscopic and

microscopic point of view. As one example, the formation of vortex lattice in super-

conducting materials in the presence of applied magnetic fields and/or transport chains

whose magnitudes when they exceed certain critical values allow the quantized magnetic

flux entrance.

With respect to the methods used in this work, we studied the behavior and charac-

teristics of a material with mesoscopic dimensions in a geometry of a superconducting

quantum interference device - "SQUID". Such device is a kind of very sensitive magneto-

meter, able to perform very faint magnetic fields measures - using a computer simulation

algorithm. The simulation basically acts to apply magnetic fields at certain time steps -

with fixed temperature - to the modeled material, allowing the study of the characteristics

of this type of material by means of output data generated by the simulator, which allows

the construction of graphs to study the order parameter behavior, current densities, vortex

formation, and other related parameters. Thus, we review the phenomenological theories

and developed a numerical solution algorithm of the involved equations, which have no

exact solution. Then we use the CUDA technology with the C language to implement

this algorithm and apply it to solving the problem of a SQUID with a small constriction.

Finally, the graphics are constructed using the data generated by the simulator for analysis



and interpretation of the behavior of the material under study.

Key words: Superconductivity, vortex, London equations, Ginzburg-Landau equati-

ons, ψ − U method, SQUID, parallelization, CUDA, C language, mesoscopic.



Lista de Figuras

1.1 Gráfico histórico de Onnes: resistência do mercúrio sólido como função

da temperatura. [1] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

1.2 Experimentos ilustrando a diferença entre um supercondutor e um con-

dutor perfeito. Experimento 1: amostras resfriadas sem a presença de

um campo magnético aplicado. Experimento 2: amostras resfriadas na

presença de um campo magnético aplicado. Fonte: Próprio autor. . . . . . 19

1.3 (a) Gráfico geral da dependência do campo crítico Hc(T ) pela tempera-

tura para diversos materiais. [1]. (b) Curva do calor específico em função

da temperatura para um material supercondutor na ausência de campo

magnético, onde verifica-se uma descontinuidade. [2] (Adaptada da fonte

original). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

2.1 (a) Perfil do campo magnético local em um supercondutor semi-infinito

na presença de uma campo aplicado. Note que, na fronteira, Bz(0) = H .

(b) Superfície curva demonstrando a penetração de campo na direção do

eixo x. Fonte: Próprio autor. . . . . . . . . . . . . . . . . . . . . . . . . 25

7



2.2 Perfis da energia livre. Em (a), (b) e (c), esses perfis não satisfazem a

teoria de Landau de transições de fase de segunda ordem, enquanto em

(d), a teoria é satisfeita, com mínimo ocorrendo em M 6= 0, indicado pela

seta apontada sobre a curva de cor verde. Fonte: Próprio autor. . . . . . . 28

2.3 Ilustração do comportamento do parâmetro de ordem próximo à interface

metal-supercondutor. Após o comprimento ξ, o material encontra-se no

estado supercondutor com valor máximo de parâmetro de ordem. . . . . . 33

2.4 Representação esquemática de perfis de campo e de parâmetro de ordem

para supercondutores do tipo-I e II, onde a configuração em (b) permite

a formação de um vórtice, enquanto a configuração em (a) não. Fonte:

Próprio autor. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35

2.5 Diagrama de fase de um supercondutor do tipo-I. Fonte: Próprio autor. . . 36

2.6 Diagrama de fase de um supercondutor do tipo-II. Fonte: Próprio autor. . 37

2.7 (a) Magnetização de um supercondutor do tipo-I e (b) Magnetização de

um supercondutor do tipo-II, ambos em função do campo aplicado H . A

linha tracejada representa um supercondutor do tipo-I. Fonte: Dissertação

de Mestrado de Mauro César Vieira Pascolati. Unesp Bauru/POSMAT -

Ano 2010. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38

2.8 Estrutura de um vórtice. Fonte: Dissertação de Mestrado de Mauro César

Vieira Pascolati. Unesp Bauru/POSMAT - Ano 2010. . . . . . . . . . . . 40

2.9 Rede de vórtices de Abrikosov: perfil do campo local. Experimento com

material supercondutor, onde partículas ferromagnéticas foram aplicadas

à superfície deste material, permanecendo nesta configuração. As regiões

mais escuras correspondem aos valores mais altos de campo (região nor-

mal), enquanto as regiões mais claras, aos valores mais baixos de campo

(região supercondutora)[1]. . . . . . . . . . . . . . . . . . . . . . . . . . 42



3.1 Representação da caixa de simulação com o supercondutor no seu interior.

Fonte: Próprio autor. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50

3.2 Malha de discretização: a legenda indica os pontos onde cada quantidade

física é calculada. Fonte: Próprio autor. . . . . . . . . . . . . . . . . . . 53

3.3 A figura ilustra a aproximação (3.35) para a primeira derivada [3]. . . . . 55

4.1 Hardware básico para computação paralela com CUDA: Com a imple-

mentação da placa de vídeo (device), o programa, antes escrito para ser

executado apenas no host ((a) Memória principal e (b) CPU), agora tam-

bém pode ser executado no device ((c) Memória e (d) GPU). Fonte: Pró-

prio autor. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65

4.2 Estrutura básica de processamento em paralelo. Fonte: Próprio autor. . . . 66

4.3 Esquema simplificado de alocação dinâmica de memória em C/C++. Os

cubos em (a), (b) e (c) representam regiões de memória que armazenam

o endereço de outras regiões de memória (por isso são denominados de

ponteiros). Em (d) encontra-se o espaço alocado para ser utilizado pelos

dados processados. Fonte: Próprio autor. . . . . . . . . . . . . . . . . . . 68

4.4 Achatamento de índice. Em (a), temos o vetor tridimensional antigo, e

em (b), o novo vetor unidimensional, ambos com quantidades idênticas de

memória alocada. Após esse processo, o acesso às posições do novo vetor

em (b) - o qual está relacionado com as posições do vetor tridimensional

em (a) - é feito por meio de um índice unidimensional. Fonte: Próprio autor. 69

4.5 A linha 1 cria um ponteiro para memória do tipo específico que se quer

trabalhar, enquanto na linha 2, além da quantidade e tamanho desses es-

paços de memória, o comando retorna um ponteiro da memória alocada

no device para vetor3D. Fonte: Próprio autor. . . . . . . . . . . . . . . . 70



4.6 (a) Transferência HostToDevice. (b) Transferência DeviceToHost. Fonte:

Próprio autor. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70

4.7 Código de execução serial. Fonte: Próprio autor. . . . . . . . . . . . . . . 71

4.8 Dimensionamento em 3D do processo paralelizado. Fonte: Próprio autor. 73

4.9 Código de execução paralela para o antigo laço “for”. Fonte: próprio autor. 74

4.10 Ilustração da geometria do supercondutor em forma de geometria SQUID

com uma constrição na lateral direita da fenda. Fonte: Próprio autor. . . . 75

4.11 Função básica para cálculo de exponencial de número do tipo cuDouble-

Complex. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77

4.12 Exemplo de alocação de memória: o vetor F, utilizado para guardar o

valor do parâmetro de ordem, e a variável VX, variável de ligação, são

alocadas no host, enquanto dev_F e dev_VX, no device. . . . . . . . . . . 77

4.13 Cópia dos dados do vetor localizado no host para o respectivo vetor no

device. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77

4.14 Declaração de variáveis do tipo dim3. . . . . . . . . . . . . . . . . . . . 78

4.15 Chamada de kernel. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78

4.16 Corpo do kernel referido. . . . . . . . . . . . . . . . . . . . . . . . . . . 78

5.1 Magnetização versus campo aplicado. . . . . . . . . . . . . . . . . . . . 83

5.2 Energia versus campo aplicado. . . . . . . . . . . . . . . . . . . . . . . . 84

5.3 (a)-(c) Intensidade ψ do parâmetro de ordem; (d) Fase φ do parâmetro de

ordem; (e) A paleta de cores ilustra a intensidade do parâmetro de ordem

variando de 0 a 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86

5.4 (a)-(b) Perfil do campo local para os campos aplicados indicados nos pai-

néis; (c)-(d) Respectivas densidades de correntes. . . . . . . . . . . . . . 87

5.5 Intensidade do parâmetro de ordem para campos intermediários (painéis

(a)-(c)) e para campo alto (painel (d)). . . . . . . . . . . . . . . . . . . . 88



Lista de Tabelas

4.1 Tempos de execução do código para alguns valores de campo aplicado H . 79

11



Sumário

1 Introdução 14
1.1 Revisão Histórica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
1.2 O Fenômeno da Supercondutividade . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

2 Teorias Fenomenológicas da Supercondutividade 21
2.1 Teoria Fenomenológica de London . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

2.1.1 Introdução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.1.2 Primeira Abordagem da Teoria de London . . . . . . . . . . . . . . . . . . . . . . 22
2.1.3 Segunda Abordagem da Teoria de London . . . . . . . . . . . . . . . . . . . . . . 24
2.1.4 Efeito Meissner . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

2.2 Teoria Fenomenológica de Ginzburg-Landau . . . . . . . . . . . . . . . . . . . . . . . . 26
2.3 Comprimentos Característicos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

2.3.1 Comprimento de Coerência . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
2.3.2 Comprimento de Penetração de London . . . . . . . . . . . . . . . . . . . . . . . 33
2.3.3 Parâmetro de Ginzburg-Landau . . . . . . . . . . . . . . . . . . . . . . . . . . . 34

2.4 Classificação dos Supercondutores . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
2.5 Estrutura de Vórtices em Supercondutores do Tipo-II . . . . . . . . . . . . . . . . . . . . 38

2.5.1 Quantização do Fluxo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
2.5.2 Vórtice Isolado . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
2.5.3 Interação entre Vórtices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41

2.6 Supercondutividade Mesoscópica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

3 O Método ψ − U 43
3.1 Introdução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
3.2 Equações TDGL . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
3.3 Transformações de Calibre . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
3.4 Estado Meissner . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
3.5 Unidades Reduzidas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
3.6 Campos Auxiliares . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
3.7 Problema a Ser Resolvido: Condições de Contorno . . . . . . . . . . . . . . . . . . . . . 49
3.8 Discretização das Equações TDGL . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51

3.8.1 Malha de Discretização . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
3.8.2 Definições . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52

12



3.8.3 Aproximações para as Derivadas . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
3.8.4 Discretização da Primeira Equação TDGL . . . . . . . . . . . . . . . . . . . . . . 55
3.8.5 Discretização da Segunda Equação TDGL (Lei de Ampère) . . . . . . . . . . . . 56
3.8.6 Densidade de Corrente . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
3.8.7 Campo Local . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
3.8.8 Condições de Contorno . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59

3.9 Evolução Temporal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
3.10 Algoritmo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61

4 Métodos de Paralelização do Algoritmo ψ − U 63
4.1 Conceitos Básicos Sobre Paralelização com Tecnologia CUDA R© e Linguagem C . . . . . 63

4.1.1 Introdução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
4.1.2 Threads, Blocks e Grids . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71

4.2 Implementação da Paralelização no Método ψ − U . . . . . . . . . . . . . . . . . . . . . 75
4.2.1 Características Gerais do Sistema Simulado . . . . . . . . . . . . . . . . . . . . . 75
4.2.2 Implementação da Paralelização do Método ψ − U . . . . . . . . . . . . . . . . . 76

5 Supercondutor em Forma de um SQUID com uma Constrição 80
5.1 Introdução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
5.2 Parâmetros Utilizados e Metodologia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
5.3 Resultados e Discussão . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82

5.3.1 Magnetização e Energia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
5.3.2 Configurações da Rede de Vórtices . . . . . . . . . . . . . . . . . . . . . . . . . 84

6 Conclusões 89

7 Apêndices 91

A Cálculo dos passos de tempo 92

B Modelagem do Simulador 94



1
Introdução

1.1 Revisão Histórica

A descoberta do fenômeno da supercondutividade ocorreu quando o físico holandês

Heike Kamerlingh Onnes [4] verificou, na data de 8 de abril de 1911, o que ele cha-

mou posteriormente de supercondutividade, um fenômeno até então desconhecido aos

olhos dos cientistas da época. Após conseguir liquefazer o hélio pela primeira vez em

1908, ele e seu colaborador Gilles Holst observaram três anos depois, que o mercúrio,

ao ser submetido a uma temperatura abaixo de 4, 2 K por essa técnica de resfriamento,

apresentava efeito Joule nulo, ou seja, a resistência do material à passagem de corrente

elétrica era aproximadamente zero, e a corrente podia fluir pelo material sem apresentar

nenhuma perda. Essa temperatura em que o material começava a apresentar o novo es-

tado foi chamada, posteriormente, de temperatura de transição de fase, ou temperatura

crítica Tc, sendo essa a transição entre o estado supercondutor e o estado normal. Experi-

mentos posteriores utilizando chumbo e estanho mostraram que a supercondutividade não

era apenas inerente ao mercúrio, mas a vários outros metais, sendo que cada um possui

uma temperatura crítica associada. Onnes investigou também se a resistência no estado

supercondutor era exatamente zero ou aproximadamente zero através de um aparato que

consistia em uma corrente fluindo por um anel supercondutor. Essa análise viria a ser

14



Introdução 15

repetida com maior precisão por Mills e Files, sendo constatado que essas supercorren-

tes tinham um tempo de vida de mais de 100000 anos. Outra característica intrínseca a

esses materiais foi descrita por dois cientistas em 1933. Walter Meissner e Robert Och-

senfeld [5] notaram que no estado supercondutor esses materiais exibiam a propriedade

de excluir campos magnéticos não muito intensos de seu interior, i.e., eram diamagnetos

perfeitos. Essa propriedade ficou conhecida como efeito Meissner, uma das assinaturas

características da supercondutividade.

Desde sua descoberta, a supercondutividade mostrou-se um desafio para diversos ci-

entistas. Só em 1935 - aproximadamente 25 anos após a sua descoberta - dois físicos

teóricos (e também irmãos), Fritz e Heinz London [6] contribuíram com a primeira te-

oria fenomenológica de característica macroscópica, que daria uma descrição simples e

sucinta sobre o funcionamento desses materiais na condição supercondutora. Utilizando

as Equações de Maxwell, a teoria era capaz de explicar como a corrente e o campo mag-

nético atuam dentro de um supercondutor, fazendo relação com os fenômenos observados

por Onnes e Meissner. Deste modo, a teoria engloba com sucesso a resistência zero e o

diamagnetismo perfeito, mostrando a existência de um comprimento fundamental λ dado

como a profundidade que o campo penetra no material supercondutor, o qual depende da

temperatura T em que este se encontra.

Posteriormente, em 1950, surgiu uma importante teoria fenomenológica, baseada nas

ideias apresentadas pelos irmãos London e inspirada pela Teoria de Landau para tran-

sições de fase de segunda ordem. Elaborada por Vitaly Ginzburg e Lev Landau, essa

teoria ficou conhecida como a teoria de Ginzburg-Landau [7], conseguindo ir além em

sua consistência e apresentando resultados de forma mais elegante que sua antecessora.

Basicamente a teoria associa uma função de onda Ψ (r) (chamada de parâmetro de or-

dem) aos elétrons supercondutores e ao potencial vetor A (r). Esta teoria apresenta duas

equações características, uma descrevendo a variação do parâmetro de ordem no espaço



Introdução 16

e a outra o campo magnético local. A primeira equação contém um comprimento fun-

damental denominado comprimento de coerência ξ(T ), que denota a distância ao longo

da qual o parâmetro de ordem tem uma variação significativa. Já a segunda equação é

essencialmente a lei de Ampère, e também possui um comprimento fundamental denomi-

nado comprimento de penetração λ(T ). Uma das características importantes dessa teoria

é que ela consegue descrever a destruição da supercondutividade através da temperatura

e do campo magnético aplicado. A teoria ainda prevê a existência de dois tipos de su-

percondutores, simplesmente denominados tipo-I e tipo-II. A diferença básica pode ser

constatada da seguinte forma: os do tipo-I apresentam uma transição de fase do estado

normal para o estado supercondutor, enquanto os do tipo-II apresentam uma transição de

fase do estado normal para o estado misto e, na sequência, para o estado supercondutor,

onde no estado misto o material é supercondutor, mas possui regiões normais, i.e., regiões

onde ocorre penetração de fluxo magnético quantizado (possibilidade proposta em 1956

pelo físico russo Alexei A. Abrikosov [8]). Se considerarmos um parâmetro adimensional

denominado κ = λ(T )/ξ(T ), como sendo um parâmetro cujos valores abrangidos descre-

vem essa diferenciação, teremos que, para supercondutores tipo-I,
(
κ < 1/

√
2
)
, enquanto

para supercondutores tipo-II,
(
κ > 1/

√
2
)
, demonstrando que, para o limite de κ → ∞,

a teoria de London era recuperada.

Segundo Abrikosov, haveria a possibilidade de alguns tipos de supercondutores acei-

tarem uma penetração de fluxo de campo magnético em seu interior, sendo que essas

linhas de campo penetradas eram quantizadas e se organizavam em um arranjo ordenado,

de forma que a minimização de energia prevalecesse. Esse estado, posteriormente, ficou

conhecido como estado misto, coexistindo no material regiões supercondutoras e normais.

Nesse estado misto, o campo magnético fica aprisionado por supercorrentes elétricas que

o circundam, o que se denomina por vórtice.

Apesar da teoria de Ginzburg-Landau ter sido recebida com certa relutância pelos



Introdução 17

cientistas da época por suas características simples, pois ela não entrava na questão de

como o material se torna supercondutor, em 1959, Lev Gorkov [9] mostrou que ela po-

deria ser obtida como um caso particular de uma outra teoria proposta em 1957 chamada

BCS (sigla abreviada dos nomes de seus contribuidores, John Bardeen, Leon Cooper e

John Robert Schrieffer [10, 11]) que é atualmente a teoria mais conhecida que explica, do

ponto de vista microscópico, como um material se torna um supercondutor.

Em meados dos anos 80, K. A. Müller e J. G. Bednorz [12] publicaram a descoberta

do fenômeno supercondutor em altas temperaturas - em relação às ligas metálicas utiliza-

das até então - após o desenvolvimento de compostos cerâmicos à base de lantânio, bário,

cobre e oxigênio. Posteriormente, M. K. Wu [13] e colaboradores desenvolveram o pri-

meiro composto de uma família denominada YBCO, compostos por ítrio, bário, cobre e

oxigênio, que podem ser resfriados através de nitrogênio líquido em substituição ao hélio

líquido, já que esses compostos do tipo YBCO possuem Tc maior que as ligas de lantânio.

Percebe-se, então, uma viabilização do processo geral, uma vez que a obtenção do nitrogê-

nio líquido possui menor custo, se comparado ao hélio líquido, potencializando um maior

leque de aplicações. Além disso, o nitrogênio líquido é obtido por destilação fracionada

do ar líquido, enquanto o hélio geralmente está associado ao gás metano, proveniente de

bolsões subterrâneos, e sua liquefação exige aumento considerável da pressão desse gás,

constatando-se, também, que o rendimento, calculado em termos de energia necessária

para obtenção desses dois tipos de gases, é maior no caso do nitrogênio líquido.

Não demorou muito para que os pesquisadores começassem a propor várias tecno-

logias novas na área, desde cabos elétricos capazes de transportar energia elétrica por

longas distâncias sem apresentar perdas, até trens que levitam sobre trilhos. Outro tipo

de possibilidade são supercomputadores, cujas unidades centrais de processamento são

constituídas por junções Josephson [14], possibilitando a ocorrência de um tipo de dados

denominado qubit e aumento significativo do poder de processamento.



Introdução 18

1.2 O Fenômeno da Supercondutividade

A primeira característica inerente desses materiais foi descoberta por Kamerlingh Onnes

ao resfriar o mercúrio a uma temperatura próxima de 4, 2 K, na qual o material super-

condutor apresentava resistência elétrica de aproximadamente 10−6 Ω, como mostra a

Figura 1.1.

Figura 1.1: Gráfico histórico de Onnes: resistência do mercúrio sólido como função da
temperatura. [1]

Esses dados ilustram a primeira observação referente ao fenômeno, sendo constatado

para diversos metais e compostos intermetálicos que, a partir de uma certa temperatura

crítica Tc, a resistência elétrica cai abruptamente para zero, ocorrendo,então, a transição

de fase entre o estado supercondutor e o estado normal.

A segunda característica inerente a esses materiais foi descoberta por W. Meissner

e R. Ochsenfeld [5], podendo ser descrita como no experimento hipotético ilustrado na

Figura 1.2.

No supercondutor, a exclusão das linhas de campo no material ocorre independente-

mente da história magnética à qual foi submetido, ou seja, no primeiro caso o supercon-



Introdução 19

Figura 1.2: Experimentos ilustrando a diferença entre um supercondutor e um condutor
perfeito. Experimento 1: amostras resfriadas sem a presença de um campo magnético
aplicado. Experimento 2: amostras resfriadas na presença de um campo magnético apli-
cado. Fonte: Próprio autor.

dutor e o condutor hipotético perfeito são submetidos a um resfriamento e posteriormente

é aplicado um campo magnético externo (não muito forte). O resultado obtido é que,

independentemente do material, essas linhas de campo sempre serão expulsas. Já para

o segundo caso, ambos materiais são submetidos a um resfriamento e ao mesmo tempo

é aplicado um campo magnético externo que não seja muito forte, sendo verificado que

apenas o supercondutor exclui as linhas de campo aplicado. O resultado dessa experi-

ência ficou conhecida como efeito Meissner [5]. É interessante mencionar que o grau

da exclusão de fluxo magnético pode depender do material ou das condições em que foi

realizada a medida. Essa combinação entre resistência nula e efeito Meissner resulta em

uma distinção clara entre um supercondutor e um hipotético condutor perfeito.

A terceira característica inerente a esses materiais está vinculada ao campo magnético

aplicado, pois esse campo não pode ser arbitrariamente alto. O seu valor está intrinse-

camente relacionado a cada valor de temperatura, fazendo uma fronteira entre o estado



Introdução 20

supercondutor e o estado normal. Um gráfico geral para vários materiais está ilustrado na

Figura ??. Qualquer campo acima do campo crítico Hc(T ) leva o supercondutor para o

estado normal.

A quarta característica da supercondutividade está relacionada à descontinuidade no

calor específico nas proximidades da temperatura crítica TC , dada por:

CS − CN = µ0TC
dHC

dT
. (1.1)

Na equação 1.1 temos que a diferença entre o calor específico no estado normal (CN ) e

o calor específico no estado supercondutor (CS) é proporcional à variação do campo mag-

nético aplicado HC com relação à temperatura. Assim, em uma curva C(T ), verifica-se

uma descontinuidade do calor específico, como mostrado na Figura 1.3. Tal descontinui-

dade só ocorre na ausência de campo magnético, o que caracteriza a transição normal-

supercondutora como de segunda ordem.

Figura 1.3: (a) Gráfico geral da dependência do campo crítico Hc(T ) pela temperatura
para diversos materiais. [1]. (b) Curva do calor específico em função da temperatura
para um material supercondutor na ausência de campo magnético, onde verifica-se uma
descontinuidade. [2] (Adaptada da fonte original).



2
Teorias Fenomenológicas da Supercondutividade

2.1 Teoria Fenomenológica de London

2.1.1 Introdução

A descrição das teorias relacionadas neste capítulo teve como base os livros de Ketterson

e Song [15], Parks [16] e Tinkham [17].

As assinaturas fundamentais da supercondutividade são:

(a) resistência elétrica nula abaixo de Tc;

(b) descontinuidade no calor específico próximo a Tc;

(c) efeito Meissner-Ochsenfeld;

(d) os materiais supercondutores são classificados em dois tipos: os do tipo-I e do tipo-

II. Quando na presença de um campo magnético externo, os supercondutores do

tipo-I apresentam apenas o estado Meissner e o estado normal, sendo que estes são

delimitados por um determinado campo crítico que depende da temperatura,Hc(T ).

Além dos estados Meissner e normal, os do tipo-II possuem mais um estado, conhe-

cido como estado misto e caracterizado pela exclusão parcial de campo magnético.

Assim, os supercondutores do tipo-II têm dois campos críticos: um campo crítico

21



Teorias Fenomenológicas da Supercondutividade 22

inferior, Hc1(T ), que separa o estado Meissner do estado misto e um campo crítico

superior, Hc2(T ), que separa o estado misto do estado normal.

Foram necessárias quatro décadas para que estas quatro características fossem com-

preendidas. As duas primeiras só podem ser explicadas em bases mecânico quânticas

(teoria BCS [10]). A terceira propriedade pode ser explicada por meio de uma teoria fe-

nomenológica combinada com o eletromagnetismo e a termodinâmica (teoria de London).

Já a quarta propriedade pode ser explicada por uma combinação do eletromagnetismo com

uma teoria semi-clássica (teoria de Ginzburg-Landau). Neste capítulo descreveremos as

teorias de London e de Ginzburg-Landau.

2.1.2 Primeira Abordagem da Teoria de London

A primeira versão desta teoria considera o supercondutor como um condutor perfeito.

Tendo como base o modelo de Drude-Lorentz, temos que a ação de um campo elétrico E

sobre um elétron gera uma aceleração constante. Assim, a equação de movimento é dada

por:

m
dv

dt
= −m

τ
v + eE , (2.1)

onde m é a massa dos portadores de carga (elétrons), e é a carga elementar e v a sua

velocidade. A variável τ trata-se do tempo, medido entre uma colisão de um elétron

contra um ponto da rede e a próxima colisão deste mesmo elétron contra um outro ponto

qualquer. Assumindo um condutor perfeito, a contribuição do atrito tende a zero, pois

τ → ∞, restando apenas o segundo termo. Considerando uma densidade de corrente

supercondutora J = nev, onde n é a densidade desses elétrons (assumida uniforme), a

equação (2.1) pode ser reescrita da seguinte forma:

dJ

dt
=
ne2

m
E . (2.2)

A equação (2.2) é conhecida como a primeira equação de London e é equivalente

à segunda lei de Newton aplicada aos elétrons do supercondutor. Partindo das leis de



Teorias Fenomenológicas da Supercondutividade 23

Ampère (2.3) e de Faraday (2.4):

∇×B =
4π

c
J , (2.3)

∇× E = −1

c

∂B

∂t
, (2.4)

e substituindo (2.3) em (2.2) e aplicando o rotacional ao resultado e, em seguida,

usando (2.4), obtemos o seguinte resultado:

∂

∂t

(
∇×∇×B +

4πne2

mc2
B

)
= 0 . (2.5)

Com efeito, o argumento da derivada temporal é rigorosamente nulo ou uma cons-

tante. Para o momento, consideraremos a primeira hipótese1 a qual será suficiente para

explicar o estado Meissner. Temos que:

∇×∇×B +
4πne2

mc2
B = 0 . (2.6)

A constante que aparece na equação 2.6 tem unidades de inverso de quadrado de

comprimento, e é denominada por:

λ2 =
mc2

4πne2
, (2.7)

onde λ é conhecido como comprimento de penetração de London. Assim, podemos es-

crever que:

λ2∇×∇×B + B = 0 . (2.8)

A equação (2.8) é conhecida como a segunda equação de London. Conforme veremos

adiante, essa equação nos fornece informação a respeito do diamagnetismo dos materiais

supercondutores, mostrando que o campo magnético não pode penetrar no interior do

supercondutor além de uma camada superficial de espessura λ.

1O lado direito da equação 2.3 é considerada nula. Este argumento é suficiente para descrever a pene-
tração superficial do campo magnético. Posteriormente veremos que esta equação não é nula quando existir
nucleação de vórtices no supercondutor.



Teorias Fenomenológicas da Supercondutividade 24

2.1.3 Segunda Abordagem da Teoria de London

A segunda abordagem proposta pelos irmãos London é baseada em um modelo de dois

fluidos. O modelo de dois fluidos parte da ideia de que um deles é um fluido sem viscosi-

dade (superfluido) e o outro, um fluido com viscosidade finita (fluido normal). Assumindo

que a energia livre associada aos dois fluidos pode ser dividida em três contribuições, po-

demos escrever:

F = FN + FKin + FMag , (2.9)

onde FN é a energia do fluido no estado normal, FKin é a energia cinética do fluido em

movimento e FMag a energia armazenada pelo campo magnético. As segunda e terceira

contribuições são dadas por:

FKin =

∫
1

2
ρ (r) v2 d3r , (2.10)

e

FMag =
1

8π

∫
B (r) ·B (r) d3r , (2.11)

onde ρ (r) é a densidade específica referente ao superfluido. Tomando ρ = nm, v = 1
ne

J e

utilizando a Lei de Ampère (2.3), podemos rearranjar a equação (2.10) da seguinte forma:

FKin =
1

8π

∫
λ2∇×B ·∇×B d3r . (2.12)

Somando todas as contribuições para a energia livre, obtemos:

F = FN +
1

8π

∫ [
B ·B + λ2∇×B ·∇×B

]
d3r (2.13)

Podemos obter a equação de London (2.8) tratando B variacionalmente substituindo

o funcional de (2.13) na equação de Euler:

∂F

∂B
+ ∇× ∂F

∂∇×B
= 0 . (2.14)



Teorias Fenomenológicas da Supercondutividade 25

2.1.4 Efeito Meissner

Consideremos o caso simples de um supercondutor semi-infinito, onde, em x > 0, esteja

o supercondutor, e fora (em x < 0), consideremos que seja o vácuo. Ao imergirmos o

mesmo em um campo magnético H paralelo ao plano yz, teríamos B = (0, 0, Bz(x)) em

seu interior e, fora dele, B = Hk, onde H é o campo magnético externo aplicado. Tendo

em vista as considerações, a segunda equação de London pode ser escrita como:

d2Bz

dx2
− 1

λ2
Bz = 0 , (2.15)

onde usamos a identidade vetorial ∇× (∇×B) = ∇ · (∇ ·B)−∇2B, e a equação de

Maxwell ∇ ·B = 0. A solução desta equação é dada por:

B(x) = He−x/λ k . (2.16)

A Figura 2.1 ilustra este resultado onde vemos claramente que o campo magnético

local penetra em uma distância λ superfície adentro do material.

Figura 2.1: (a) Perfil do campo magnético local em um supercondutor semi-infinito na
presença de uma campo aplicado. Note que, na fronteira, Bz(0) = H . (b) Superfície
curva demonstrando a penetração de campo na direção do eixo x. Fonte: Próprio autor.

Substituindo o campo magnético local na lei de Ampère-Maxwell (2.3) encontramos

a seguinte expressão para a densidade de corrente induzida:

J =
c

4πλ
He−x/λ j (2.17)



Teorias Fenomenológicas da Supercondutividade 26

Essa equação nos permite visualizar como o campo magnético externo é impedido de

entrar no interior do supercondutor, sendo esta conhecida também como corrente de blin-

dagem.

Deste modo, a teoria de London consegue descrever, para temperaturas abaixo de Tc,

o caráter diamagnético perfeito de um material supercondutor e a penetração superficial

de campo magnético. Vale lembrar, ainda, que diferentemente da teoria de Ginzburg-

Landau, a teoria de London não descreve a maneira como podem surgir os vórtices, de

modo que, para fins de cálculos de minimização de energia - considerando esse fator e

a teoria de London - é necessário, primeiramente, considerar a presença do vórtice no

interior do material e, na sequência, proceder com a minimização de energia.

2.2 Teoria Fenomenológica de Ginzburg-Landau

Embora muito tempo depois da supercondutividade ter sido descoberta, não havia, ainda,

uma explicação satisfatória do ponto de vista microscópico com bases mecânico quân-

ticas. Vitaly Lazarevich Ginzburg e Lev Davidovich Landau propuseram uma teoria

semi-clássica que era capaz de descrever o fenômeno. Ginzburg e Landau notaram que a

transição do estado supercondutor nas proximidades da temperatura crítica Tc, podia ser

descrita em termos da teoria de Landau para transições de fase de segunda ordem baseada

em um parâmetro de ordem. Esse parâmetro apresenta um comportamento equivalente

à magnetização estudada na teoria de Landau (podendo ser tratado como um fenômeno

quântico macroscópico) se anulando na fase desordenada (T > Tc) ou sendo diferente de

zero na fase ordenada (T < Tc).

Na sequência, iremos partir da magnetização estudada na teoria de Landau para pos-

teriormente fazermos analogia com a supercondutividade.

Para a magnetização, o ponto de transição de fase é chamado de temperatura de Curie

Tc, sendo M diferente de zero abaixo dessa temperatura. Landau observou que o valor da



Teorias Fenomenológicas da Supercondutividade 27

magnetização era pequeno perto da transição de fase e apresentava variações suaves no

espaço até zero, de modo que a função poderia ser considerada contínua. A partir dessa

observação, a energia livre poderia ser escrita em uma série de potências da magnetização:

F = F0 + αM2 +
β

2
M4 , (2.18)

onde F0 é a energia para a magnetização igual a zero. Olhando para (2.18) nota-se que

a constante β é responsável pela concavidade da curva gerada pela função e que, se seu

valor for negativo, a energia será mínima para um valor arbitrariamente alto de magneti-

zação, fato esse que acarretaria em divergência. Portanto, o valor de β deve ser positivo

para que possamos ter um valor finito de energia mínima, gerando os gráficos (c) e (d) da

Figura 2.2. O mesmo já não ocorre para a constante α, que pode assumir tanto valores po-

sitivos quanto valores negativos, sendo α responsável pelas formas que a curva da função

(2.18) pode assumir, não implicando em divergência. Assim, temos que:

• se α > 0, o mínimo ocorre em M = 0;

• se α < 0, o mínimo ocorre em M 6= 0.

Na sequência, precisamos relacionar o parâmetro de ordem com a temperatura para

termos uma transição de fase de segunda ordem. Ainda, devemos considerar que o valor

de α deve ser negativo para termos um mínimo de magnetização diferente de zero, con-

forme ilustrado na Figura 2.2 (d). Assim, uma relação para α que satisfaz estes requisitos

pode ser escrita como:

α(T ) =

{
α0(T − Tc) , para T < Tc ,
0 , para T ≥ Tc ,

(2.19)

onde α0 é uma constante positiva.

Os valores de M que correspondem aos mínimos de energia podem ser obtidos por

meio da equação:
∂F

∂M
= [α(T ) + βM2]M = 0. (2.20)



Teorias Fenomenológicas da Supercondutividade 28

Figura 2.2: Perfis da energia livre. Em (a), (b) e (c), esses perfis não satisfazem a teoria
de Landau de transições de fase de segunda ordem, enquanto em (d), a teoria é satisfeita,
com mínimo ocorrendo em M 6= 0, indicado pela seta apontada sobre a curva de cor
verde. Fonte: Próprio autor.

Resolvendo a equação 2.20 obtemos duas soluções:

α > 0⇒M2 = 0 para T ≥ Tc ,

α ≤ 0⇒M2 = −α
β

= −α0(T − Tc)
β

para T < Tc . (2.21)

Notamos em (2.21) que, conforme a temperatura aumenta, o módulo de M diminui

até se anular em T = Tc. Ginzburg e Landau adotaram esta mesma ideia utilizada na

teoria de transição de fase de Landau, agora, no contexto da supercondutividade. Assim,

o parâmetro de ordem ψ(r), uma função de onda complexa, tal que o quadrado de seu

módulo representa a densidade de pares de Cooper [11] ns(r) = |ψ(r)|2 = ψ̄ψ, foi utili-

zada em substituição de M2. Ela é interpretada como uma função de onda dos portadores

de carga do supercondutor e a densidade de pares de Cooper, ns(r), não necessariamente

é homogênea no espaço (como no caso da teoria de London). Semelhantemente à magne-

tização da teoria de Landau, o parâmetro de ordem é tomado com valores pequenos perto

da transição de fase e com variações suaves no espaço até zero, de modo que a energia

livre possa ser expandida em uma série de potências do |ψ(r)|2 e do potencial vetor A.



Teorias Fenomenológicas da Supercondutividade 29

Portanto, podemos escrever:

F = F0 + α|ψ|2 +
β

2
|ψ|4 . (2.22)

Minimizando com relação a ψ̄, encontramos:

α > 0⇒ |ψ|2 = 0 para T > Tc ,

α ≤ 0⇒ |ψ|2 = −α
β

= −α0(T − Tc)
β

para T ≤ Tc . (2.23)

Devido às características de inomogeneidade do estado supercondutor, a quantidade

F pode ser tratada como energia livre e ser escrita da seguinte forma:

F = F0 +

∫
FC [ψ(r)] d3r , (2.24)

onde F0 corresponde à energia do estado normal, e FC = α|ψ|2 + β
2
|ψ|4, a energia do

condensado.

Efeitos relacionados ao comprimento de coerência2 não estão sendo levados em con-

sideração na equação acima, ou seja, isso significa que não está sendo considerado um

aumento de energia em relação a uma distorção espacial no parâmetro de ordem. Ginz-

burg e Landau consideraram esses efeitos e fizeram uma pequena correção na relação

acima, acrescentando um termo que seria o gradiente do parâmetro de ordem:

FG =
~2

2m∗
|∇ψ|2 , (2.25)

onde o termo ~2/2m∗ seria o termo associado à energia cinética na mecânica quântica e

m∗ a massa efetiva do par de Cooper.

A fim de levar em conta efeitos de aplicação de um campo magnético externo, Ginz-

burg e Landau acrescentaram na equação (2.25) o momento linear cinético:

p = −i~∇− e∗

c
A . (2.26)

2Este comprimento fundamental da teoria de Ginzburg-Landau será posteriormente estudado em maio-
res detalhes.



Teorias Fenomenológicas da Supercondutividade 30

onde e∗ é a carga efetiva do par de Cooper e A é o potencial vetor, que se relaciona com

o campo local por meio da equação B = ∇×A.

Introduzindo esta modificação na equação (2.25), temos que:

FG =
~2

2m∗

∣∣∣∣(∇− ie∗

~c
A(r)

)
ψ

∣∣∣∣2 . (2.27)

Por fim, é considerada a contribuição do campo magnético para a densidade de energia

como se segue:

FB =
1

8π
B2(r) , (2.28)

Portanto a energia total pode ser escrita como:

F = F0 +

∫
(FC + FG + FB) d3r

= F0 +

∫ {
α|ψ|2 +

β

2
|ψ|4 +

~2

2m∗

∣∣∣∣[∇− ie∗

~c
A

]
ψ

∣∣∣∣2 +
1

8π
B2

}
d3r .

(2.29)

Para minimizar o funcional de energia livre em (2.29) com relação a ψ̄ e A, podemos

utilizar as equações de Euler-Lagrange:

∂F

∂ψ̄
−∇ ·

[
∂F

∂(∇ψ̄)

]
= 0 ,

∂F

∂A
−∇×

[
∂F

∂(∇×A)

]
= 0 . (2.30)

De posse das relações acima, podemos finalmente escrever as duas equações de Ginzburg-

Landau:

− ~2

2m∗

(
∇− ie∗

~c
A

)2

ψ + αψ + β|ψ|2ψ = 0 , (2.31)

∇×∇×A =
e∗

2m∗

[
ψ̄

(
−i~∇− e∗

c
A

)
ψ + ψ

(
+i~∇− e∗

c
A

)
ψ̄

]
. (2.32)

O membro direito da segunda equação de Ginzburg-Landau (2.32) também pode ser

escrito como:

Js =
e∗

m∗
Re

[
ψ̄

(
~
i
∇− e∗

c
A

)
ψ

]
, (2.33)



Teorias Fenomenológicas da Supercondutividade 31

onde Js é a chamada de densidade de corrente supercondutora. Assim, a equação (2.32)

nada mais é do que a lei de Ampère ∇×B = 4π
c

Js.

2.3 Comprimentos Característicos

No contexto da teoria de Ginzburg-Landau, existem duas escalas de comprimentos fun-

damentais. O primeiro deles, é o comprimento de coerência ξ(T ), cujo valor é a variação

espacial do parâmetro de ordem. O segundo é o comprimento de penetração de London,

λ(T ), responsável pela variação do campo magnético no interior do supercondutor. No

que segue, mostramos que estes comprimentos são inerentes à teoria de Ginzburg-Landau.

2.3.1 Comprimento de Coerência

Consideraremos para esta primeira análise um caso envolvendo uma inomogeneidade do

parâmetro de ordem gerada pela presença de um contorno, na ausência de campos e cor-

rentes. Tomamos um supercondutor semi-infinito preenchendo o espaço x > 0, e escre-

vemos a primeira equação de Ginzburg-Landau (2.31) em sua forma unidimensional:

− ~2

2m∗
d2ψ

dx2
+ αψ + βψ3 = 0 , (2.34)

onde consideramos ψ real.

Para o estado supercondutor a constante α é negativa, então podemos escrever α =

−|α|. Substituindo na equação anterior, obtemos:

− ~2

2m∗
d2ψ

dx2
− |α|ψ + βψ3 = 0. (2.35)

É mais conveniente escrever o parâmetro de ordem na forma adimensional por meio da

transformação: ψ =
√
|α|/βf , onde α|/β é dado na equação (2.23). Além disso, tendo

em vista que trata-se de uma amostra semi-infinita, utilizamos as seguintes condições de

contorno: para x→∞, a função f 2 → 1, pois em uma posição longe da fronteira a amos-

tra é perfeitamente supercondutora, e ainda, para x → ∞, podemos também considerar



Teorias Fenomenológicas da Supercondutividade 32

f ′ = df/dx→ 0. Assim, dadas tais condições, temos:

− ~2

2m∗|α|
d2f

dx2
− f + f 3 = 0 , (2.36)

onde o coeficiente da derivada segunda tem unidades de quadrado de comprimento e é

conhecido como comprimento de coerência:

ξ2 =
~2

2m∗|α|
. (2.37)

Na sequência iremos multiplicar a equação (2.36) por f ′, integrar seu resultado e, em

seguida, derivar com relação a x, obtendo:

d

dx

[
− ξ2(f ′)2

2
− f 2

2
+
f 4

4

]
= 0 . (2.38)

Consequentemente, temos que:

− ξ2(f ′)2

2
− f 2

2
+
f 4

4
= C . (2.39)

Usando as condições de contorno, encontramos que C = −1/4. Com efeito,

ξ2(f ′)2 =
1

2
(1− f 2)2 , (2.40)

que resolvida por integração direta nos fornece a seguinte relação:

f = tanh

(
x√
2ξ

)
. (2.41)

A Figura 2.3 mostra a variação do parâmetro de ordem nas vizinhanças da interface

metal-supercondutor conforme a equação (2.41).

Finalmente, podemos obter a relação do comprimento de coerência com a tempera-

tura. Uma vez que α(T ) = −α0Tc[1− (T/Tc)], substituindo na equação (2.37) encontra-

mos:

ξ(T ) =

√
~2

2m∗α0Tc

(
1− T

Tc

)−1/2

= ξ(0)

(
1− T

Tc

)−1/2

. (2.42)



Teorias Fenomenológicas da Supercondutividade 33

Figura 2.3: Ilustração do comportamento do parâmetro de ordem próximo à interface
metal-supercondutor. Após o comprimento ξ, o material encontra-se no estado supercon-
dutor com valor máximo de parâmetro de ordem.

2.3.2 Comprimento de Penetração de London

Para esta segunda análise, com fins de encontrarmos o valor do comprimento de penetra-

ção de London, utilizaremos outras condições de contorno, quais sejam: admitir a existên-

cia de um campo magnético aplicado, mas uma variação espacial desprezível do módulo

do parâmetro de ordem. Em outras palavras, assumiremos que ns = |ψ|2 = |α|/β em

todo espaço. A partir dessas considerações, a equação (2.33) para densidade de corrente

supercondutora torna-se:

Js = −(e∗)2

m∗c
|ψ|2A. (2.43)

Usando a lei de Ampère ∇ × B = 4π
c

Js, e substituindo esta na equação (2.43),

encontramos:
m∗c2β

4π(e∗)2|α|
∇×B + A = 0. (2.44)

Por fim, aplicando o rotacional a ambos os lados da equação (2.44), encontramos:

λ2∇×∇×B + B = 0 (2.45)

Esta equação nada mais é do que a segunda equação de London, com a diferença de



Teorias Fenomenológicas da Supercondutividade 34

que, agora, o comprimento de penetração depende da temperatura e é escrito como segue:

λ2(T ) =
m∗c2β

4π(e∗)2|α|

=
m∗c2β

4π(e∗)2|α|

(
1− T

Tc

)−1

= λ2(0)

(
1− T

Tc

)−1

. (2.46)

2.3.3 Parâmetro de Ginzburg-Landau

Embora os comprimentos fundamentais dependam da temperatura, a razão entre eles re-

sulta em uma constante que independe de T . Esta constante é conhecida como parâmetro

de Ginzburg-Landau, e é dada por:

κ =
λ(T )

ξ(T )
=
λ(0)

ξ(0)
=
m∗c

e∗~

√
β

2π
. (2.47)

O valor κ é de grande importância na teoria de Ginzburg-Landau no que diz respeito à

classificação dos materiais supercondutores. Discutiremos este ponto em maiores detalhes

na próxima seção.

2.4 Classificação dos Supercondutores

A classificação dos supercondutores pode ser feita dependendo de seu comportamento na

presença de um campo magnético externo aplicado. Esta divisão é baseada no fato de que

a energia de superfície é proporcional à diferença (ξ − λ), considerando a fronteira entre

um supercondutor e um material normal. Cálculos aproximados mostram que a diferença

entre a energia do supercondutor e a do material normal é dada por γ = α2

2β
(ξ − λ).

Sendo assim, se ξ > λ, a energia superficial é sempre positiva, indicando que o material

supercondutor não é propenso à formação de domínios normais, permanecendo no estado

Meissner até um determinado campo crítico, acima do qual se torna normal. Por outro



Teorias Fenomenológicas da Supercondutividade 35

lado, se ξ < λ, a energia superficial é negativa, e o supercondutor é propenso à formação

de domínios normais em seu interior, e o fluxo magnético penetra em pequenos tubos

(vórtices) em quantidades quantizadas de fluxo dadas por Φ0 = hc/2e = 2, 07× 10−7 G ·

cm2. No primeiro caso, classificamos o supercondutor como sendo do tipo-I, e na segunda

situação como sendo do tipo-II. Cálculos mais exatos mostram que materiais do tipo-I

exibem κ < 1/
√

2 , enquanto materiais do tipo-II exibem κ > 1/
√

2. A Figura 2.4 ilustra

os perfis de campo e parâmetro de ordem para ambos os tipos de supercondutor.

Figura 2.4: Representação esquemática de perfis de campo e de parâmetro de ordem para
supercondutores do tipo-I e II, onde a configuração em (b) permite a formação de um
vórtice, enquanto a configuração em (a) não. Fonte: Próprio autor.

Os supercondutores do tipo-I de dimensões muito maiores que os comprimentos fun-

damentais ξ e λ permanecem no estado Meissner para campos aplicados até atingir o

campo crítico termodinâmico Hc(T ) dado por:

Hc(T ) =
Φ0

2
√

2πλ(T )ξ(T )
. (2.48)

Neste estado de diamagnetismo perfeito todo campo magnético é expelido do interior

da amostra. Acima de Hc(T ) a supercondutividade não se mantém e a amostra vai para o

estado normal (ver Figura 2.5).

Supercondutores do tipo-II, também de dimensões muito maiores que os comprimen-

tos fundamentais, permanecem no estado Meissner para campos menores que o primeiro



Teorias Fenomenológicas da Supercondutividade 36

Figura 2.5: Diagrama de fase de um supercondutor do tipo-I. Fonte: Próprio autor.

campo crítico Hc1(T ), cuja expressão é dada por:

Hc1(T ) =
Φ0

4πλ2(T )
lnκ . (2.49)

Para campos maiores que Hc1(T ), vórtices começam a nuclear no supercondutor, até

que a densidade destes torna-se muito alta, a ponto de destruir a supercondutividade. Este

valor de H é conhecido como segundo campo crítico, dado por:

Hc2(T ) =
Φ0

2πξ2(T )
. (2.50)

A região Hc1(T ) < H < Hc2(T ) é conhecida como estado misto, onde as regiões

normal e supercondutora coexistem no interior do material. Para H = Hc2(T ) a super-

condutividade é destruída, mas ainda há regiões supercondutoras remanescentes em um

fina camada próxima à superfície do material. Aumentando-se o campo aplicado H até

atingir o terceiro campo crítico Hc3(T ), a supercondutividade é então totalmente supri-

mida. Para um supercondutor semi-infinito,Hc3(T ) = 1.69Hc2(T ). Já para um filme com

o campo aplicado paralelamente ao plano do filme, Hc3(T ) = 2Hc2(T ).

Do ponto de vista experimental, as diferentes fases de um supercondutor podem ser



Teorias Fenomenológicas da Supercondutividade 37

Figura 2.6: Diagrama de fase de um supercondutor do tipo-II. Fonte: Próprio autor.

identificadas por meio da medida da magnetização do estado de equilíbrio. Desprezando

efeitos de desmagnetização3, podemos escrever:

M =
B̄−H

4π
. (2.51)

onde B̄ é a indução magnética média4.

A Figura 2.7 (a) mostra a magnetização, −M, como função do campo aplicado, H,

para um supercondutor do tipo-I. No estado de Meissner o fluxo é expelido (B̄ = 0) do

interior da amostra, de tal modo que −M = H/4π tendo um comportamento linear com

o campo aplicado. Tão logo o campo aplicado alcance o valor de Hc(T ), a magnetização

decai abruptamente para zero (ver Figura 2.7).

A Figura 2.7 (b) exibe a magnetização de um supercondutor do tipo-II. No intervalo

H < Hc1(T ) temos o estado de Meissner completo (B̄ = 0). Acima do primeiro campo

crítico Hc1(T ) observa-se penetração parcial de fluxo magnético (B̄ 6= 0) até atingir um

campo crítico superior Hc2(T ), valor a partir do qual o material volta ao estado normal

3Provocados pela curvatura das linhas de campo magnético em torno do supercondutor.
4Média espacial do campo local B.



Teorias Fenomenológicas da Supercondutividade 38

anulando a magnetização (B̄ = H).

Figura 2.7: (a) Magnetização de um supercondutor do tipo-I e (b) Magnetização de um
supercondutor do tipo-II, ambos em função do campo aplicado H . A linha tracejada
representa um supercondutor do tipo-I. Fonte: Dissertação de Mestrado de Mauro César
Vieira Pascolati. Unesp Bauru/POSMAT - Ano 2010.

2.5 Estrutura de Vórtices em Supercondutores do Tipo-
II

2.5.1 Quantização do Fluxo

A quantização do fluxo pode ser facilmente determinada por meio da segunda equação da

densidade de corrente (2.33), a qual pode ser escrita introduzindo o parâmetro de ordem

em termos de sua magnitude e fase, ψ = |ψ|eiφ:

Js =
2e~
m∗
|ψ|2∇φ− 4e2

m∗c
A|ψ|2 . (2.52)

Esta equação também pode ser reescrita na forma:

A =
~c
2e

∇φ− m∗cJs
4e2|ψ|2

. (2.53)

Agora, usando o teorema de Stokes, vem que:∮
C

A · dr =

∫
S

∇×A · ndS =

∫
S

B · ndS = Φ , (2.54)



Teorias Fenomenológicas da Supercondutividade 39

onde S é uma superfície bilateral restrita ao caminho fechado C. Assim, integrando o

potencial vetor da equação (2.53), encontramos:

~c
2e

∮
C

∇φ · dr− m∗c

4e2

∮
C

Js
|ψ|2
· dr = Φ. (2.55)

Uma vez que a fase φ do parâmetro de ordem complexo é uma função unívoca (em inglês,

single-values function), isto impõe que a fase deve mudar em valores de múltiplos inteiros

de 2π em um círculo fechado C. Esta condição é conhecida como quantização de Bohr-

Sommerfeld. Assim, temos que:

Φ = n
hc

2e
− m∗c

4e2

∫
C

J

|ψ|2
· dr. (2.56)

onde n é um número inteiro e conhecido como número quântico de fluxóide e hc/2e = Φ0

é o quantum de fluxo magnético. Desconsiderando a penetração superficial de fluxo,

é possível encontrarmos um caminho C tal que a integral na equação (2.55) se anule.

Assim, temos que:

Φ = n
hc

2e
= nΦ0 , (2.57)

A equação (2.57) mostra que o fluxo penetrado no supercondutor é quantizado.

2.5.2 Vórtice Isolado

Conforme mencionamos anteriormente, em supercondutores do tipo-II, paraH > Hc1(T ),

vórtices penetram no interior do material. Cada vórtice tem um núcleo constituído de ma-

terial predominantemente normal de raio ξ(T ), distância na qual a densidade supercondu-

tora de superelétrons decai a zero no seu centro (ver Figura 2.8). Já o campo local B(r)

é máximo no centro do vórtice, e decai monotonicamente em uma distância λ(T ) devido

às supercorrentes de blindagem Js(r) ao redor do núcleo. Supondo que a separação en-

tre os vórtices seja grande se comparada com λ(T ), a superposição das supercorrentes é

desprezível. Assim, podemos desconsiderar a interação entre os vórtices. Neste limite,

podemos considerar o vórtice como um objeto isolado, e |ψ| pode ser considerada como



Teorias Fenomenológicas da Supercondutividade 40

uma constante em todo espaço, exceto no centro do núcleo. Contudo, a variação espacial

da fase será considerada. Considerando, ainda, a lei de Ampère (2.3) e a expressão da

densidade de corrente supercondutora (2.33), e usando a condição topológica da fase 5,

vem:

∇×∇φ = Φ0δ(r) , (2.58)

obtemos a equação de London:6

− λ2∇2B + B = kΦ0δ(r) . (2.59)

Figura 2.8: Estrutura de um vórtice. Fonte: Dissertação de Mestrado de Mauro César
Vieira Pascolati. Unesp Bauru/POSMAT - Ano 2010.

A solução exata da equação (2.59) pode ser obtida por meio da técnica da transfor-

mada de Fourier. A solução é dada por:

B =
Φ0

2πλ2
K0(r/λ)k , (2.60)

onde K0(x) é função de Bessel [18, 19] modificada de ordem zero; K0 decresce como

(
√
πλ/2r)e−r/λ para r � λ e diverge logaritmicamente como ln(λ/r) para r → 0.

Usando a lei de Ampère (2.3) na equação (2.60), obtemos a densidade de corrente que

flui ao redor do vórtice:

Js =
cΦ0

4π2λ2
K1(r/λ)θ , (2.61)

5Ocorrência de entrada de um quantum de fluxo a cada mudança 2π da fase
6Conforme havíamos observado no rodapé da página 23, o lado direito da equação (2.8) corresponde

justamente a esta singularidade da fase do parâmetro de ordem.



Teorias Fenomenológicas da Supercondutividade 41

ondeK1 é a função de Bessel modificada de primeira ordem a qual diverge como 1/r para

r → 0, e daí, como (
√
πλ/2r)e−r/λ para r � λ. É possível ainda determinar a energia

livre (por unidade de comprimento l do sistema) de um vórtice isolado resultando em:

F

L
=

(
Φ0

4πλ

)2

lnκ . (2.62)

Como podemos observar, a energia depende do quadrado de Φ0. Assim, não é favorável

a nucleação de vórtices com mais de um quantum de fluxo.

2.5.3 Interação entre Vórtices

Consideremos, agora, um conjunto de vórtices, cada qual na posição {ri , i = 1, . . . , N}.

A equação de London (2.59) então pode ser escrita como:

− λ2∇2B + B = Φ0

N∑
i=1

δ(r− ri)k . (2.63)

Usando técnicas de transformada de Fourier, podemos mostrar que a solução desta

equação é dada por:

B =
Φ0

2πλ2

N∑
i=1

K0(|r− ri|/λ)k . (2.64)

Substituindo este resultado na energia livre de London (2.13) encontramos:

F

L
=
FN
L

+

(
Φ0

4πλ

)2 N∑
i=1

N∑
j=1

K0(|ri − rj|/λ) . (2.65)

Nesta equação estão contabilizadas tanto as autoenergias de cada vórtice individual-

mente, como a energia de interação entre eles. Em 1957, foi primeiramente demonstrado

por Abrikosov que as posições dos vórtices que minimizam a energia de interação entre

eles corresponde a uma rede triangular conforme ilustrado na Figura 2.9. A presente aná-

lise restringe-se ao caso em que os vórtices estão suficientemente distantes um do outro.

Em outras palavras, o limite de campo aplicado de validade é Hc1(T ) < H � Hc2(T ).

Abrikosov também mostrou que esta rede é a que minimiza a energia livre para campos



Teorias Fenomenológicas da Supercondutividade 42

próximos de Hc2(T ). As previsões teóricas de Abrikosov foram posteriormente compro-

vadas experimentalmente usando-se técnicas de decoração magnética com o auxílio de

microscopia eletrônica.

Figura 2.9: Rede de vórtices de Abrikosov: perfil do campo local. Experimento com ma-
terial supercondutor, onde partículas ferromagnéticas foram aplicadas à superfície deste
material, permanecendo nesta configuração. As regiões mais escuras correspondem aos
valores mais altos de campo (região normal), enquanto as regiões mais claras, aos valores
mais baixos de campo (região supercondutora)[1].

2.6 Supercondutividade Mesoscópica

Denominamos por mesoscópico os materiais supercondutores com tamanho da ordem de

λ(T ) e/ou ξ(T ). As propriedades eletrônicas e estruturais da rede de vórtices nestes ma-

teriais são consideravelmente influenciadas pelos efeitos de confinamento. Com efeito, a

nucleação do estado de vórtices depende fortemente das condições de contorno impostas

ao material, isto é, da topologia do sistema. Tem-se observado que a matéria de vórtices

comporta-se de maneira muito diferente da mesma em supercondutores de dimensões ma-

croscópicas. Nos supercondutores mesoscópicos há uma competição entre a rede de vór-

tices triangular de Abrikosov e a geometria da amostra. Tem sido observado também que

é possível a ocorrência de vórtices gigantes (um vórtice com vorticidade superior a um).

Do ponto de vista teórico, os sistemas mesoscópicos têm sido estudados resolvendo-se as

equações de Ginzburg-Landau numericamente. No próximo capítulo iremos desenvolver

um algoritimo de solução numérica destas equações e, posteriormente, iremos aplicá-lo

para estudar um supercondutor de geometria SQUID com uma constrição.



3
O Método ψ − U

3.1 Introdução

Neste Capítulo apresentamos um método de discretização das equações de Ginzburg-

Landau dependentes do tempo (TDGL, time dependent Ginzburg-Landau). Em primeiro

lugar, escreveremos estas equações em unidades arbitrárias por meio de um sistema de

unidades reduzidas. O segundo passo será introduzir nas equações TDGL os campos au-

xiliares [20, 21], os quais tornam a primeira destas equações muito semelhante à equação

de difusão. Finalmente, apresentamos a malha de discretização das equações TDGL e

deduzimos fórmulas de aproximação, as quais podem ser implementadas computacional-

mente. O arcabouço de aproximações destas equações é conhecido como método ψ − U

ou também como método das variáveis de ligação [20–23]. Este Capítulo baseia-se fun-

damentalmente no trabalho de Gropp e outros [20].

3.2 Equações TDGL

As propriedades físicas dos materiais supercondutores com baixo Tc são descritas pelas

equações TDGL. As funções incógnitas dessas equações são o parâmetro de ordem com-

plexo ψ, e o potencial vetor A. O valor absoluto |ψ|2 representa a densidade de pares

de Cooper, e o campo local B é encontrado através da relação: B = ∇ ×A. Em 1966,

43



O Método ψ − U 44

Schmid [24] propôs a seguinte extensão das equações de Ginzburg-Landau estáticas (con-

forme vistas no Capítulo 2), para problemas dinâmicos:

~
2m∗D

(
∂

∂t
+ i

e∗

~
ϕ

)
ψ = − 1

2m∗
Π2ψ − αψ − β|ψ|2ψ , (3.1)

4πσ

c

(
1

c

∂A

∂t
+ ∇ϕ

)
=

4π

c
Js −∇×∇×A , (3.2)

onde ϕ é o potencial escalar, Π = (−i~∇− e∗

c
A) é o operador derivada covariante e Js

é a densidade de corrente supercondutora dada por:

Js =
e∗

m∗
Re
[
ψ̄Πψ

]
, (3.3)

onde Re denota a parte real de uma quantidade complexa.

A equação (3.2) é equivalente à lei de Ampère.

As constantes das equações (3.1), (3.2) e (3.3) têm os seguintes significados:

• c é a velocidade da luz;

• m∗ = 2m é a massa efetiva;

• e∗ = 2e é a carga efetiva dos super elétrons;

• D é o coeficiente de difusão;

• σ é a condutividade elétrica;

• β e α são constantes fenomenológicas, sendo que β não depende da temperatura, e

α = α0(T − Tc) para toda T < Tc e α = 0 para T > Tc;

• Tc é a temperatura crítica.

A constante β é definida positiva e α negativa.



O Método ψ − U 45

3.3 Transformações de Calibre

As equações TDGL são invariantes sob transformações de simetria (transformações de

calibre, ou "gauge invariance"), o que é importante, pois significa que descrevem uma

simetria global. Essas transformações são as seguintes:

ψ′ = ψeiχ ,

A′ = A +
~c
e∗

∇χ ,

ϕ′ = ϕ− ~
e∗
∂χ

∂t
. (3.4)

Pode-se mostrar que sob estas transformações, as equações TDGL permanecem inva-

riantes com as novas funções em lugar das antigas [23].

Para exemplificar como demonstramos a invariância de calibre, vamos tomar a equa-

ção (3.2). Usando a propriedade vetorial de que o rotacional do gradiente é nulo, ∇ ×

∇χ = 0, temos que B′ = ∇×A′ = ∇×A = B. Para a equação da densidade de cor-

rente usamos a identidade Πx (e−iχf) = e−iχΠ′xf , onde Π′x =
(
−i~ ∂

∂x
− e∗

c
A′x
)
. Assim,

por exemplo, temos que ψ̄Πxψ = eiχψ̄′Πx (e−iχψ′) = ψ̄′Π′xψ
′. Então, a densidade de cor-

rente também é um invariante de calibre, J′s = Js. Para o lado esquerdo da equação (3.2),

temos que(
1
c
∂A
∂t

+ ∇ϕ
)

=
(

1
c
∂
∂t

(
A′ − ~c

e∗
∇χ
)

+ ∇
(
ϕ′ + ~

e∗
∂χ
∂t

))
=
(

1
c
∂A′

∂t
+ ∇ϕ′

)
. Isto finaliza a

demonstração da invariância de calibre.

3.4 Estado Meissner

Suponhamos que não haja campo aplicado ao sistema supercondutor. Então, temos a

situação de equilíbrio:

ψ
(
α + β|ψ|2

)
= 0 . (3.5)



O Método ψ − U 46

Podemos ter duas possíveis soluções para a equação (3.5): (a) ψ = 0 ou (b) |ψ|2 =

−α
β

. Consideremos apenas a segunda solução. Considerando o parâmetro de ordem real,

obtemos ψ2 = ψ2
0 = −α

β
. Esta solução para o parâmetro de ordem é conhecida como

estado Meissner. Note que as condições α < 0 e β > 0 implicam ψ2
0 > 0, solução a qual

garante um mínimo de energia para o estado Meissner, conforme visto no Capítulo 2.

3.5 Unidades Reduzidas

Para resolvermos as equações TDGL numericamente é mais conveniente escrevermos es-

tas em unidades reduzidas, pois teremos um menor número de parâmetros para trabalhar.

As unidades que usaremos estão especificadas em (3.6). Em seguida, exemplificamos

um caso em que elas aparecem naturalmente nas equações TDGL quando usamos como

referência o estado Meissner. Temos que:

ψ = ψ0ψ̃ ,

T = TcT̃ ,

∇ =
1

ξ(0)
∇̃ ,

t = t0t̃ ,

A = Hc2(0)ξ(0)Ã ,

ϕ =
Hc2(0)D

c
ϕ̃ . (3.6)

onde t0 = ξ2(0)/D.

Para determinarmos a primeira equação TDGL (3.1) em unidade reduzidas, procede-

mos da seguinte forma. Primeiramente, multiplicamos ambos os lados desta equação por

ψ0 e consideramos cada um dos termos separadamente.

1. Termo do condensado: substituindo as componentes do parâmetro de ordem temos



O Método ψ − U 47

que:

−ψ0

(
αψ + βψ|ψ|2

)
= −ψ0

(
αψ0ψ̃ + β[ψ0ψ̃|ψ0ψ̃|2]

)
= −ψ2

0

(
αψ̃ + β[ψ2

0ψ̃|ψ̃|2]
)

= −ψ2
0

(
αψ̃ − ��β

α

��β
ψ̃|ψ̃|2

)
= ψ2

0α0(Tc − T )ψ̃
(

1− |ψ̃|2
)

= α0Tcψ
2
0︸ ︷︷ ︸(1− T̃ )ψ̃

(
1− |ψ̃|2

)
(3.7)

2. Termo de energia cinética:

ψ0
1

2m∗

(
−i~∇− e∗

c
A

)2

ψ = ψ2
0

~2

2m∗

(
−i∇− e∗

~c
A

)2

ψ̃

= ψ2
0

~2

2m∗
1

ξ2(0)

(
−i∇̃− e∗ξ(0)

~c
A

)2

ψ̃ ,

ψ2
0

~2

2m∗
1

ξ(0)2
= ψ2

0
��~2

���2m∗
1
��~2

��2m∗α0Tc

= α0Tcψ
2
0 , (3.8)

e∗ξ(0)

~c
A =

2e
hc
2π

ξ(0)A =
2π
hc
2e

ξ(0)A

=
2π

Φ0

ξ(0)2 1

ξ(0)
A

=
1
Φ0

2πξ(0)2

1

ξ(0)
A =

1

Hc2(0)ξ(0)
A = Ã , (3.9)

ψ0
1

2m∗

(
−i~∇− e∗

c
A

)2

ψ = α0Tcψ
2
0

(
−i∇̃− Ã

)2

ψ̃

= α0Tcψ
2
0︸ ︷︷ ︸Π2ψ̃ . (3.10)

3. Derivada temporal:

ψ0
~2

2m∗D

∂ψ

∂t
= α0Tcψ

2
0︸ ︷︷ ︸ ∂ψ̃∂t̃ .



O Método ψ − U 48

4. Potencial:

ψ0
~

2m∗D

(
i
e∗

~
ϕ

)
ψ = α0Tcψ

2
0︸ ︷︷ ︸ iϕ̃ψ̃ , (3.11)

Colecionando todos os termos obtemos:

∂ψ̃

∂t̃
+ iϕ̃ψ̃ = −Π̃

2
ψ̃ + (1− T̃ )ψ̃

(
1− |ψ̃|2

)
. (3.12)

onde, agora, a derivada covariante em unidades adimensionais é dada por Π̃ = −i∇̃−Ã.

Procedendo de forma análoga, podemos escrever para as outras equações:

β

(
∂Ã

∂t̃
+ ∇̃ϕ̃

)
= J̃s − κ2∇̃× B̃ , (3.13)

J̃s = Re
(

¯̃ψΠ̃ψ̃
)
, (3.14)

onde β = 4πσDκ2/c2.

Daqui em diante desconsideraremos os til’s em todas as funções.

3.6 Campos Auxiliares

Agora iremos introduzir o vetor campo auxiliar U = (Ux,Uy,Uz), cujas componentes são

definidas por:

Ux(x, y, z, t) = exp

(
−i
∫ x

x0

Ax(ν, y, z, t) dν

)
, (3.15)

Uy(x, y, z, t) = exp

(
−i
∫ y

y0

Ay(x, η, z, t) dη

)
, (3.16)

Uz(x, y, z, t) = exp

(
−i
∫ z

z0

Az(x, y, γ, t) dγ

)
, (3.17)

onde (x0, y0, z0) é um ponto de referência arbitrário. Daqui em diante omitiremos a de-

pendência temporal e recuperaremos essa dependência quando necessário.

Derivando as equações (3.15-3.17) com relação à coordenada x para (3.15), y para

(3.16) ou z para (3.17), temos que ∂Uζ
∂ζ

= −iAζUζ , onde {ζ = x, y, z}. Usando o fato que



O Método ψ − U 49

os campos auxiliares são funções unimodulares, ŪζUζ = 1, encontramos:

Πζf = iŪζ
∂(Uζf)

∂ζ
, (3.18)

onde f é qualquer função complexa.

Analogamente, podemos mostrar que:

Π2
ζf = Ūζ

∂2 (Uζf)

∂ζ2
. (3.19)

Usando a propriedade (3.18) e a identidade Re(−iw) = Im(w), onde w é uma função

complexa, as componentes da densidade de corrente (3.14) tornam-se:

Jsζ = Im

[
Ūζψ̄

∂(Uζψ)

∂ζ

]
. (3.20)

Em seguida usamos a propriedade (3.19) nos termos cinéticos da equação (3.12), e

podemos escrever:

∂ψ

∂t
+ iϕψ = Ūx

∂2 (Uxψ)

∂x2
+ Ūy

∂2 (Uyψ)

∂y2
+ Ūz

∂2 (Uzψ)

∂z2
+ (1−T )ψ

(
1− |ψ|2

)
. (3.21)

A vantagem de escrever as equações TDGL em termos dos campos auxiliares é que

elas não perdem a invariância de calibre quando são discretizadas.

3.7 Problema a Ser Resolvido: Condições de Contorno

Denotemos por Ωsc a região supercondutora. O meio exterior pode ser o vácuo ou um

metal, embora consideremos apenas o primeiro caso. Seja ∂Ωsc a interface separando os

dois meios. Consideremos agora um domínio maior Ω tal que contenha a região super-

condutora, isto é, Ωsc ⊂ Ω. A Figura 3.1 ilustra a geometria do problema. A interface

vácuo-vácuo mais externa é denotada por ∂Ω. Assim, as equações a serem resolvidas são

dadas por:

∂ψ

∂t
+ iϕψ =

z∑
ζ=x

Ūζ
∂2 (Uζψ)

∂ζ2
+ (1− T )ψ

(
1− |ψ|2

)
, em Ωsc , (3.22)



O Método ψ − U 50

Figura 3.1: Representação da caixa de simulação com o supercondutor no seu interior.
Fonte: Próprio autor.

β

(
∂A

∂t
+ ∇ϕ

)
=

{
Js − κ2∇×B , em Ωsc ,
−κ2∇×B , em Ω\Ω̄sc ,

(3.23)

com Jsζ dada pela equação (3.20).

A solução das equações TDGL precisam das condições de contorno na interface sepa-

rando os meios para serem resolvidas. As condições físicas que vamos assumir são as de

que as componentes perpendiculares de Js à interface são nulas. Isto significa que os pa-

res de Cooper estão restritos ao supercondutor. Em outras palavras, não há corrente entre

o domínio ΩSC e Ω. Seja n um vetor normal à interface supercondutor-vácuo apontando

para fora do supercondutor. Então temos que:

Js · n = 0 . (3.24)

Portanto, em se tratando de um problema com simetria retangular, as condições de



O Método ψ − U 51

contorno (3.24) tornam-se do tipo von Neumann:

∂(Uxψ)

∂x
= 0 , faces da esquerda e da direita , (3.25)

∂(Uyψ)

∂y
= 0 , faces de frente e de fundo , (3.26)

∂(Uzψ)

∂z
= 0 , faces inferior e superior , (3.27)

onde usamos a equação (3.20). Como veremos mais adiante, estas condições fornecem os

valores de ψ em ∂Ωsc.

Precisamos ainda das condições para o campo local. As interfaces ∂Ωsc e ∂Ω serão

consideradas suficientemente distantes, de tal modo que o campo magnético local B seja

igual ao campo aplicado H em ∂Ω. Assim, nesta interface assumimos a continuidade do

campo:

Bx = Hx , By = Hy , Bz = Hz , (3.28)

ou ainda,

∇×A = H , (3.29)

em ∂Ω.

Como um exemplo, podemos imaginar um supercondutor em forma de um paralele-

pípedo dentro de um domínio de mesma geometria e concêntricos.

Nas próximas seções vamos discretizar as equações TDGL (3.22,3.23,3.20), e as con-

dições de contorno (3.25-3.27,3.28).

3.8 Discretização das Equações TDGL

3.8.1 Malha de Discretização

Para discretizar as equações de Ginzburg-Landau, substitui-se as variáveis (x, y, z) por

suas discretas equivalentes (xi, yj, zk). Considerando (A ,B,C ) os lados de um parale-

lepípedo, dividimos o mesmo em intervalos de tamanhos (∆x,∆y,∆z). Assim, os pontos



O Método ψ − U 52

de vértice da malha de discretização são dados por xi = (i − 1)∆x, yj = (j − 1)∆y e

zk = (k − 1)∆k, para todos i = 1, . . . , Nx + 1, j = 1, . . . , Ny + 1, k = 1, . . . , Nz + 1,

onde Nx = A /∆x, Ny = B/∆y e Nz = C /∆z são os números de células unitárias nas

direções x, y e z, respectivamente. A Figura 3.2 ilustra a malha de discretização.

O domínio Ω encontra-se dento da caixa de simulação de dimensões (A ,B,C ) com

as suas faces a meia distância (∆x/2,∆y/2,∆z/2) de cada face do paralelepípedo nas

direções (x, y, z), respectivamente. Temos que:

Ω =

{
r = (x, y, z) : x1 +

∆x

2
< x < xNx +

∆x

2
, y1 +

∆y

2
< y < yNy +

∆y

2
,

z1 +
∆z

2
< z < zNz +

∆z

2

}
. (3.30)

O domínio supercondutor, de lados (a, b, c) é limitado por faces cada qual a uma

distância correspondente a um múltiplo inteiro de espaçamentos da malha:

Ωsc =

{
r = (x, y, z) : xnix +

∆x

2
< x < xnfx +

∆x

2
, yniy +

∆y

2
< y < ynfy +

∆y

2
,

zniz +
∆z

2
< z < znfz +

∆z

2

}
, (3.31)

onde, por exemplo, nix = (A − a)/2 + 1 e nfx = (A + a)/2 + 1, uma vez que os

paralelepípedos são concêntricos. Os índices (ix, iy, iz) e (fx, fy, fz) indicam os pontos

iniciais e finais nas direções (x, y, z), respectivamente.

As condições de contorno (3.25-3.26) são aplicadas em ∂Ωsc e as (3.28) em ∂Ω.

3.8.2 Definições

Antes de discretizarmos as equações TDGL, precisamos introduzir algumas definições.

A principal delas é a variável de ligação. Consideremos o produto do campo auxiliar Ux

entre dois pontos consecutivos de uma malha (xi, yj, zk) e (xi+1, yj, zk) ao longo do eixo



O Método ψ − U 53

Figura 3.2: Malha de discretização: a legenda indica os pontos onde cada quantidade
física é calculada. Fonte: Próprio autor.

x:

Ux,i,j,k ≡ Ūx(xi, yj, zk)Ux(xi+1, yj, zk)

= exp

(
−i
∫ xi+1

xi

Axdx

)
= exp (−iAx,i,j,k∆x) , (3.32)

onde usamos a regra do ponto médio para integração. Analogamente, temos que:

Uy,i,j,k ≡ Ūy(xi, yj, zk)Uy(xi, yj+1, zk)

= exp (−iAy,i,j,k∆y) , (3.33)

Uz,i,j,k ≡ Ūz(xi, yj, zk)Uy(xi, yj, zk+1)

= exp (−iAz,i,j,k∆z) , (3.34)



O Método ψ − U 54

ondeAx,i,j,k ≡ Ax(xi+
∆x
2
, yj, zk),Ay,i,j,k ≡ Ay(xi, yj+

∆y
2
, zk) eAz,i,j,k ≡ Az(xi, yj, zk+

∆z
2

). Assim, fica claro a denominação "variável de ligação", pois as variáveis Ux,i,j,k,

Uy,i,j,k e Uz,i,j,k ligam pontos adjacentes da malha através do potencial no ponto médio.

Outras definições são as seguintes:

• ψi,j,k ≡ ψ(xi, yj, zk),

• ϕi,j,k ≡ ϕ(xi, yj, zk),

• Jsx,i,j,k ≡ Jsx(xi + ∆x
2
, yj, zk),

• Jsy,i,j,k ≡ Jsy(xi, yj + ∆y
2
, zk),

• Jsz,i,j,k ≡ Jsy(xi, yj, zk + ∆z
2

),

• Bx,i,j,k ≡ Bz(xi, yj + ∆y
2
, zk + ∆z

2
),

• By,i,j,k ≡ Bz(xi + ∆x
2
, yj, zk + ∆z

2
),

• Bz,i,j,k ≡ Bz(xi + ∆x
2
, yj + ∆y

2
, zk).

A Figura 3.2 ilustra os pontos onde cada quantidade física é calculada.

3.8.3 Aproximações para as Derivadas

Usando expansão em série de Taylor, obtemos:

f(x+
∆x

2
) = f(x) +

∆x

2
f ′(x) +

∆x2

8
f ′′(x) + · · · ,

f(x− ∆x

2
) = f(x)− ∆x

2
f ′(x) +

∆x2

8
f ′′(x) + · · · .

Subtraindo uma da outra, obtemos a aproximação:

f ′(x) =
f(x+ ∆x

2
)− f(x− ∆x

2
)

∆x
. (3.35)

A interpretação geométrica desta fórmula encontra-se na figura (3.3).



O Método ψ − U 55

Figura 3.3: A figura ilustra a aproximação (3.35) para a primeira derivada [3].

Usando a aproximação (3.35) obtemos:

f ′′(x) =
f ′(x+ ∆x

2
)− f ′(x− ∆x

2
)

∆x

=
f(x+∆x)−f(x)

∆x
− f(x)−f(x−∆x)

∆x

∆x

=
f(x+ ∆x)− 2f(x) + f(x−∆x)

∆x2
. (3.36)

3.8.4 Discretização da Primeira Equação TDGL

Usando a aproximação para a segunda derivada (3.36) e as de variáveis de ligação (3.32-

3.33) encontramos:

Ūx
∂2

∂x2
(Uxψ)

∣∣∣∣
(xi,yj ,zk)

= Ūx(xi, yj, zk) [Ux(xi+1, yj, zk)ψ(xi+1, yj, zk)

−2Ux(xi, yj, zk)ψ(xi, yj, zk)

+Ux(xi−1, yj, zk)ψ(xi−1, yj, zk)] /∆x
2

=
Ux,i,j,kψi+1,j,k − 2ψi,j,k + Ūx,i−1,j,kψi−1,j,k

∆x2
. (3.37)

De forma análoga, podemos obter uma aproximação para a segunda derivada com relação

às coordenadas y e z fazendo uma permutação cíclica dos índices (i, j, k).

Consequentemente, a primeira equação TDGL (3.22) pode ser aproximada por:

∂ψi,j,k
∂t

= Di,j,k , (3.38)



O Método ψ − U 56

onde,

Di,j,k =
Ux,i,j,kψi+1,j,k − 2ψi,j,k + Ūx,i−1,j,kψi−1,j,k

∆x2

+
Uy,i,j,kψi,j+1,k − 2ψi,j,k + Ūy,i,j−1,kψi,j−1,k

∆y2

+
Uz,i,j,kψi,j,k+1 − 2ψi,j,k + Ūz,i,j,k−1ψi,j,k−1

∆z2

−iϕi,j,kψi,j,k + (1− T )ψi,j,k(1− |ψi,j,k|2) , (3.39)

para todo {i = nix + 1, . . . , nfx− 1, j = niy + 1, . . . , nfy − 1, k = niz + 1, . . . , nfz − 1}.

3.8.5 Discretização da Segunda Equação TDGL (Lei de Ampère)

Para finalizar, precisamos discretizar a segunda equação TDGL (3.23). Temos que:

β

(
∂Ax
∂t

+
∂ϕ

∂x

)
= τJsx − κ2

(
∂Bz

∂y
− ∂By

∂z

)
, (3.40)

β

(
∂Ay
∂t

+
∂ϕ

∂y

)
= τJsy − κ2

(
∂Bx

∂z
− ∂Bz

∂x

)
, (3.41)

β

(
∂Az
∂t

+
∂ϕ

∂z

)
= τJsz − κ2

(
∂By

∂x
− ∂Bx

∂y

)
, (3.42)

onde τ(x, y, z) é uma função que vale um em Ωsc e zero em Ω\Ω̄sc.

Agora, usando a aproximação (3.35) para a primeira derivada, encontramos:

∂Bz

∂y

∣∣∣∣
(xi+

∆x
2
,yj ,zk)

=
Bz(xi + ∆x

2
, yj + ∆y

2
, zk)−Bz(xi + ∆x

2
, yj − ∆y

2
, zk)

∆y

=
Bz,i,j,k −Bz,i,j−1,k

∆y
, (3.43)

∂By

∂z

∣∣∣∣
(xi+

∆x
2
,yj ,zk)

=
By(xi + ∆x

2
, yj, zk + ∆z

2
)−By(xi + ∆x

2
, yj, zk − ∆z

2
)

∆z

=
By,i,j,k −By,i,j,k−1

∆z
, (3.44)

∂ϕ

∂x

∣∣∣∣
(xi+

∆x
2
,yj ,zk)

=
ϕ(xi + ∆x, yj, zk)− ϕ(xi, yj, zk)

∆x

=
ϕi+1,j,k − ϕi,j,k

∆x
. (3.45)



O Método ψ − U 57

Substituindo estes resultados em (3.40) encontramos:

∂Ax,i,j,k
∂t

=
1

β

[
τx,i,j,kJsx,i,j,k −

κ2

∆y
(Bz,i,j,k −Bz,i,j−1,k) +

κ2

∆z
(By,i,j,k −By,i,j,k−1)

]
−ϕi+1,j,k − ϕi,j,k

∆x
, (3.46)

onde τx,i,j,k = τ(xi + ∆x
2
, yj, zk), para todo {i = 1, 2, . . . , Nx, j = 2, 3, . . . , Ny, k =

2, 3, . . . , Nz}.

Podemos encontrar aproximações para as outras componentes do potencial vetor por

meio de permutação cíclica dos índices (i, j, k) e das coordenadas (x, y, z), de modo que:

∂Ay,i,j,k
∂t

=
1

β

[
τy,i,j,kJsy,i,j,k −

κ2

∆z
(Bx,i,j,k −Bx,i,j,k−1) +

κ2

∆x
(Bz,i,j,k −Bz,i−1,j,k)

]
−ϕi,j+1,k − ϕi,j,k

∆y
, (3.47)

onde τy,i,j,k = τ(xi, yj + ∆y
2
, zk), para todo {i = 2, 3, . . . , Nx, j = 1, 2, . . . , Ny, k =

2, 3, . . . , Nz}, e,

∂Az,i,j,k
∂t

=
1

β

[
τz,i,j,kJsz,i,j,k −

κ2

∆x
(By,i,j,k −By,i−1,j,k) +

κ2

∆y
(Bx,i,j,k −Bx,i,j−1,k)

]
−ϕi,j,k+1 − ϕi,j,k

∆z
, (3.48)

onde τz,i,j,k = τ(xi, yj, zk + ∆z
2

), para todo {i = 2, 3, . . . , Nx, j = 2, 3, . . . , Ny, k =

1, 2, . . . , Nz}.

3.8.6 Densidade de Corrente

A componente x da densidade de corrente (3.20) envolve o produto dos termos Ūxψ̄ e

∂Uxψ
∂x

. Usando a aproximação para a primeira derivada temos que:

∂(Uxψ)

∂x

∣∣∣∣
(xi+

∆x
2
,yj ,zk)

=
Ux(xi+1, yj, zk)ψ(xi+1, yj, zk)− Ux(xi−1, yj, zk)ψ(xi−1, yj, zk)

∆x
.

(3.49)

O outro termo aproximamos pelo seu valor médio entre os pontos adjacentes na malha:

Ūxψ̄
∣∣
(xi+

∆x
2
,yj ,zk)

=
Ūx(xi+1, yj, zk)ψ̄(xi+1, yj, zk) + Ūx(xi−1, yj, zk)ψ̄(xi−1, yj, zk)

2
.

(3.50)



O Método ψ − U 58

Multiplicando as equações (3.49) e (3.50) pela identidade

Ux(xi, yj, zk)Ūx(xi, yj, zk) = 1 e usando a definição de variável de ligação (3.32), ob-

temos:

Ūxψ̄
∂(Uxψ)

∂x

∣∣∣∣
(xi+

∆x
2
,yj ,zk)

=
1

2∆x

[
Ūx,i,j,kψ̄i+1,j,k + ψ̄i,j,k

]
[Ux,i,j,kψi+1,j,k − ψi,j,k] .

(3.51)

Usando a equação (3.20) e a aproximação obtida com a equação (3.51), a densidade

de corrente pode ser escrita como:

Jsx,i,j,k =
1

∆x
Im
[
ψ̄i,j,kUx,i,j,kψi+1,j,k

]
, (3.52)

para todo {i = nix, . . . , nfx − 1, j = niy + 1, . . . , nfy − 1, k = niz + 1, . . . , nfz − 1}.

Para as componentes y e z temos que:

Jsy,i,j,k =
1

∆y
Im
[
ψ̄i,j,kUy,i,j,kψi,j+1,k

]
, (3.53)

para todo {i = nix + 1, . . . , nfx − 1, j = niy, . . . , nfy − 1, k = niz + 1, . . . , nfz − 1}, e,

Jsz,i,j,k =
1

∆z
Im
[
ψ̄i,j,kUz,i,j,kψi,j,k+1

]
, (3.54)

para todo {i = nix + 1, . . . , nfx − 1, j = niy + 1, . . . , nfy − 1, k = niz, . . . , nfz − 1}.

3.8.7 Campo Local

Usando a definição:

B = Bxi +Byj +Bzk

=

(
∂Az
∂y
− ∂Ay

∂z

)
i +

(
∂Ax
∂z
− ∂Az

∂x

)
y +

(
∂Ay
∂x
− ∂Ax

∂y

)
k , (3.55)

e a aproximação para a primeira derivada (3.35) encontramos:

Bx,i,j,k =
Az,i,j+1,k − Az,i,j,k

∆y
− Ay,i,j,k+1 − Ay,i,j,k

∆z
,

By,i,j,k =
Ax,i,j,k+1 − Ax,i,j,k

∆z
− Az,i+1,j,k − Az,i,j,k

∆x
,

Bz,i,j,k =
Ay,i+1,j,k − Ay,i,j,k

∆x
− Ax,i,j+1,k − Ax,i,j,k

∆y
, (3.56)



O Método ψ − U 59

para todo {i = 2, 3, . . . , Nx − 1, j = 2, 3, . . . , Ny − 1, k = 2, 3, . . . , Nz − 1}.

3.8.8 Condições de Contorno

Supondo continuidade do campo na interface ∂Ω podemos escrever para a continuidade

do campo local nas faces da esquerda e da direita, respectivamente:

Bz,1,j,k = Hz , (3.57)

Bz,Nx,j,k = Hz , (3.58)

para todo {j = 1, 2, . . . , Ny, k = 1, 2, . . . , Nz + 1}; para as faces de frente e de fundo,

respectivamente:

Bz,i,1,k = Hz , (3.59)

Bz,i,Ny,k = Hz , (3.60)

para todo {i = 1, 2, . . . , Nx, k = 1, 2, . . . , Nz + 1}; para as faces inferior e superior,

respectivamente:

Bz,i,j,1 = Hz , (3.61)

Bz,i,j,Nz+1 = Hz , (3.62)

para todo {i = 1, 2, . . . , Nx, j = 1, 2, . . . , Ny}. As outras componentes do campo local

podem ser obtidas por permutação cíclica dos índices (i, j, k) e das coordenadas (x, y, z).

Usando novamente a aproximação para a primeira derivada (3.35), e a definição de

variáveis de ligação (3.32-3.34), as expressões (3.25-3.26) podem ser aproximadas por:

ψ1,j,k = Ux,1,j,kψ2,j,k , (3.63)

ψNx+1,j,k = Ūx,Nx,j,kψNx,j,k , (3.64)

para as faces da esquerda e da direita respectivamente, para todo



O Método ψ − U 60

j = niy + 1, . . . , nfy − 1, k = niz + 1, . . . , nfz − 1,

ψi,1,k = Uy,i,1,kψi,2,k , (3.65)

ψi,Ny+1 = Ūy,i,Ny ,kψi,Ny ,k , (3.66)

para as faces da frente e do fundo respectivamente, para todo

i = nix + 1, . . . , nfx − 1, k = niz + 1, . . . , nfz − 1,

ψi,j,1 = Uz,i,j,1ψi,j,2 , (3.67)

ψi,j,Nz+1 = Ūz,i,j,Nzψi,j,Nz , (3.68)

para as faces inferior e superior respectivamente, para todo

i = nix + 1, . . . , nfx − 1, j = niy + 1, . . . , nfy − 1.

Conhecendo-se o parâmetro de ordem e as variáveis de ligação no interior da malha,

determinamos os valores de ψ nas faces mais próximas a ∂Ωsc e externas ao domínio Ωsc.

Neste ponto podemos ver que os valores do parâmetro de ordem não dependem de

seus valores fora do domínio Ωsc. Por exemplo, substituindo (3.63) no primeiro termo da

equação (3.39), e considerando i = 1, encontramos (Ux,1,j,kψ2,j,k − ψ2,j,k)/∆x
2. Estes

valores de ψ são usados apenas para efeitos de ilustração do parâmetro de ordem.

3.9 Evolução Temporal

Para determinar a evolução temporal das quantidades físicas, usamos o método de Euler

explícito. Considerando a aproximação de primeira ordem [3]:∫ t+∆t

t

f(t)dt = f(t)∆t , (3.69)

as equações (3.38) e (3.46-3.48) podem ser integradas e obtemos os resultados:



O Método ψ − U 61

ψi,j,k(t+ ∆t) = ψi,j,k(t) + ∆tDi,j,k(t) , (3.70)

Ax,i,j,k(t+ ∆t) = Ax,i,j,k(t) +
∆t

β

{
τx,i,j,kJsx,i,j,k(t)−

κ2

∆y
[Bz,i,j,k(t)−Bz,i,j−1,k(t)]

+
κ2

∆z
[By,i,j,k(t)−By,i,j,k−1(t)]

}
−∆t

[ϕi+1,j,k(t)− ϕi,j,k(t)]
∆x

, (3.71)

Ay,i,j,k(t+ ∆t) = Ay,i,j,k(t) +
∆t

β

{
τy,i,j,kJsy,i,j,k(t)−

κ2

∆z
[Bx,i,j,k(t)−Bx,i,j,k−1(t)]

+
κ2

∆x
[Bz,i,j,k(t)−Bz,i−1,j,k(t)]

}
−∆t

[ϕi,j+1,k(t)− ϕi,j,k(t)]
∆y

, (3.72)

Az,i,j,k(t+ ∆t) = Az,i,j,k(t) +
∆t

β

{
τz,i,j,kJsz,i,j,k(t)−

κ2

∆x
[By,i,j,k(t)−By,i−1,j,k(t)]

+
κ2

∆y
[Bx,i,j,k(t)−Bx,i,j−1,k(t)]

}
−∆t

[ϕi,j,k+1(t)− ϕi,j,k(t)]
∆z

, (3.73)

3.10 Algoritmo

A implementação do algorítimo funciona da seguinte forma: suponhamos que conhece-

mos todas as quantidades em todos os pontos da malha em um determinado instante de

tempo t. Então, podemos determiná-las em t + ∆t. Em outras palavras, se conhecemos

(ψi,j,k, Ax,i,j,k, Ay,i,j,k, Az,i,j,k) em t, das equações (3.39), (3.52-3.54), e (3.56-3.56), po-

demos calcular todas as quantidades necessárias para determinar o lado direito de (3.70-

3.73) nos pontos interiores da malha. Finalmente, as condições de contorno são aplicadas.

Antes de ir para o próximo tempo, atualizamos todas as quantidades relevantes como as

variáveis de ligação, densidade de corrente e o campo local em todos os pontos da ma-

lha. Então vamos para o nível de tempo acima. Repetindo este procedimento podemos

determinar todas as quantidades físicas relevantes em todos os níveis de tempo até que



O Método ψ − U 62

o sistema alcance o estado estacionário. Infelizmente, este método não converge para

qualquer passo de tempo ∆t. Usamos a seguinte regra prática [3]:

∆t ≤ min

{
δ2

4
,
δ2β

4κ2

}
, (3.74)

δ2 =
2

1
∆x2 + 1

∆y2 + 1
∆z2

. (3.75)

(Ver Apêndice A).

O algorítimo descrito até este ponto pressupõe que o potencial esteja fixo, ou seja, é

conhecido. Com efeito, podemos determiná-lo em todos os pontos da malha. Do contrá-

rio, teríamos de recorrer a alguma outra lei física para determiná-lo. Uma vez que iremos

considerar o supercondutor apenas na presença de campo magnético aplicado, podemos

trabalhar no calibre de Coulomb ϕ = 0.



4
Métodos de Paralelização do Algoritmo ψ − U

4.1 Conceitos Básicos Sobre Paralelização com Tecnolo-
gia CUDA R© e Linguagem C

4.1.1 Introdução

A descrição das teorias deste capítulo tiveram como base as videoaulas: “Introduction to

Parallel Programming” do site www.udacity.com e o livro de Sanders e Kandrot[25].

Na Universidade Estadual Paulista Júlio de Mesquita Filho, Campus de Bauru, o mé-

todo ψ − U foi originalmente implementado em linguagem Fortran pelo professor André

Luiz Malvezzi e posteriormente aprimorado pelo professor Edson Sardella. Neste traba-

lho, traduzimos este código para a linguagem C juntamente com técnicas de paralelização.

Antes de descrever como paralelizamos o algoritmo ψ−U , iremos descrever a tecnologia

CUDA.

O significado da palavra “paralelização” em computação está relacionado com a pos-

sibilidade de que trechos de código de um programa sejam executados por vários proces-

sadores. De fato, uma das possibilidades de se alcançar tal feito é direcionando o trecho

de código para vários computadores uniprocessados em uma rede denominada cluster,

com utilização de MPI (Message Passing Interface) - um sistema de transferência de pro-

cessamentos entre os computadores conectados aos nós do cluster - a fim de se otimizar o

63



Métodos de Paralelização do Algoritmo ψ − U 64

tempo de execução (caso de códigos com tempo de processamento e custo computacional

altos). Uma outra possibilidade de implementação de paralelização surgiu por meio da

empresa NVIDIA R©, famosa por fabricar placas de vídeo cuja GPU (Graphics Processing

Unit - Unidade de Processamento Gráfico), quando utilizada em conjunto com a tecno-

logia CUDA R© (Compute Unified Device Architecture), é capaz de comportar a tarefa de

paralelização semelhantemente ao cluster mencionado. As principais diferenças em se

trabalhar com esses dois tipos de sistemas, cluster com MPI e CUDA, podem ser enu-

meradas por suas vantagens e desvantagens. As vantagens de se trabalhar com CUDA

são: a diminuição do espaço físico e de equipamentos, uma vez que, via de regra, cluster

com MPI necessita mais de um computador para execução dos processamentos em pa-

ralelo, enquanto CUDA pode ser implementada utilizando-se apenas uma máquina com

o hardware da NVIDIA; substituição do intercâmbio de processamentos via MPI, o qual

enfrenta um gargalo - comumente denominado bottleneck - pelo intercâmbio de proces-

samentos via CUDA, já que a maioria das topologias e protocolos de rede apresentam

uma velocidade de transferência de dados inferior à do barramento de comunicação entre

memória e processador; a utilização de memória compartilhada ao se utilizar CUDA é

um método de compartilhamento de dados mais eficaz que a utilização de memória dis-

tribuída ao se utilizar cluster com MPI, pois neste último caso, um módulo de memória

pode ser acessado diretamente por apenas um dos nós do cluster, ou seja, não há acesso

simultâneo. Já as principais desvantagens de se trabalhar com CUDA em vez de clus-

ter com MPI seriam: menor escalabilidade, significando que, dado um aumento da carga

de processamentos em um cluster, é possível adicionar outros nós para subdivisão dos

processos, enquanto na paralelização com CUDA - quando fazemos referência à utiliza-

ção de uma única máquina - a escalabilidade depende da quantidade de slots livres nesta

máquina. Devemos considerar, ainda, que existe a possibilidade de se trabalhar com am-

bas tecnologias, CUDA e cluster com MPI, caso haja várias máquinas com software e



Métodos de Paralelização do Algoritmo ψ − U 65

hardware necessários instalados.

Para que a paralelização possa ser implementada nos moldes propostos neste tralha-

lho, é necessária, além da GPU, também da linguagem CUDA, responsável por modelar

a paralelização, que também será responsável pelo intercâmbio da comunicação entre

os conjuntos CPU + memória principal (referenciado como host) e GPU + memória da

placa de vídeo (referenciado como device), conforme ilustrado pela Figura 4.1. Esse

conjunto forma a plataforma de computação paralela da NVIDIA. Outros itens necessá-

rios, como: compilador, drivers específicos para comunicação entre o device e o sistema

operacional utilizado, e o pacote CUDA development toolkit - que contém as bibliotecas

de modelagem dos problemas paralelizáveis - podem ser obtidos no site da fabricante

www.nvidia.com gratuitamente.

Figura 4.1: Hardware básico para computação paralela com CUDA: Com a implementa-
ção da placa de vídeo (device), o programa, antes escrito para ser executado apenas no
host ((a) Memória principal e (b) CPU), agora também pode ser executado no device ((c)
Memória e (d) GPU). Fonte: Próprio autor.

Diferentemente de uma CPU (Central Processing Unit - Unidade Central de Proces-

samento) composta geralmente por um núcleo com ULA (Unidade Lógica Aritmética)



Métodos de Paralelização do Algoritmo ψ − U 66

e uma unidade de controle, as GPUs da NVIDIA possuem várias ULAs denominadas

CUDA cores, que são os núcleos CUDA para onde são direcionadas as partes do algo-

ritmo consideradas críticas em termos de tempo de processamento (geralmente devido

ao grande volume de dados e/ou laços de repetição multi dimensionados contendo fun-

ções de tamanho e complexidade consideráveis para cálculo, o que é o caso da simulação

computacional de um material supercondutor com o método ψ − U implementado).

A Figura 4.2 resume os passos para implementação da paralelização.

Figura 4.2: Estrutura básica de processamento em paralelo. Fonte: Próprio autor.

Como exemplo de código a ser processado pela GPU, um laço “for” de C/C++ com

índices i, j, k, pode ser subdividido em partes denominadas threads, cada qual sendo

executada por uma ULA (core) da GPU. O esquema simplificado de funcionamento con-

siste no seguinte: Um algoritmo sendo executado serialmente no host, em dado momento

executa a alocação de memória no device e transfere os dados necessários para o processa-

mento paralelo por meio da linguagem CUDA - que consegue acessar ambas as memórias

e intercambiar o controle do processamento. Na sequência, ocorrem as chamadas de fun-

ções que serão executadas no device, seguida da passagem de controle para a GPU, que se

encarregará do processamento. Ao término dos processamentos e atualização da memó-



Métodos de Paralelização do Algoritmo ψ − U 67

ria do device, ocorre a transferência para a memória principal e a passagem do controle

do processamento da GPU para a CPU. Cabe salientar, ainda, que enquanto não ocorre

a transferência do controle da GPU para a CPU, esta última não retoma o controle e não

executa outros trechos de código do restante do programa.

Para que o item 2 da Figura 4.2 seja implementado, é necessário estabelecer um con-

ceito fundamental para se trabalhar com paralelização e relacionado com o acesso das

regiões de memória, denominado flatten (achatamento). Antes, vamos analizar como a

memória é alocada em C/C++ e, posteriormente, com CUDA. Tomando como exemplo

um vetor 3D (aqui denominado “vetor3D”), o acesso das posições de memória deste vetor

pode ser realizado pelos usuais índices i, j, k na forma: vetor3D[i][j][k]. Este vetor 3D

poderia ser alocado dinamicamente conforme Figura 4.3 utilizando uma função malloc()

em conjunto com os laços de repetição.

Assim, pela Figura 4.3, vemos que em (a) uma região de memória guarda como dado

um endereço de memória (por isso recebe o nome de ponteiro). A região de memória

apontada por esse ponteiro é o primeiro bloco de (b). Caso uma função do tipo malloc()

seja utilizada, semelhantemente ao que ocorre na linha 2 do algoritmo desta figura, me-

mória é alocada a partir do primeiro bloco apontado por (a) e a quantidade total desta

alocação (largura*sizeof(double)) é a contida na chave indicada por 2 na figura. Como

todas essas posições de memória alocada são também do tipo ponteiro, estas têm igual-

mente a capacidade de apontar para outras regiões de memória e, quando utilizadas a

função malloc() em conjunto com um laço de repetição (conforme linhas 3 a 5 do algo-

ritmo) esses ponteiros também podem alocar mais regiões de memória, e assim sucessi-

vamente, até que os últimos ponteiros apontem para regiões de memória cuja alocação

servirá para armazenamento dos dados do programa em (d). Todo esse processo leva à

seguinte conclusão: vetores multi dimensionados em C/C++ na verdade são conjuntos de

vetores unidimensionais. Por fim, o acesso a um dado dessa região que se encontre, por



Métodos de Paralelização do Algoritmo ψ − U 68

Figura 4.3: Esquema simplificado de alocação dinâmica de memória em C/C++. Os
cubos em (a), (b) e (c) representam regiões de memória que armazenam o endereço de
outras regiões de memória (por isso são denominados de ponteiros). Em (d) encontra-se
o espaço alocado para ser utilizado pelos dados processados. Fonte: Próprio autor.

exemplo, na posição [4][6][4] de vetor3D, na verdade é um comando para que este vetor

acesse o ponteiro 4 de (b), ponteiro 6 de (c) e região 4 da memória alocada em (d).

Cabe salientar que também é possível a alocação de memória em C/C++ utilizando

apenas um comando de criação de variável, do tipo:

double teste[quantidade], (4.1)

mas a explicação do procedimento com utilização de ponteiros, além de ilustrar como esse

processo é executado em nível mais baixo nesta linguagem, esclarece melhor a alocação

unidimensional em CUDA, e as linhas de execução - threads - também têm semelhança



Métodos de Paralelização do Algoritmo ψ − U 69

com a forma como a memória é alocada tridimensionalmente. Além disso, para os que

trabalham com a modelagem de problemas e estrutura de dados, esses conceitos são de

suma importância.

Em CUDA, a alocação dinâmica de memória no device - utilizando a função “cuda-

Malloc()” - deve ser realizada em uma única etapa, pois esta função não permite multi

dimensionamento da memória no device (como descrito na abordagem envolvendo mal-

loc, memória principal e laços de repetição). Assim, devemos utilizar o conceito de acha-

tamento de índice a fim de trabalharmos com alocação e transferência de dados entre

device e host. A Figura 4.4 exemplifica o achatamento de índice, enquanto as linhas de

programação da Figura 4.5 são um exemplo de alocação de memória no device.

Figura 4.4: Achatamento de índice. Em (a), temos o vetor tridimensional antigo, e em
(b), o novo vetor unidimensional, ambos com quantidades idênticas de memória alocada.
Após esse processo, o acesso às posições do novo vetor em (b) - o qual está relacionado
com as posições do vetor tridimensional em (a) - é feito por meio de um índice unidimen-
sional. Fonte: Próprio autor.

Pela Figura 4.4, as dimensões do paralelepípedo em (a) são dadas por X (largura), Y

(profundidade) e Z (altura). Como (b) é unidimensional mas possui a mesma quantidade

de posições de memória que (a), sua dimensão tem tamanho X ∗ Y ∗ Z. Então, para

acessar cada uma das posições de memória do vetor em (b), devemos usar a relação:

i + j ∗X + k ∗X ∗ Y ou i + X ∗ (j + k ∗ Y ) . (4.2)



Métodos de Paralelização do Algoritmo ψ − U 70

Figura 4.5: A linha 1 cria um ponteiro para memória do tipo específico que se quer tra-
balhar, enquanto na linha 2, além da quantidade e tamanho desses espaços de memória, o
comando retorna um ponteiro da memória alocada no device para vetor3D. Fonte: Próprio
autor.

Como exemplo, considerando que o paralelepípedo (a) da Figura 4.4 tem dimensões

X = 5, Y = 7 e Z = 5, se todas as posições de (a) forem copiadas na sequência, desde

[0][0][0] até [4][6][4] para (b), o índice equivalente para, por exemplo, a antiga posição

[4][6][4], seria: vetor3D[4 + 6*5 + 4*5*7] = vetor3D[174] neste novo vetor.

Após estabelecido o “achatamento” dos índices dos vetores, é possível iniciar o pro-

cedimento de cópia dos dados do host para o device. A função utilizada para isso é dada

por:

cudaMemcpy(campo 1, campo 2, campo 3, campo 4) , (4.3)

a qual efetua a transferência de dados entre host e device, e atua da seguinte forma: campo

1 é o nome do ponteiro para região de memória destino, campo 2 é o nome do ponteiro

para região de memória fonte, campo 3 é o tamanho da transferência e campo 4 é o

comando que indica o sentido do fluxo desses dados. A Figura 4.6 exemplifica as trans-

ferências: (a) do Host para o Device e (b) do Device para o Host.

Figura 4.6: (a) Transferência HostToDevice. (b) Transferência DeviceToHost. Fonte:
Próprio autor.

Após a transferência, o device está preparado para a execução do trecho crítico do

programa. Este trecho - geralmente um laço de repetição com tempo de execução crítico -



Métodos de Paralelização do Algoritmo ψ − U 71

deve ser transformado em uma função - referenciada como kernel - cujos parâmetros são

os vetores a serem passados por referência. Na Figura 4.7, temos um exemplo de trecho

de código que é executado apenas em série, enquanto a Figura 4.9 representa o mesmo

código paralelizado.

Figura 4.7: Código de execução serial. Fonte: Próprio autor.

A seguir, vamos tratar de outros conceitos relacionados às linhas de processo paralelo,

as threads, blocks e grids.

4.1.2 Threads, Blocks e Grids

Os threads blocks e grids são elementos utilizados para subdivisão dos processos a serem

paralelizados. Esses termos têm relação com as capacidades de processamento das GPUs

e também com a maneira como o processamento é executado, se uni, bi ou tridimensi-

onalmente (na verdade, os problemas a serem paralelizados podem ter até mais que três

dimensões, dependendo das especificações das GPUs utilizadas). Na Seção 4.2, empre-



Métodos de Paralelização do Algoritmo ψ − U 72

garemos esses conceitos na simulação de um material supercondutor levando em conta as

três dimensões.

As threads constituem o elemento mais básico da computação paralela. Para processos

grandes e que ultrapassam a capacidade das GPUs de execução simultânea de threads,

estas podem ser subdivididas em conjuntos denominados blocks e grids - blocos e grades

- que também possuem limites, dependendo das especificações da GPU. Para processos

menores, mas cujo problema possui mais de uma dimensão, tais recursos de subdivisão

também são importantes, pois é possível relacionar cada dimensão do problema a ser

tratado com as dimensões dos blocos de threads e da grade de blocos.

A quantidade de processos que podem ocorrer em paralelo depende das especificações

da GPU, e uma delas é a capacidade computacional. Esta capacidade está relacionada à

quantidade de dimensões qu