Miguel Dias Costa
Propriedades Magnéticas de
Sistemas de Partículas Acopladas
por Interacção Magnetostática
Dissertação orientada pelo Prof. Dr. Yuriy Pogorelov e pelo Prof. Dr. João Viana
Lopes, submetida à Faculdade de Ciências da Universidade do Porto para a obtenção do
grau de Doutor em Física
Departamento de Física
Faculdade de Ciências da Universidade do Porto
Junho 2010
Agradecimentos
Ao longo de todos estes anos, foram muitas as pessoas que me apoiaram e criaram
condições para que eu pudesse levar a cabo este desafio. Por tudo o que fizeram por
mim, estou grato à minha família, aos meus orientadores e professores e a todos os amigos
na Universidade do Porto (dos Departamentos de Física e de Química da Faculdade de
Ciências, na Escola de Gestão do Porto e na Faculdade de Economia), na Escola Superior
de Tecnologia a Gestão do Instituto Politécnico de Viana do Castelo, na GMS Consulting
e na Neoscopio.
Esta dissertação foi apoiada pela Fundação para a Ciência e Tecnologia, através da Bolsa
de Doutoramento SFRH/BD/7003/2001, com o co-financiamento do Programa Operacional Ciência e Tecnologia 2010 e do Fundo Social Europeu.
3
Agradecimentos
4
Resumo
Este trabalho dedica-se ao estudo, através de métodos tanto analíticos como numéricos,
de propriedades termodinâmicas de sistemas planares constituídos por partículas magnéticas nanoscópicas acopladas pela interacção magnetostática entre os seus momentos
magnéticos. Do ponto de vista prático, estes sistemas servem como modelos de modernos
dispositivos nano-estruturados. Do ponto de vista teórico, apresentam comportamentos
interessantes, quando comparados com sistemas mais tradicionais de física estatística,
devido ao carácter de longo alcance e anisotrópico das interacções.
Na introdução são apresentados alguns desenvolvimentos experimentais e tecnológicos
que continuam a sustentar o interesse nestes sistemas, bem como os modelos utilizados
para os estudar. Segue-se uma exposição de métodos utilizados durante este trabalho.
No capítulo 2 é apresentado um método de expansão de altas temperaturas que, pela sua
generalidade, pode ser usado como mecanismo de controlo de métodos numéricos como,
no caso particular, simulações de Monte Carlo. Estes métodos numéricos são apresentados
no capítulo 3, bem como os problemas de convergência em sistemas magnéticos frustrados
e formas de ultrapassar estes problemas. É também incluído um levantamento de algumas
ferramentas de software utilizadas para a realização deste trabalho, na esperança de vir
a ser útil para quem se esteja a iniciar neste ramo.
No capítulo 4, os métodos apresentados anteriormente são aplicados a redes regulares
de dipolos pontuais. A medição dos tempos de autocorrelação na energia e no parâmetro
de ordem de uma rede quadrada revelou não só o problema típico de critical slowing down
mas também um abrandamento abaixo da temperatura crítica relacionado com aspectos
particulares da interacção dipolar. Através da utilização de técnicas como o método de
Parallel Tempering e optimizações da taxa de aceitação das simulações, estes tempos de
autocorrelação foram reduzidos significativamente relativamente às abordagens canónicas,
o que permitiu melhorar a precisão das simulações e caracterizar melhor o comportamento
crítico destes sistemas. No capítulo 5, alguns destes métodos foram aplicados a um sistema
desordenado de dipolos pontuais, com o objectivo de clarificar o efeito da anisotropia da
interacção no comportamento das correlações magnéticas entre partículas.
No capítulo 6, é abandonada a aproximação de dipolos pontuais e são calculadas correcções de tamanho finito para diferentes formas das partículas. Estes resultados são utiliza-
5
Resumo
dos para explicar a anisotropia encontrada em experiências de ressonância ferromagnética
numa rede quadrada de discos magnéticos, quando é variada a orientação, relativamente
aos eixos da rede, do campo aplicado. É demonstrado que a inclusão destas correcções de
tamanho finito, mantendo a magnetização uniforme e considerando as interacções entre
partículas, leva a uma distribuição de campos internos que exibe a anisotropia verificada
experimentalmente. É também apresentado um procedimento variacional que permite
estimar a deformação da magnetização no interior da partícula devida tanto à forma das
partículas como à interacção entre partículas.
Finalmente, no capítulo 7, são apresentadas as conclusões e questões em aberto relacionadas com este trabalho.
6
Abstract
This works focuses on the analytical and numerical study of thermodynamic properties
of planar systems with nanoscopic magnetic particles coupled by the magnetostatic interaction between their magnetic moments. From a practical point of view, these systems
act as models for modern nanostructured devices. From a theoretical point of view, they
present interesting behaviours, when compared to more traditional statistical physics systems, due to the long-range and anisotropic character of the interactions.
In the introduction, some relevant experimental and technological developments are
presented, along with the models used to study these systems. This is followed by an
exposition of methods that were used during this work. In Chapter 2, a method of high
temperature series expansion is presented, which, thanks to its generality, can be used
as a control mechanism for numerical methods such as, in our present case, Monte Carlo
simulations. These numerical methods are presented in Chapter 3, as well as the convergence problems in frustrated magnetic systems and ways of overcoming those problems.
A survey of some software tools used during this work is also included, with the hope that
it can be useful to someone coming into computational physics.
In Chapter 4, these methods are applied to studying regular lattices of point dipoles.
Measuring the autocorrelation times of the energy and order parameter revealed not only
the typical problem of critical slowing down but also an increase of these times below
the critical temperature, which is related to relevant specific aspects of the magnetostatic interaction. Through the use of techniques like the method of Parallel Tempering
and optimizations of the acceptance rate in simulations, these autocorrelation times were
significantly reduced when compared with canonical approaches, which allowed for an increased accuracy of the simulations and a better characterization of the critical behaviour
of these systems. In chapter 5, some of these methods were applied to a disordered array of point dipoles, in order to clarify the effect of the interaction’s anisotropy in the
behaviour of the magnetic correlations between particles.
In Chapter 6, the point dipole approximation is dropped and finite size corrections for
different particle shapes are considered. These results are used to explain the anisotropy
found in ferromagnetic resonance experiments on a square lattice of magnetic disks, when
the orientation of the applied field is varied with relation to the lattice axis. It is shown
7
Abstract
that the inclusion of these finite size corrections, mantaining an uniform magnetization
and considering interactions between disks, leads to a distribution of internal fields that
already exhibits the anisotropy that had been found in the experiments. A variational
procedure is presented to estimate the deformation of the magnetization inside each disk,
due not only to the shape but also to the interaction with other disks.
Finally, in Chapter 7, we present our conclusions and open questions related to this
work.
8
Conteúdo
Agradecimentos
3
Resumo
5
Abstract
7
1 Introdução
13
1.1
Aplicações . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
1.2
Modelos teóricos e computacionais
1.3
Notação . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
1.4
Unidades . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
1.5
Conteúdo original e material externo . . . . . . . . . . . . . . . . . . . . . 22
2 Métodos analíticos
. . . . . . . . . . . . . . . . . . . . . . 18
23
2.1
Propriedades termodinâmicas . . . . . . . . . . . . . . . . . . . . . . . . . 24
2.2
Momentos livres . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
2.3
Interacções como perturbação . . . . . . . . . . . . . . . . . . . . . . . . . 26
2.4
Expansão de altas temperaturas . . . . . . . . . . . . . . . . . . . . . . . . 27
2.5
Conclusões . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
3 Métodos numéricos
41
3.1
Revisão de métodos de Monte Carlo . . . . . . . . . . . . . . . . . . . . . . 42
3.2
Observáveis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
3.3
Ferramentas de computação científica . . . . . . . . . . . . . . . . . . . . . 55
3.4
Conclusões . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
4 Redes regulares de dipolos pontuais
65
4.1
Rede quadrada . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
4.2
Cadeia unidimensional . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
4.3
Conclusões . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
9
Conteúdo
5 Plano desordenado de dipolos pontuais
5.1 Modelo . . . . . . . . . . . . . . . . .
5.2 Efeitos de tamanho finito . . . . . . .
5.3 Simulações . . . . . . . . . . . . . . .
5.4 Magnetização . . . . . . . . . . . . .
5.5 Correlações magnéticas . . . . . . . .
5.6 Conclusões . . . . . . . . . . . . . . .
6 Partículas de tamanho finito
6.1 Função de forma . . . .
6.2 Magnetização uniforme .
6.3 Rede de discos . . . . .
6.4 Conclusões . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
79
79
81
83
85
87
89
.
.
.
.
91
92
95
99
108
7 Considerações finais
111
A Médias livres
115
B Amplitudes de forma
121
C Lista de publicações
125
Bibliografia
131
10
Lista de Figuras
1.1
Métodos de fabrico de sistemas experimentais . . . . . . . . . . . . . . . . 15
1.2
Esquema STM
1.3
Esquema MFM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
1.4
Esquema FMR . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
1.5
Aproximação de dipolos pontuais . . . . . . . . . . . . . . . . . . . . . . . 19
2.1
Representação diagramática de partições de conjuntos . . . . . . . . . . . . 37
2.2
Representação diagramática da correcção de segunda ordem à energia livre. 38
2.3
Representação diagramática de derivadas de partições . . . . . . . . . . . . 38
2.4
Representação diagramática da correcção de segunda ordem à magnetização. 38
2.5
Representação diagramática da correcção de segunda ordem à correlação. . 38
3.1
Peso estatístico em função da energia . . . . . . . . . . . . . . . . . . . . . 42
3.2
Representação esquemática da selecção uniforme . . . . . . . . . . . . . . . 44
3.3
Representação esquemática da selecção com campo externo . . . . . . . . . 46
3.4
Representação esquemática da selecção com componente planar . . . . . . 47
3.5
Histogramas em função da energia . . . . . . . . . . . . . . . . . . . . . . . 49
3.6
Histograma multicanónico em função da energia . . . . . . . . . . . . . . . 50
3.7
Escolha de temperaturas para Parallel Tempering . . . . . . . . . . . . . . 50
3.8
Ferramentas de computação científica. . . . . . . . . . . . . . . . . . . . . 56
4.1
Estado de micro-vórtice e os casos particulares de Luttinger e Tisza. . . . . 67
4.2
Representação esquemática das condições de fronteira periódicas. . . . . . 68
4.3
Parâmetro de ordem e calor específico na rede quadrada
4.4
Susceptibilidade e cumulante de Binder na rede quadrada . . . . . . . . . . 70
4.5
Exemplo de cálculo de tempos de autocorrelação . . . . . . . . . . . . . . . 71
4.6
Tempos de autocorrelação com Metropolis convencional . . . . . . . . . . . 71
4.7
Comparação entre dinâmicas na rede quadrada . . . . . . . . . . . . . . . . 72
4.8
Comparação entre métodos de selecção na rede quadrada . . . . . . . . . . 72
4.9
Comparação de tempos de autocorrelação com diferentes dinâmicas . . . . 73
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
. . . . . . . . . . 69
11
Lista de Figuras
4.10
4.11
4.12
4.13
Cumulante de Binder na rede quadrada . . . . . . . . . . . . . . . . . . . .
Colapsos das curvas de magnetização, cumulante de Binder e susceptibilidade.
Magnetização e cumulante de Binder para uma cadeia unidimensional . . .
Tempos de autocorrelação para uma cadeia unidimensional . . . . . . . . .
74
75
77
77
5.1
5.2
5.3
5.4
5.5
5.6
5.7
5.8
5.9
5.10
5.11
5.12
5.13
Representação esquemática de cada amostra desordenada. . . .
Distribuição de partículas num plano aleatório . . . . . . . . . .
Distribuição da magnetização com condições de fronteira abertas
Representação esquemática da “esfera de Lorentz” . . . . . . . .
Convergência do campo de Lorentz . . . . . . . . . . . . . . . .
Distribuição da magnetização com condições de fronteira abertas
Taxas de aceitação no sistema desordenado . . . . . . . . . . . .
Magnetização em função do campo aplicado . . . . . . . . . . .
Magnetização em função do campo aplicado . . . . . . . . . . .
Magnetização em função da temperatura e concentração . . . .
Exemplo da distribuição da correlação magnética . . . . . . . .
Correlação em função do campo aplicado . . . . . . . . . . . . .
Correlação das componentes paralelas ao campo . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
80
80
81
82
82
83
84
85
86
86
87
88
88
6.1
6.2
6.3
6.4
6.5
6.6
6.7
6.8
6.9
6.10
Representação esquemática do teorema da reciprocidade. . . . . . . . . .
Representação esquemática da função de forma. . . . . . . . . . . . . . .
Representação esquemática de um corpo com magnetização não uniforme.
Representação esquemática de um corpo com magnetização uniforme. . .
Representação esquemática da magnetização distribuída . . . . . . . . .
Representação esquemática de uma rede quadrada de discos . . . . . . .
Campo de ressonância em função do ângulo do campo aplicado . . . . . .
Campo de ressonância e anisotropia em função do parâmetro de rede . .
Intensidade de absorção em função do campo . . . . . . . . . . . . . . . .
Deformação da magnetização num disco . . . . . . . . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
91
92
93
95
97
99
100
101
105
108
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
A.1 Ângulo entre o momento e o campo externo . . . . . . . . . . . . . . . . . 115
B.1 Esquema de esfera com magnetização uniforme . . . . . . . . . . . . . . . . 121
B.2 Esquema das configurações de um cilindro . . . . . . . . . . . . . . . . . . 122
12
1 Introdução
À escala atómica, a interacção magnetostática entre os momentos magnéticos atómicos
é fraca quando comparada com a interacção de troca, principal responsável pela ordem
magnética dos materiais. Ainda assim, o carácter de longo alcance da interacção magnetostática pode alterar as propriedades magnéticas de materiais com a ordem microscópica
já assegurada pela interacção de troca, como a formação de domínios, anisotropias de
forma e modos de inversão magnética.
Já entre os sistemas constituídos por partículas de domínio único (normalmente de dimensões nanométricas), com um número significativo, mas contável, de momentos atómicos ordenados ferromagneticamente, a interacção magnetostática assume uma relevância
adicional, tanto dentro das partículas, onde pode ter uma magnitude comparável à da
interacção de troca, como entre as partículas, podendo inclusivamente ser responsável
pela ordem magnética do sistema composto à temperatura ambiente.
Estas nano-estruturas magnéticas constituídas por partículas separadas possuem propriedades estáticas e dinâmicas muito distintas das propriedades de amostras macroscópicas dos materiais que as constituem. A dimensão e a forma de pequenas partículas
magnéticas afectam não só a distribuição e dinâmica internas da magnetização (e.g. Cowburn [1]) como também a interacção entre as partículas (e.g. Guslienko [2]).
Compreender estes efeitos, para além do interesse do ponto de vista fundamental, permite a concepção de novos materiais cujo comportamento pode ser ajustado à medida de
cada aplicação particular, e é esta a principal motivação de trabalhos nesta área como os
que são apresentados nesta dissertação.
Algumas destas aplicações são apresentadas na secção que se segue. Posteriormente,
serão apresentados os modelos teóricos e computacionais utilizados para estudar estes
sistemas.
1.1 Aplicações
As aplicações tecnológicas de materiais magnéticos nano-estruturados encontram-se em
áreas tão distintas como o fabrico de ímans mais fortes e a criação de sensores biomédicos,
13
1 Introdução
passando por meios de armazenamento digital de maior densidade, memórias magnéticas,
dispositivos de lógica programável, sistemas micro-electromecânicos, etc. (e.g. Moriarty [3], Bader [4, 5], Skomski [6]).
Por exemplo, no domínio do armazenamento de dados digitais em meios magnéticos, o
desenvolvimento de sensores com maior resolução graças aos efeitos de magneto-resistência
gigante (e.g. Prinz [7]) permitiu um aumento exponencial nas densidades de gravação dos
discos que usamos nos nossos computadores.
No entanto, esse aumento está limitado por instabilidades térmicas no meio de armazenamento, e tem sido referida a possibilidade de usar sistemas regulares de nanopartículas
de domínio único (patterned media) como uma forma de ultrapassar este limite (e.g.
Ross [8], Moser et al. [9], Adeyeye e Singh [10], Chien et al. [11]).
As memórias magnéticas de acesso aleatório (MRAM, e.g. Tehranie et al. [12, 13]) têm
como principal característica diferenciadora das memórias actuais, baseadas em semicondutores, a não-volatilidade, sem prejudicar, ou até melhorando, as velocidades de acesso
e escritas convencionais. Isto permitiria, por exemplo, o desenvolvimento de computadores que podem ser ligados e desligados restaurando instantaneamente o estado anterior,
sem o procedimento de gravação e leitura no disco que actualmente permite um resultado
semelhante.
O conceito de dispositivos lógicos com nano-estruturas magnéticas (e.g. Cowburn e
Welland [14], Csaba et al. [15], Allwood et al. [16], Imre et al. [17]) foi proposto como
alternativa aos dispositivos actuais, que processam a informação de forma electrónica.
Usando redes de partículas magnéticas para registar estados lógicos e excitações magnéticas para propagar informação através da rede, alimentada por um campo magnético
oscilatório que serve também como relógio, estas propostas prometem aumentos ambiciosos na densidade de integração e reduções na dissipação de energia quando comparados
com os sistemas microelectrónicos actuais.
No entanto, estas eventuais aplicações requerem o conhecimento detalhado dos efeitos
colectivos das interacções magnéticas na sua forma mais relevante, a interacção dipolar,
com as particularidades de serem anisotrópicas e de longo alcance, o que não acontece
com a interacção de troca.
Para além da compreensão teórica do comportamento destas estruturas, tem sido necessário desenvolver métodos de preparação e observação destes sistemas regulares em grande
escala, bem como sistemas de leitura, gravação e processamento diferentes dos actuais, o
que constitui por si só uma área de investigação activa. Estes métodos de preparação e
observação são apresentados nas sub-secções seguintes.
14
1.1 Aplicações
Sistemas experimentais
Os métodos de fabrico são normalmente [18] divididos entre os métodos top-down (e.g.,
litografia, deposição) e bottom-up (e.g., self-assembly), e tanto podem produzir sistemas
ordenados, por exemplo para sistemas de armazenamento magnético (e.g. Ross [8]), como
desordenados, cujas aplicações se encontram mais frequentemente em sensores. Estes
sensores exploram a ligação entre as propriedades magnéticas e eléctricas destes materiais,
nomeadamente a magneto-resistência gigante (GMR, Berkowitz et al. [19]) semelhante à
inicialmente observada em multicamadas magnéticas (Baibich et al. [20]).
Top-down
Bottom-up
(e.g. litografia)
(e.g. self-assembly)
Figura 1.1: Ilustração de métodos de fabrico de sistemas experimentais.
Métodos top-down Os métodos top-down são tipicamente evoluções, em termos de resolução, dos processos convencionais de microfabricação utilizados pela indústria
da microelectrónica e processadores. Não são (ainda) os métodos imaginados por
Feynman [21], na famosa palestra There’s plenty of room at the bottom, que sugeriu
uma recursão de máquinas que construiriam sucessivamente máquinas mais pequenas até à escala atómica, mas, ainda assim, já demonstraram um controlo e precisão
considerável. O principal desafio nesta área é o fabrico de estruturas regulares em
grande escala.
Métodos bottom-up Os métodos bottom-up utilizam propriedades físicas e químicas de
átomos e moléculas, nomeadamente as taxas anisotrópicas de difusão na superfície
atomicamente lisa do substrato e as interacções entre os átomos no substrato, para
fazer com que estes se auto-organizem para formar uma determinada configuração
estrutural e/ou posicional (Drexler [22]). A ideia é emular os sistemas biológicos,
onde “dispositivos” baseados em proteínas se auto-organizam para formar “máquinas” moleculares. A capacidade de replicação que encontramos nos sistemas
biológicos poderia permitir a produção em paralelo de estruturas macroscópicas.
15
1 Introdução
Aqui, a principal dificuldade é manter o controlo preciso das estruturas moleculares
à medida que as dimensões e complexidade aumentam.
Métodos híbridos Existem também métodos de preparação híbridos, que tipicamente
combinam uma primeira fase top-down para definir uma máscara onde elementos
atómicos e moleculares são depositados, seguidos de uma fase de auto-organização
confinada à estrutura preparada inicialmente.
Imagiologia
A engenharia de sistemas à escala nanométrica exige novas formas de observar estes sistemas, novos nanoscópios que vão além do uso comum da luz visível, dado que esta não
permite resoluções muito inferiores a um micrómetro. Estas novas formas de observação
exploram efeitos eléctricos e magnéticos à escala atómica para produzir imagens com a resolução necessária. À invenção do microscópio de efeito de túnel (STM, Binnig et al. [23],
figura 1.2), que permite não só a observação mas também, em certas condições, a manipulação à escala atómica, seguiu-se o desenvolvimento de tecnologias relacionadas como
a microscopia de força atómica (AFM, Binnig et al. [24]) e a microscopia de força magnética (MFM, Martin e Wickramasinghe [25], figura 1.3), com as suas várias modificações,
tecnologias estas que marcaram o nascimento da nanociência e nanotecnologia.
VT
IT
Figura 1.2: Representação esquemática da experiência de STM.
Figura 1.3: Representação esquemática da experiência de MFM.
16
1.1 Aplicações
Ressonância Ferromagnética Existem muitos outros mecanismos de caracterização das
propriedades estruturais, eléctricas e magnéticas de sistemas experimentais. Pela
sua relevância na parte final desta dissertação, as medições de ressonância ferromagnética (FMR, Griffiths [26], Kittel [27]), desenvolvidas muito antes do surgimento
das nanotecnologias, merecem aqui um destaque especial. A experiência de ressonância ferromagnética consiste em colocar a amostra na parede de uma cavidade
rectangular que termina um guia de ondas alimentado por um gerador de microondas, de tal forma que o campo magnético das micro-ondas esteja sempre no plano
da amostra (figura 1.4). É aplicado um campo magnético constante, podendo situarse tanto no plano da amostra como sob um ângulo φ com este, mas perpendicular ao
campo magnético das micro-ondas, campo este que alinha os momentos magnéticos
da amostra induzindo uma precessão a uma frequência que depende da intensidade
do campo magnético aplicado e do campo desmagnetizante da amostra, relacionado
com as interacções magnetostáticas dentro desta. À medida que se varia o campo
aplicado, verifica-se um aumento abrupto na perda de energia com a passagem da
onda através da cavidade, correspondente à condição de ressonância entre a frequência de micro-ondas do campo alternado e a precessão dos momentos magnéticos, o
que nos dá informação sobre os campos locais no interior da amostra.
Figura 1.4: Representação esquemática da experiência de FMR.
17
1 Introdução
1.2 Modelos teóricos e computacionais
Para estudar o efeito das interacções entre partículas magnéticas de domínio único, objectivo que constitui a maior componente desta dissertação, o modelo mais simples passa por
substituir as partículas pelo seu momento magnético, colocando um momento magnético
pontual, de magnitude igual à soma dos seus momentos atómicos, no centro de cada partícula. Em tudo o que se segue, assumimos que a magnitude deste momento magnético
total é suficiente para justificar uma abordagem puramente clássica ao problema, dado
que os sistemas experimentais que pretendemos modelar são constituídos por partículas
com, pelo menos, milhares de átomos ordenados ferromagneticamente.
Momentos magnéticos
Como não existem cargas magnéticas livres, a entidade mais simples que podemos usar
para descrever materiais magnéticos é o dipolo magnético pontual, ou seja, considerar
a existência dos dois pólos magnéticos no mesmo ponto do espaço. Excepto no caso
das partículas elementares, isto é uma aproximação. Como iremos ver, o tamanho finito
das partículas, mantendo a magnetização uniforme, pode ser tratado considerando uma
separação entre os pólos magnéticos, dando origem a correcções de ordens superiores
à interacção dipolar, ou tendo em conta toda a forma da amostra. Isto dá origem a
correcções à energia de interacção entre dois momentos magnéticos. Já uma distribuição
não uniforme da magnetização dentro das partículas pode obrigar-nos a abandonar por
completo este modelo de interacção entre partículas magnéticas.
Dipolos pontuais
No caso de considerarmos as partículas como dipolos pontuais, e se designarmos por µi
o momento magnético da partícula i e por r ij a posição relativa entre os centros das
partículas i e j, a energia de interacção entre duas partículas é dada por
E=
µi · µj
−3
3
rij
(µi · r ij ) µj · r ij
5
rij
.
(1.1)
Esta aproximação é exacta no caso de partículas esféricas uniformemente magnetizadas
(e.g. Jackson [28]), mas só é razoável assumir esta uniformidade em partículas muito
pequenas e/ou campos muito fortes. Para partículas não esféricas, mesmo que uniformemente magnetizadas, existem correcções à interacção entre as partículas que podem alterar
as simetrias e o decaimento da interacção, particularmente a distâncias comparáveis com
o tamanho da partícula.
18
1.2 Modelos teóricos e computacionais
Figura 1.5: Dipolos pontuais como modelo da interacção entre partículas.
A dimensionalidade do momento magnético também pode ser relevante (e.g. [29]),
considerando-se momentos com apenas duas orientações possíveis (correspondente ao modelo de Ising), livres para rodar no plano (modelo XY) ou no espaço (modelo Heisenberg).
Em sistemas planares, dado que o estado fundamental será também planar, os casos XY e
Heisenberg são muitas vezes considerados equivalentes, mas o grau de liberdade adicional
fora do plano afecta as excitações a temperatura finita e pode ser relevante, pelo que
iremos considerar este último caso.
Expansão de altas temperaturas
Dado que a interacção entre as partículas pode assumir diferentes formas, no capítulo 2
desenvolvemos uma expansão de altas temperaturas que assume apenas uma interacção
da forma
− Uijαβ µαi µβj
(1.2)
ou seja, que é possível escrever um tensor de interacção que contém toda a informação
excepto a direcção e magnitude dos momentos magnéticos das partículas. O sinal negativo
é apenas uma convenção que nos permite, se definirmos um campo
Hjβ =
X
Uijαβ µαi
(1.3)
i,α
escrever a energia em termos do campo local na forma canónica
E=−
X
Hjβ µβj .
(1.4)
j,β
No caso da interacção U corresponder à interacção dipolar, este tensor será dado explicitamente por
Uijαβ
3uαij uβij − δαβ
,
= δ ij
3
rij
(1.5)
19
1 Introdução
com δ ij = (1 − δij ) e uij = r ij /rij , r ij = r i − r j .
Como já foi mencionado, as duas características mais importantes desta matriz de
interacção são o facto de ser anisotrópica e de longo alcance. Ao contrário de outros
Hamiltonianos magnéticos, a anisotropia não vem simplesmente do facto de a interacção
ter magnitudes diferentes para componentes diferentes - a interacção depende da posição
relativa das partículas, de tal forma que tanto pode ser ferromagnética como antiferromagnética.
Simulações de Monte Carlo
O objectivo do nosso desenvolvimento da expansão de altas temperaturas foi validar, no
limite correspondente a este regime, o resultado de simulações de Monte Carlo. Desde
que mantenhamos a magnetização uniforme e de igual magnitude em todas as partículas,
tanto a expansão de altas temperaturas como as simulações de Monte Carlo são genéricas,
bastando alterar o tensor de interacção para estudar outros sistemas, e com ambas as
técnicas tratamos as interacções de longo-alcance sem qualquer cut-off. Este aspecto
contrasta com diversos trabalhos, tanto analíticos como numéricos, que ou limitam a
interacção a próximos vizinhos ou têm em conta o decaimento da interacção para, no
limite termodinâmico, excluir correcções de ordem superior. Se por um lado o facto de
as simulações de Monte Carlo tratarem sempre sistemas finitos já justifica o tratamento
completo das interacções, o facto de estarmos a estudar sistemas nano-estruturados, com
um número contável de partículas, põe em causa a utilização do limite termodinâmico
mesmo em estudos analíticos, pelo que se torna ainda mais relevante esta abordagem.
No capítulo 3, é feita uma revisão sobre métodos de Monte Carlo em geral e sobre
as técnicas adoptadas nesta dissertação em particular. São discutidos os problemas de
convergência e baixas taxas de aceitação das simulações, e as formas utilizadas para os
ultrapassar. Neste capítulo, será ainda apresentado um levantamento de ferramentas de
computação científica utilizadas.
Redes regulares de dipolos pontuais
No caso de a interacção dipolar ser a única interacção entre partículas colocadas numa rede
regular, a combinação das simetrias da rede e da interacção produz comportamentos não
triviais, reconhecidos pela primeira vez por Luttinger e Tisza [30], começando pelo facto de
o estado ferromagnético não ser o estado fundamental para uma simples rede quadrada,
e que continuam a surpreender. Para estudar as consequências termodinâmicas destes
efeitos, os métodos de Monte Carlo descritos anteriormente serão aplicados, no capítulo
20
1.3 Notação
4, ao estudo de uma rede quadrada de dipolos pontuais, onde poderemos também verificar
a importância da escolha de diferentes dinâmicas de simulação.
Sistemas desordenados
Como vimos, alguns dos sistemas experimentais com relevância tecnológica são constituídos por distribuições desordenadas de grânulos magnéticos, cuja presença define as
propriedades eléctricas e magnéticas do sistema composto. Nomeadamente, a interacção
magnetostática entre os grânulos introduz correlações magnéticas que se manifestam na
resistência eléctrica dos materiais. No capítulo 5, estudamos as correlações magnéticas
em planos com distribuições aleatórias de dipolos pontuais, com vista a caracterizar estes
efeitos.
Tamanho finito das partículas
Em geral, a dimensão e a forma de uma partícula produz um campo desmagnetizante que
provoca a deformação da magnetização no interior desta. Para além disso, as interacções
magnetostáticas entre partículas podem dar origem a importantes correcções a este campo
desmagnetizante, e é este o tópico da última parte da dissertação. Será apresentado
um mecanismo genérico de introdução da forma das partículas no cálculo da energia
magnetostática. Este mecanismo é utilizado para calcular a distribuição de campos e
a deformação da magnetização no interior de partículas com, neste caso em particular,
a forma de discos e dispostas numa rede quadrada, com o objectivo de explicar uma
anisotropia de quarta ordem verificada em experiências de ressonância ferromagnética.
1.3 Notação
Ao longo da dissertação, adoptamos a seguinte notação:
• o momento magnético de uma partícula é designado por µ
• os índices que identificam partículas são letras do alfabeto latino, colocadas em
subscripto, e.g. µi
• os índices que identificam componentes geométricas são letras do alfabeto grego,
colocadas em sobrescrito, e.g. rα
• vectores são representados a negrito, e.g. r = (rx , ry , rz )
21
1 Introdução
• índices de componentes geométricas repetidos no mesmo termo definem implicitaP
P
mente uma soma nesses índices, e.g. i Uijαβ µαi ≡ i,α Uijαβ µαi
• o valor esperado de uma variável aleatória A é designado por hAi
• o momento estatístico de ordem n de uma variável aleatória A é designado por hAn i
• o cumulante estatístico de ordem n de uma variável aleatória A é designado por [An ]
Outras notações mais específicas serão introduzidas ao longo do texto, quando apropriado.
1.4 Unidades
Na apresentação de resultados de simulações, nos capítulos 4 e 5, são usadas as unidades
naturais dos modelos estudados. Esta escolha é feita, por um lado, porque os resultados só
dependem de relações entre os parâmetros do modelo, como entre o momento magnético
µ de cada partícula e uma distância característica d do modelo. Por exemplo, no caso de
uma rede quadrada, de lado a, de dipolos pontuais com momento magnético µ, a energia
é apresentada em unidades de µ2 /a3 . Por outro lado, é preferível usar nas simulações
grandezas da ordem da unidade, para evitar perdas de precisão numérica, pelo que consideramos, na temperatura ou no calor específico, que kB , a constante de Boltzmann, é
igual à unidade. Os tempos de simulação são sempre medidos em passos de Monte Carlo.
Já no capítulo 6, onde o objectivo é modelar um sistema experimental em particular, são
usadas unidades do Sistema Internacional.
1.5 Conteúdo original e material externo
Excepto quando especificamente indicado, todo o conteúdo é de autoria própria, incluindo
textos, ilustrações e programas de computador (sem prejuízo da utilização das ferramentas
e bibliotecas indicadas).
22
2 Métodos analíticos
Os sistemas de partículas magnéticas de que nos vamos ocupar podem ser descritos por um
modelo clássico de momentos magnéticos pontuais, colocados no centro de cada partícula,
acoplados por uma interacção de longo alcance e eventualmente sujeitos a um campo
magnético externo. O Hamiltoniano deste modelo pode ser escrito como
H = HI + HZ ,
(2.1)
onde HI representa a interacção entre os momentos magnéticos e o termo de Zeeman HZ
representa a interacção com o campo externo.
O termo de interacção HI pode ser apresentado, sem perda de generalidade, na forma
de uma soma sobre todos os pares de momentos magnéticos
HI = −
1 X χλ χ λ
U µ µ ,
2 k,l kl k l
(2.2)
onde, como foi indicado na secção 1.3, estão implícitos os somatórios sobre índices de
componentes geométricas χ e λ1 . Escolhemos fazer a soma sobre toda a matriz (e dividir
por 2 para não contar duplamente os pares) em vez de somar sobre pares para que, no
que se segue, os índices k e l sejam equivalentes.
De momento, vamos assumir apenas que não existem termos de auto-interacção e que a
interacção de k com l é igual à interacção de l com k (ou seja, que Ukl é simétrica e que os
elementos diagonais são nulos). Qualquer interacção binária satisfaz estas condições, seja
de curto alcance ou longo alcance, isotrópica ou anisotrópica. Os momentos magnéticos,
por simplicidade com igual magnitude, podem estar confinados a um eixo, a um plano ou
livres para rodar no espaço, pelo que engloba qualquer modelo (clássico) de Ising, XY ou
Heisenberg.
Estando nós interessados em sistemas de partículas acopladas por interacção magne1
Vamos reservar os índices i, j e α, β para designar partículas específicas para as quais queremos calcular
médias termodinâmicas. Os índices k, l, m, n... e χ, λ, µ, ν... serão índices mudos dos somatórios
presentes nos Hamiltonianos.
23
2 Métodos analíticos
tostática, a matriz de interacção irá depois ser dada por
Uklχλ = δ kl
3uχkl uλkl − δχλ
,
3
rkl
(2.3)
com δ kl = (1 − δkl ) e ukl = r kl /rkl , r kl = r k − r l , mas no que se segue vamos considerar
uma matriz arbitrária.
Quanto ao termo de Zeeman, HZ , vamos escrevê-lo também de uma forma generalizada,
como
HZ = −
X
µk · hk = −
X χ χ
k
µk hk .
(2.4)
k
onde está implícito o somatório sobre o índice de componente cartesiana χ. Não assumimos nada em relação ao campo local, não só para manter a generalidade mas, principalmente, porque ter um campo local por partícula vai permitir-nos calcular funções
de correlação entre momentos magnéticos como derivadas da energia livre em ordem aos
campos locais. Após este cálculo, os campos locais podem ser igualados a zero para obter as grandezas a campo nulo ou identificados com um campo externo para estudar o
comportamento em função do campo aplicado.
2.1 Propriedades termodinâmicas
As propriedades em que estamos interessados são médias termodinâmicas em sistemas com
um número fixo de partículas em equilíbrio térmico com o ambiente (ensemble canónico,
e.g. Lage [31]). Neste caso, a termodinâmica do sistema é completamente descrita pela
sua função de partição
Z=
onde
dµ e−βH ,
(2.5)
dµ representa o integral sobre todas as configurações possíveis do sistema e
β=
1
,
kB T
(2.6)
ou pela energia livre
1
F = − ln Z.
β
(2.7)
A média termodinâmica de qualquer observável O, tipicamente uma função dos momentos magnéticos µi , é dada pela média dessa grandeza sobre todas as configurações
24
2.2 Momentos livres
possíveis, ponderada pelo peso de Boltzmann de cada configuração,
dµ O e−βH
.
hOi =
−βH
dµ e
(2.8)
Note-se que, se o Hamiltoniano H incluir o termo HZ como definido na equação 2.4,
então as grandezas que dependem dos momentos magnéticos µi podem ser escritas à custa
de derivadas da energia livre em ordem aos campos locais hi . Por exemplo,
hµαi i
∂
F=
=
∂hαi
dµ µαi e−βH
.
dµ e−βH
(2.9)
Na presença de interacções, ainda para mais de longo alcance, estas grandezas tornam-se
difíceis ou mesmo impossíveis de calcular, à medida que aumenta o número de partículas.
No entanto, o sistema sem interacções pode ser resolvido exactamente, como veremos a
seguir.
2.2 Momentos livres
No caso de momentos magnéticos sem interacção (referidos no que se segue como momentos livres), o Hamiltoniano é constituído apenas pelo termo de Zeeman, que é um
somatório de termos que envolvem uma única partícula, pelo que a função de partição
Z0 =
(2.10)
dµ e−βHZ
factoriza, podendo escrever-se como um produtório sobre todas as partículas
Z0 =
N
Y
χ
χ
dµk eβµk Hk .
(2.11)
k=1
O índice “0” aqui introduzido será usado para designar todas as grandezas calculadas com
um Hamiltoniano sem interacções.
De forma análoga à eq. 2.7, a energia livre deste sistema sem interacções apresenta-se
como
1
F0 = − ln Z0
β
(2.12)
e para a análise que se segue, vai ser particularmente relevante definir a média termodinâmica livre,
dµ O e−βHZ
hOi0 = ,
−βH
dµ e
Z
(2.13)
25
2 Métodos analíticos
onde O é uma função arbitrária dos momentos magnéticos µi . De forma análoga à equação 2.9, podemos também apresentar estas médias livres como
hµαi i0
∂
=
F0 =
∂hαi
dµ µαi e−βHZ
.
dµ e−βHZ
(2.14)
Note-se que os termos de HZ que envolverem momentos magnéticos que não estejam
presentes em O vão cancelar no numerador e denominador, pelo que
α α
dµi Oi eβhi µi
.
hOi i0 = α α
dµi eβhi µi
(2.15)
A propriedade mais importante desta média livre é o facto de que uma média que
envolva partículas diferentes pode ser separada num produto de médias que envolvem
apenas uma partícula, ou seja, se i 6= j,
α α
hOi Oj i0 =
=
α α
dµi dµj Oi Oj eβhi µi eβhj µj
α α
α α
dµi dµj eβhi µi eβhj µj
α
βhα
i µi
dµi Oi e
α α
dµi eβhi µi
(2.16)
α
βhα
j µj
dµj Oj e
α α
dµj eβhj µj
ou
hOi Oj i0 = hOi i0 hOj i0 .
(2.17)
Este facto irá ser usado para simplificar as expressões que ocorrem nos termos de ordem
superior de uma expansão de altas temperaturas.
2.3 Interacções como perturbação
Se, na função de partição, separarmos os dois termos do Hamiltoniano da equação 2.1
Z=
−βH
dµ e
=
dµ e−βHI e−βHZ
(2.18)
verificamos que podemos tratar o primeiro termo exponencial como uma “observável” O
da eq. 2.13 no sistema sem interacções, e escrever
D
Z = Z0 e−βHI
e, consequentemente,
26
E
0
(2.19)
2.4 Expansão de altas temperaturas
E
1
1 D
1
F = − ln Z = − ln Z0 − ln e−βHI .
0
β
β
β
(2.20)
A análise que se segue centra-se no último termo desta equação.
2.4 Expansão de altas temperaturas
Quando a temperatura é alta por comparação com a interacção, ou seja, quando βHI 1,
podemos expandir a exponencial na equação 2.20,
D
ln e
−βHI
1
' ln 1 − βHI + (βHI )2 + ...
2
0
!
2 D
E
β
= ln 1 − β hHI i0 +
HI2 + ...
0
2
E
0
(2.21)
e o logaritmo, no mesmo limite,
β 2 D 2E
1
β 2 D 2E
' − β hHI i0 −
HI
−
β hHI i0 −
HI
0
0
2
2
2
h
i
β 2 D 2 E
−β hHI i0 +
HI − hHI i20 + O (βHI )3
0
2
!
!2
+ ...
(2.22)
onde mantivemos apenas os termos até segunda ordem em βHI .
D
A expressão e−βHI
E
0
é a função geradora dos momentos estatísticos (e.g. Kubo [32]),
D
−βHI
e
E
0
(−β)n
hHIn i0
=
n!
n=0
∞
X
(2.23)
e o seu logaritmo é a função geradora dos cumulantes estatísticos (Fisher e Wishart [33]),
D
ln e−βHI
E
0
=
(−β)n n
[HI ]0 ,
n!
n=1
∞
X
(2.24)
o que nos permite escrever a energia livre do sistema com interacções como uma expansão
em cumulantes livres,
∞
1
1X
(−β)n n
[HI ]0 .
F = − ln Z0 −
β
β n=1 n!
(2.25)
Os primeiros termos desta expansão serão
27
2 Métodos analíticos
β 2 h 3i
1
β h 2i
HI +
HI + ...
F ≈ − ln Z0 + [HI ]0 −
0
0
β
2
6
(2.26)
ou, alternativamente, em termos dos momentos estatísticos,
1
β D 2 E
F ≈ − ln Z0 + hHI i0 −
HI − hHI i20 + ... .
0
β
2
(2.27)
As grandezas que vamos estar interessados em calcular podem escrever-se à custa de
derivadas da energia livre em ordem aos campos locais hk . Por exemplo, a componente
α da magnetização da partícula i é dada por
hµαi i = −
∂
∂
∂
F = − α F0 − α F1 − ....,
α
∂hi
∂hi
∂hi
(2.28)
onde Fn = (−β)n−1 [HIn ]0 /n!, para n > 0, e a correlação entre os momentos magnéticos
das partículas i e j escreve-se como
D
E
D
E
µαi µβj − hµαi i µβj = −
1 ∂ ∂
1 ∂ ∂
1 ∂ ∂
F1 − ...
βF = −
β F0 −
α
α
β ∂hi ∂hj
β ∂hi ∂hj
β ∂hαi ∂hβj
(2.29)
No que se segue, vamos calcular explicitamente as correcções de primeira e segunda ordem
à magnetização e à função de correlação, designando por Miα (n) a correcção de ordem n
ao valor médio da componente α da magnetização da partícula i
Miα (n) = −
∂
Fn
∂hαi
(2.30)
e por Cijαβ (n) a correcção de ordem n à correlação média entre as componentes α e β dos
momentos magnéticos das partículas i e j,
Cijαβ (n) = −
1 ∂ ∂
Fn .
β ∂hαi ∂hβj
(2.31)
Naturalmente,
hµαi i
=
∞
X
Miα (n)
(2.32)
n=0
e
D
E
D
E
µαi µβj − hµαi i µβj =
∞
X
n=0
de acordo com as expressões 2.28 e 2.29.
28
Cijαβ (n)
(2.33)
2.4 Expansão de altas temperaturas
Notação simplificada
Para simplificar a notação, vamos introduzir definições específicas para os momentos estatísticos relevantes:
hµχk i0 ≡ mχk
D
µχk µλk
D
E
µχk µλk µµk
0
E
0
≡ mχλ
k
≡ mχλµ
k
(2.34)
(2.35)
(2.36)
e assim sucessivamente. Estas grandezas dependem dos graus de liberdade dos momentos
magnéticos (Apêndice A), mas são sempre função dos campos locais hαi e da temperatura.
Quanto às derivadas em ordem aos campos locais, é importante verificar que a média livre
de qualquer observável de uma partícula só depende do campo local nessa partícula, pelo
que as derivadas de médias de observáveis da partícula k em ordem a observáveis da
partícula i só não são nulas se k = i. Sendo assim, podemos introduzir as notações para
as derivadas correspondentes,
∂
hµχk i0 ≡ δik ∂α mχi
∂hαi
(2.37)
∂ D χ λE
µk µk ≡ δik ∂α mχλ
i
α
0
∂hi
(2.38)
∂ D χ λ µE
µk µk µk ≡ δik ∂α mχλµ
i
α
0
∂hi
(2.39)
e assim sucessivamente.
2.4.1 Termo de ordem zero
O termo de ordem zero da expansão 2.25 corresponde ao sistema sem interacções, logo a
energia livre é simplesmente dada por
1
F0 = − ln Z0 .
β
(2.40)
A média de qualquer momento magnético é dada pela equação 2.28
Miα (0) = hµαi i0 = −
∂
F = mαi
α 0
∂hi
(2.41)
29
2 Métodos analíticos
e a correlação entre os momentos magnéticos de duas partículas diferentes resulta nula
D
Cijαβ (0) = µαi µβj
E
D
0
− hµαi i0 µβj
E
= −
0
1 ∂ ∂
1 ∂
mα = 0,
β F0 =
α
β ∂hi ∂hj
β ∂hβj i
(2.42)
dado que mαi não depende de hβj .
2.4.2 Correcção de primeira ordem
Segundo a equação 2.25, a correcção de primeira ordem à energia livre é dada pela média
livre do Hamiltoniano de interacção
F1 = hHI i0 = −
1 X χλ D χ λ E
U
µk µl ,
0
2 k,l kl
(2.43)
onde, mais uma vez, está implícita a soma sobre os índices de componentes χ e λ. Como
a matriz de interacção obriga a que k seja diferente de l, obtemos
F1 = −
1 X χλ χ λ
U m m .
2 k,l kl k l
(2.44)
Este resultado pode ser usado para calcular a correcção de primeira ordem à magnetização média da partícula i, que é dada pela derivada em ordem ao campo local,
hµαi i1 = −
∂
1 ∂ X χλ χ λ
F
=
U m m
1
∂hαi
2 ∂hαi k,l kl k l
(2.45)
Derivando cada um dos termos do produto de magnetizações,
Miα (1) = hµαi i1 =
1 X χλ χ 1 X χλ λ
Ukl mk δli ∂α mλl +
U m (δik ∂α mχk )
2 k,l
2 k,l kl l
(2.46)
e notando que os índices k e l são equivalentes, podemos escrever
hµαi i1 =
X
χλ χ
Uki
mk ∂α mλi ,
(2.47)
k
onde agora o somatório sobre partículas é só sobre o índice k.
A correcção de primeira ordem à correlação entre as partículas i e j obtém-se directamente a partir da correcção de primeira ordem à magnetização,
30
2.4 Expansão de altas temperaturas
Cijαβ (1) = −
1 ∂ ∂
1 ∂
α
F
=
1
β
β hµi i1
α
β ∂hi ∂hj
β ∂hj
1 ∂
=
β ∂hβj
(
X
χλ χ
Uki
mk
∂α mλi
(2.48)
)
.
k
A derivada só não é nula quando k = j, logo
Cijαβ (1) =
1 χλ Uji ∂β mχj ∂α mλi ,
β
(2.49)
reduzindo-se portanto a uma forma separável de médias termodinâmicas livres. Esta
forma separável será sempre o objectivo também no tratamento dos termos de ordens
superiores. Note-se que, em primeira ordem, a correlação entre i e j não envolve nenhum
somatório sobre as restantes partículas, apenas a interacção entre i e j. Nos termos de
ordem superior, esperamos encontrar o efeito da presença das restantes partículas.
2.4.3 Correcção de segunda ordem
A correcção de segunda ordem à energia livre envolve médias de quatro momentos magnéticos, na seguinte forma
1 D E
F2 = − β HI2 − hHI i20
0
2
D E
1 X X χλ µν D χ λ µ ν E
µν
= − β
Ukl Umn µk µl µm µn − Uklχλ hµχk i0 µλl Umn
hµµm i0 hµνn i0
0
0
8 k,l m,n
D E
1 X χλ µν D χ λ µ ν E
Ukl Umn µk µl µm µn − hµχk i0 µλl hµµm i0 hµνn i0 .
= − β
0
0
8 k,l,m,n
(2.50)
Nos casos em que todos os índices k, l, m e n são diferentes, os dois termos são iguais
e a expressão é nula, pelo que só temos de considerar os casos restantes. Este procedimento pode ser automatizado particionando o conjunto {k, l, m, n} em blocos de tamanho
b = 1 até 4 (e.g., com a função KSetP artitions do Mathematica), mas, numa primeira
exposição, vamos tratar este termo de segunda ordem manualmente.
Em primeiro lugar, verificamos que nunca ocorrem três nem quatro índices iguais,
porque as matrizes de interacção obrigam a que k 6= l e m 6= n. Portanto, as únicas
composições não nulas são as de três índices diferentes (4 casos):
31
2 Métodos analíticos
k 6= m, l = n
k 6= n, l = m
l 6= m, k = n
l 6= n, k = m
e dois índices diferentes (2 casos, dado que são proibidos os que implicam k=l e/ou
m=n):
k = m, l = n
k = n, l = m
Como em F2 os índices são todos equivalentes, podemos fazer uma rotação cíclica dos
índices e escrever:
D E
D
E
1 X χλ µν h
Ukl Umn 4δlm δ kn hµχk i0 hµνn i0 µλl µµm − µλl hµµm i0 (2.51)
F2 = − β
0
0
8 k,l,m,n
D
+2δlm δkn hµχk µνn i0 µλl µµm
E
D
0
− hµχk i0 µλl
E
0
hµµm i0 hµνn i0
i
.
Contraindo os índices iguais (nas partículas) e, usando a notação da secção 2.4, chegamos à forma
1 X χλ µν
λ µ
Ukl Uln δ kn mχk mνn mλµ
F2 = − β
l − ml ml
2 k,l,n
1 X χλ µν χν λµ
− β
Ukl Ulk mk ml − mχk mλl mµl mνk .
4 k,l
(2.52)
Esta equivalência entre as equações 2.50 e 2.52 foi verificada numericamente para um
sistema de três partículas com interacção dipolar.
Magnetização
A correcção de segunda ordem à magnetização média da partícula i é dada por
∂
F2
∂hαi
1 ∂ X χλ µν
λ µ
−
m
m
=
β α
Ukl Uln δ kn mχk mνn mλµ
l
l
l
2 ∂hi k,l,n
Miα (2) = hµαi i2 = −
+
32
1 ∂ X χλ µν χν λµ
β α
Ukl Ulk mk ml − mχk mλl mµl mνk .
4 ∂hi k,l
(2.53)
2.4 Expansão de altas temperaturas
Separando nos casos em que cada um dos índices na soma é igual a i, temos a série de
contribuições

X χλ µν
1
λ µ
Ukl Uln δ kn (δki ∂α mχk ) mνn mλµ
hµαi i2 = β 
l − ml ml
2 k,l,n
+
X
µν
λ µ
Uklχλ Uln
δ kn mχk (δni ∂α mνn ) mλµ
l − ml ml
k,l,n
+
X
µν
Uklχλ Uln
δ kn mχk mνn
n
δli ∂α mλµ
− δli ∂α mλl mµl
l
o
(2.54)
k,l,n
+
o
1 X χλ µν n
λµ
χ ν
λ µ
Ukl Ulk (δki ∂α mχν
k ) ml − ml ml (δki ∂α mk mk )
2 k,l
+
1
2

X
n
Uklχλ Ulkµν mχν
δli ∂α mλµ
− mχk mνk ∂α mλl mµl
k
l
o
.
k,l
Contraindo os índices iguais e renumerando índices equivalentes, podemos simplificar
para

hµαi i2
1  X χλ µν
λ µ
Ukl Uli δ ki mχk (∂α mνi ) mλµ
=
β 2
l − ml ml
2
k,l
+
X
χλ µν
Uki
Uin δ kn mχk mνn
n
∂α mλµ
− ∂α mλi mµi
i
o
(2.55)
k,n
+
X
χλ µν
Uik
Uki
n
mχν
k
∂α mλµ
i
−
mχk mνk
∂α mλi mµi
o
#
,
k
que pode ser interpretado como, no primeiro termo, as interacções entre todas as partículas k com a partículas i mediadas por todas as partículas l, no segundo termo a
interacção entre todas as partículas k e l mediadas pela partículas i e, por último, uma
auto-interacção de i mediada por todas as partículas k.
Correlação
Mais uma vez, a correcção de segunda ordem à correlação entre as partículas i e j, é dada
pela derivada em ordem ao campo local de j da magnetização de i,
Cijαβ (2) = −
1 ∂ ∂
1 ∂
hµα i .
β F2 =
α
β ∂hi ∂hj
β ∂hβj i 2
(2.56)
Substituindo o termo de segunda ordem da magnetização, já calculada na eq. 2.55,
33
2 Métodos analíticos
obtemos

1 ∂  X χλ µν
χ
λµ
ν
λ µ
(∂
)
2
U
U
δ
m
m
m
−
m
m
Cijαβ (2) =
α i
l
kl
li ki k
l
l
2 ∂hβj
k,l
+
X
χλ µν
Uki
Uin δ kn mχk mνn
n
∂α mλµ
− ∂α mλi mµi
i
o
(2.57)
k,n
+
X
χλ µν
Uki
Uik
n
mχν
k
∂α mλµ
i
−
mχk mνk
∂α mλi mµi
#
o
.
k
Separando nos casos em que cada um dos índices é igual a j,

1  X χλ µν
λ µ
Cijαβ (2) =
2
Ukl Uli δ ki (δkj ∂β mχk ) (∂α mνi ) mλµ
−
m
m
l
l
l
2
k,l
+ 2
X
Uklχλ Uliµν δ ki mχk (∂α mνi )
n
o
o
δlj ∂β mλµ
− δlj ∂β mλl mµl
l
o
k,l
+
X
χλ µν
Uin δ kn (δkj ∂β mχk ) mνn
Uki
n
χλ µν
Uin δ kn mχk (δnj ∂β mνn )
Uki
n
∂α mλµ
− ∂α mλi mµi
i
(2.58)
k,n
+
X
∂α mλµ
− ∂α mλi mµi
i
k,n
+
X
χλ µν
Uki
Uik
n
(δkj ∂β mχν
k )
∂α mλµ
i
−
∂α mλi mµi
#
o
(δkj ∂β mχk mνk )
k
e, contraindo os índices iguais e rodando índices equivalentes, podemos simplificar para
Cijαβ (2) =
X
+
X
λ µ
Ujlχλ Uliµν δ ji ∂β mχj (∂α mνi ) mλµ
l − ml ml
l
χλ µν
Uji δ ki mχk (∂α mνi )
Ukj
n
∂β mλµ
− ∂β mλj mµj
j
o
(2.59)
k
+
X
χλ µν
Uki
Uij δ kj mχk ∂β mνj
n
∂α mλµ
− ∂α mλi mµi
i
o
k
o
1 χλ µν n
+
Uji Uij
∂β mχν
∂α mλµ
− ∂α mλi mµi ∂β mχj mνj ,
j
i
2
expressão que pode ser interpretada de forma semelhante à da secção anterior.
2.4.4 Generalização
Ao longo do cálculo dos diferentes termos da expansão, o que fizemos foi expandir os
cumulantes em momentos estatísticos e, posteriormente, separar os momentos estatísticos
34
2.4 Expansão de altas temperaturas
em todas as possibilidades de os índices, somados sobre todas as partículas, serem a)
diferentes, caso em que as médias livres factorizam, ou b) iguais, dando origem a médias
de produtos de componentes. Em ambos os casos, o que fizemos é conhecido em Análise
Combinatória como partições de conjuntos (e.g. Kreher e Stinson [34]).
Cumulantes como partições
Em primeiro lugar, notamos que os cumulantes de HI que aparecem na expansão 2.25
podem ser genericamente escritos como (Fukushima [35])
h
n
X
X
D
(1)
(n)
HI ...HI
(1)
(j )
(−1)k−1 (k − 1)! HI ...HI 1
k=1 P (n,k)
i
E D
0
=
0
(j +1)
HI 1
(2.60)
(j )
...HI 2
E
0
D
(j
... HI k−1
+1)
(j )
...HI k
n
(1)
E
0
(n)
onde a soma sobre P (n, k) é a soma sobre todas as partições do conjunto HI ...HI
em k blocos. Alguns casos particulares são
[Hi ]0 = hHI i0
D
[Hi2 ]0 = Hi2
D
[Hi3 ]0 = Hi3
E
0
E
0
D
− 3 HI2
(2.61)
(2.62)
− hHI i20
E
0
o
(2.63)
hHI i0 + 2 hHI i30
onde já tiramos partido do facto de que os n Hamiltonianos no cumulante são todos
equivalentes.
Somas de momentos como partições
Quando, na exposição anterior, escrevemos explicitamente os Hamiltonianos que figuram
nos momentos estatísticos, encontrámos termos que se podem escrever genericamente
como
hHin i0 = −
1
2
n X
X
i1 ...i2n α1 ...αn
α
α
D
α
2n−1 2n
2n−1 α2n
Uiα11i2α2 ...Ui2n−1
µαi11 µαi22 ...µi2n−1
µi2n
i2n
E
0
(2.64)
e separámos estes momentos estatísticos contabilizando o número e o tipo de formas como
os índices de partícula i podem ser iguais ou diferentes. Para generalizar este cálculo,
vamos necessitar de introduzir algumas noções de Análise Combinatória.
35
2 Métodos analíticos
Relação com partições de conjuntos
O número de partições do conjunto {1, 2, ..., n} é conhecido como número de Bell, B(n), e
o número de partições desse conjunto em k subgrupos é o número de Stirling do segundo
tipo, S2 (n, k). Naturalmente,
B(n) =
n
X
S2 (n, k)
(2.65)
k=1
Como os Hamiltonianos no cumulante são todos equivalentes, podemos agrupar termos
do mesmo tipo (o mesmo número de Hamiltonianos em cada média). Os tipos possíveis
são dadas pelas partições inteiras de n, P (n), ou seja, as formas de escrever n como uma
soma de inteiros.
Para exemplificar, o conjunto {1, 2, 3, 4} tem B(4) = 15 partições diferentes, S2 (4, 1) =
1 com um bloco
{1, 2, 3, 4}
S2 (4, 2) = 7 com dois blocos
{{1, 2, 3}, {4}} {{1, 2, 4}, {3}} {{1, 3, 4}, {2}} {{2, 3, 4}, {1}}
{{1, 2}, {3, 4}} {{1, 3}, {2, 4}} {{1, 4}, {2, 3}}
S2 (4, 3) = 6 com três blocos
{{1, 2}, {3}, {4}} {{1, 3}, {2}, {4}} {{1, 4}, {2}, {3}}
{{2, 3}, {1}, {4}} {{2, 4}, {1}, {3}} {{3, 4}, {1}, {2}}
e S2 (4, 4) = 1 com quatro blocos
{{1}, {2}, {3}, {4}}.
De seguida, tentaremos sistematizar a contagem dos tipos de partições relevantes.
Tipos de partições
As partições do inteiro 4 são dadas por [4], [3, 1], [2, 2], [2, 1, 1] e [1, 1, 1, 1], que correspondem aos tipos de partições equivalentes que encontramos em cima. Note-se que existem
dois tipos de partições com dois blocos, mas com diferentes tamanhos de bloco, pelo que
o número de Stirling do segundo tipo não é o coeficiente que procuramos para enumerar
36
2.4 Expansão de altas temperaturas
partições equivalentes. O número de vezes que cada tipo aparece é dado pelos coeficientes
de di Bruno [36],
Cn (m1 , ..., mn ) =
n!
m1 !m2 !m3 ! · · · 1!m1 2!m2 3!m3 · ··
(2.66)
onde mi é o número de blocos de tamanho i. Ou seja, existem C4 (, 0, 0, 0, 4) = 1 partições
do tipo [4], C4 (1, 0, 1, 0) = 4 partições do tipo [3, 1], C4 (0, 2, 0, 0) = 3 partições do tipo
[2, 2], C4 (2, 1, 0, 0) = 6 do tipo [2, 1, 1] e C4 (4, 0, 0, 0) = 1 do tipo [1, 1, 1, 1].
Representação Diagramática
Existe uma representação de partições de conjunto que é particularmente útil no nosso
contexto. Trata-se de representar um conjunto de n elementos por um polígono regular
com n vértices onde elementos pertencentes ao mesmo bloco são unidos por um segmento de recta e blocos de apenas um elemento são representados por um ponto isolado
(figura 2.1). As Cn (m1 , ..., mn ) partições de cada tipo correspondem simplesmente a diferentes posições para os pontos e segmentos de recta, pelo que basta representar uma
delas.
Figura 2.1: Representação diagramática dos tipos de partições de um conjunto de 4
elementos.
No nosso caso, em que os momentos magnéticos vêm de Hamiltonianos de interacção,
existem restrições devido ao facto de não existir auto-interacção, ou seja, os índices de
momentos magnéticos provenientes do mesmo Hamiltoniano nunca podem ser iguais. Isto
quer dizer que num conjunto de n elementos provenientes de n/2 Hamiltonianos, não
podem existir blocos com mais de n/2 elementos, pelo que esses tipos podem ser eliminados
à partida (no exemplo, os tipos [4] e [3, 1] não são permitidos). Esta restrição significa
também que se k for um número ímpar, não são permitidas partições em ik e ik+1 apareçam
no mesmo bloco.
Usando esta representação, verificamos rapidamente que as partições de momentos
estatísticos de F2 na equação 2.50 podem ser apresentados como na figura 2.2
37
2 Métodos analíticos
Figura 2.2: Representação diagramática da correcção de segunda ordem à energia livre.
o que corresponde à forma separada da equação 2.52. Esta notação permite calcular
a expansão da energia livre de forma mais simples e, consequentemente chegar a termos
de ordem superior. No entanto, para calcular a magnetização e funções de correlação é
necessário introduzir uma notação para as derivadas em ordem a cada índice.
Uma forma de o fazer é marcar os pontos ou linhas com o índice em ordem aos quais
esse bloco foi derivado, e multiplicar o número de partições equivalentes pelo número de
blocos equivalentes a esse, como exemplificado na figura 2.3.
Figura 2.3: Representação diagramática das derivadas de tipos de partições num conjunto
de 4 elementos.
Usando este esquema, e sabendo que M α (2) = −∂α F2 , a correcção de segunda ordem à
magnetização pode ser representada como na figura 2.4, o que corresponde à equação 2.55.
Figura 2.4: Representação diagramática da correcção de segunda ordem à magnetização.
Quanto à correcção de segunda ordem à correlação, C αβ (2) =
representada como na figura 2.5, correspondente à equação 2.59.
1
∂ M α (2),
β β
pode ser
Figura 2.5: Representação diagramática da correcção de segunda ordem à correlação.
Esta representação, embora não substitua o cálculo analítico por completo, permite
verificar rapidamente quais as partições de conjuntos que contribuem para cada grandeza,
e o respectivo peso.
38
2.5 Conclusões
2.5 Conclusões
Neste capítulo, revimos noções elementares de Física Estatística à medida que apresentámos os modelos que servem de base ao estudo que se segue. Para melhor enquadrar
e controlar os resultados de métodos numéricos como os que veremos no capítulo seguinte, calculámos analiticamente médias termodinâmicas de um sistema sem interacções
e sucessivas correcções no âmbito de uma expansão de altas temperaturas.
39
2 Métodos analíticos
40
3 Métodos numéricos
No capítulo 2, vimos como podemos obter uma aproximação analítica para as médias
termodinâmicas
dµ O(µ) e−βH
hOi =
(3.1)
dµ e−βH
de um sistema com interacções arbitrárias, através de uma expansão de altas temperaturas. Neste capítulo, vamos explorar uma outra forma de obter estimativas de médias
termodinâmicas, desta feita para qualquer temperatura.
Se projectarmos na energia o integral sobre todas as configurações, definindo a média
microcanónica
dµ O(µ) δ(E − Eµ )
(3.2)
O(E) =
dµ δ(E − Eµ )
e a densidade de estados
n(E) = log S(E) =
dµ δ(E − Eµ )
(3.3)
podemos escrever a média termodinâmica como um integral sobre a energia em vez das
configurações
hOi =
=
dE O(E) n(E) e−βE
dE n(E) e−βE
(3.4)
dE O(E) eS(E)−βE
,
dE eS(E)−βE
e notar que o termo eS(E)−βE faz com que para cada temperatura, apenas configurações
dentro de um intervalo estreito de valores da energia tenham uma contribuição relevante
para o valor do integral (figura 3.1), o que sugere a utilização de um método de integração
de Monte Carlo com amostragem por importância.
41
3 Métodos numéricos
Figura 3.1: O peso estatístico das configurações em função da sua energia é aproximadamente uma gaussiana.
3.1 Revisão de métodos de Monte Carlo
Quando uma função é conhecida analiticamente, a sua integração numérica por Monte
Carlo com amostragem por importância é trivial. No entanto, no caso do cálculo de
médias termodinâmicas, não sabemos escrever a densidade de estados e as médias microcanónicas das observáveis de interesse como função da energia, pelo que é necessário gerar
as configurações do sistema para calcular as suas propriedades. A dificuldade reside, naturalmente, em como gerar configurações do sistema com uma determinada distribuição
de amostragem e em saber qual a distribuição que permite uma melhor convergência, na
prática, ao valor da função.
3.1.1 Cadeias de Markov
Para resolver a primeira questão, de gerar configurações com uma determinada distribuição assimptótica, utiliza-se o formalismo de cadeias de Markov (e.g. Newman e Barkema [37]), uma forma de, a partir de um determinado estado inicial i, gerar um estado
final f . A probabilidade de transição entre os dois estados designa-se por P (i → f ) e,
para se gerar uma cadeia de Markov, é necessário que estas probabilidades não variem
com o tempo e que dependam apenas das propriedades dos estados i e f . Naturalmente,
é necessário que se cumpra a condição de normalização
X
P (i → f ) = 1,
(3.5)
f
onde não é necessário que P (i → i) seja igual a zero, ou seja, pode existir uma probabilidade de o sistema ficar no mesmo estado. Este facto é normalmente utilizado na
construção dos algoritmos, e de facto, em geral, P (i → i) não será igual a zero.
Para garantir que esta forma de escolher estados converge para uma amostragem com
42
3.1 Revisão de métodos de Monte Carlo
a distribuição assimptótica desejada, é necessário garantir que as probabilidades de transição, por um lado, são tais que nenhum estado do sistema é completamente inacessível
(ergodicidade) e que, por outro lado, satisfazem a condição de equilíbrio detalhado,
pi P (i → f ) = pf P (f → i),
(3.6)
onde pi é a probabilidade assimptótica de encontrar o sistema no estado i (no ensemble
canónico, pi = e−βEi /Z).
3.1.2 Selecção e aceitação
Como vimos, não é necessário que a probabilidade de ficar no mesmo estado seja nula, o
que sugere introduzir alguma liberdade na escolha das probabilidades de transição, sem
violar a condição de equilíbrio detalhado, separando a transição em dois passos, onde
primeiro seleccionamos um novo estado e, de seguida, decidimos se o aceitamos ou não.
Vamos então escrever a probabilidade de transição como o produto
P (i → f ) = S(i → f )A(i → f )
(3.7)
onde S(i → f ) é a probabilidade de selecção, a probabilidade de propor um estado f a
partir do estado i, e A(i → f ) é a probabilidade de aceitação, a probabilidade com que
aceitamos essa proposta. A condição de equilíbrio detalhado (eq. 3.6) fixa apenas a razão
entre as probabilidades de transição
S(i → f )A(i → f )
pf
=
pi
S(f → i)A(f → i)
(3.8)
Dado que pretendemos amostrar o maior número de estados diferentes para calcular
os estimadores das grandezas relevantes, é conveniente garantir que a probabilidade de
aceitação não é demasiado baixa.
Algoritmo de Metropolis
O algoritmo de Metropolis [38] consiste em seleccionar os estados acessíveis com uma
probabilidade uniforme (figura 3.2) e fazer com que uma das probabilidades de aceitação
seja igual à unidade, nomeadamente a aceitação no sentido em que a energia diminui.
43
3 Métodos numéricos
Para tal, a probabilidade de aceitação é dada então por


1,
Ef 6 Ei
pf /pi ,
Ef > Ei
A(i → f ) = 
.
(3.9)
Esta escolha, apesar de garantir a condição de equilíbrio detalhado, pode levar a probabilidades de aceitação pouco satisfatórias, principalmente a baixas temperaturas e/ou
campos elevados, dado que a maior parte das propostas implica um custo energético demasiado elevado. O ideal para uma exploração eficiente do espaço de fase seria seleccionar
os estados propostos com a probabilidade de transição correcta, de tal forma que a probabilidade de aceitação fosse sempre igual à unidade. De seguida, vamos procurar formas
de propor configurações que nos permitam aproximar deste objectivo.
S
Figura 3.2: Representação esquemática da selecção uniforme.
Selecção Optimizada
Designando por Pf i a probabilidade de transição do estado i para o estado f
Pf i = P (i → f ) = Sf i Af i
(3.10)
onde Sf i é a probabilidade de selecção do estado f e Af i é a probabilidade de aceitação
da transição i → f , a condição de equilíbrio detalhado diz-nos que
Pf i
Af i Sf i
=
= e−β(Ef −Ei ) .
Pif
Aif Sif
(3.11)
Escrevendo a energia à custa dos campos locais H k (que podem incluir tanto o campo
externo como a interacção entre as partículas), e tendo em conta que a diferença de
energia se deve apenas a uma partícula k e que os campos nessa partícula dependem do
44
3.1 Revisão de métodos de Monte Carlo
momento dessa partícula (mk,f na configuração f e mk,i na configuração i), verificamos
que a diferença de energias que aparece na exponencial é dada por
Ef − Ei = −H k · (mk,f − mk,i ).
(3.12)
Quando a selecção é uniforme, é a probabilidade de aceitação que tem de satisfazer a
condição da eq. 3.11. Para melhorar a probabilidade de aceitação, é necessário incluir
parte desta diferença de energia na proposta da nova orientação do momento.
Energia de Zeeman na probabilidade de selecção
Vamos então separar o campo local H k no campo devido às interacções, H I,k , e o campo
externo, H ext ,
H k = H I,k + H ext
(3.13)
e escrever o produto escalar à custa do ângulo que o momento que estamos a trocar faz
com o campo externo,
Ef − Ei = −H I,k · (mk,f − mk,i ) − Hext (cos θk,f − cos θk,i ).
(3.14)
Se escolhermos a nova posição com uma probabilidade proporcional a eβHext cos θ , ou
seja,
Sf i ∝ eβHext cos θk,f
e
Sif ∝ eβHext cos θk,i
então a probabilidade de aceitação vai depender apenas das interacções entre partículas.
Af i
= eβH I,k ·(mk,f −mk,i )
Aif
(3.15)
Na prática, deixamos de propor configurações cujo custo energético devido ao campo
externo seria demasiado alto, o que melhora a taxa de aceitação a campos altos.
Energia de interacção na probabilidade de selecção
Para melhorar a taxa de aceitação a campos fracos e baixas temperaturas, podemos notar
que, nos sistemas planares de que nos vamos ocupar, e se o campo externo estiver também
45
3 Métodos numéricos
S
Figura 3.3: Representação esquemática da selecção optimizada tendo em conta o campo
externo.
k
aplicado no plano, H ext = H ext , então a componente dos campos locais perpendicular
k
⊥
ao plano, Hk⊥ = HI,k
, é muito pequena. Se a separarmos da componente planar H k =
k
H ext + H I,k , podemos escrever a diferença de energia de uma forma semelhante ao que
fizemos com o termo de Zeeman,
k
⊥
e
Ef − Ei = −Hk⊥ (m⊥
k,f − mk,i ) − Hk (cos θk,f − cos θ̃k,i )
(3.16)
onde agora θe é o ângulo que o momento faz com a componente planar do campo local
k
k
H k = H ext + H I,k .
Se, mais uma vez, escolhermos a nova orientação do momento com uma probabilidade
tal que
k
e
Sf i ∝ eβHk cos θk,f
e
k
Sif ∝ eβHk cos θ̃k,i
então a probabilidade de aceitação vai depender apenas da componente da energia de
interacção devida às componentes perpendiculares ao plano dos momentos,
⊥
⊥
⊥
Af i
= eβHk (mk,f −mk,i ) .
Aif
(3.17)
Como iremos ver mais à frente, esta simples modificação pode ter efeitos muito relevantes na convergência e precisão das simulações.
46
3.1 Revisão de métodos de Monte Carlo
S
~
~
Figura 3.4: Representação esquemática da selecção optimizada tendo em conta o campo
externo e a componente planar do campo dipolar local.
Números aleatórios com distribuição exponencial
Para implementar esta selecção optimizada, é necessário gerar os novos valores de u = cos θ
de acordo com a seguinte distribuição de probabilidade
P (u) = 1
eβHu
−1
eβHu du
=
2
βH
eβHu
.
sinh(βH)
(3.18)
Para gerar um número u com esta probabilidade, geramos um número aleatório r com
distribuição uniforme entre 0 e 1 e usamos a função cumulativa para inverter,
U
P (u)du =
r=
−1
U=
eβHU − eβH
2 sinh(βH)
1
log 2 sinh(βH)r + e−βH .
βH
(3.19)
(3.20)
No caso de u = cos θ, sendo θ o ângulo que o momento faz com o eixo x, então
mx = u,
my =
mz =
√
1 − u2 sin φ,
(3.21)
√
1 − u2 cos φ
e sendo θe o
onde o ângulo φ é gerado uniformemente entre 0 e 2π. No caso de u = cos θ,
ângulo que o momento faz com a componente planar do campo local, então as componentes
cartesianas mx,y,z são calculadas a partir de u e φ através de uma rotação adicional,
47
3 Métodos numéricos
descrita pela matriz


cos ω − sin ω 

sin ω cos ω
(3.22)
onde ω é o ângulo que a componente planar do campo local faz com o eixo x, definido por
cos ω = r
sin ω = r
x
Hext + HI,k
Hext +
x
HI,k
2
+
+
2
y
HI,k
y
HI,k
Hext +
x
HI,k
2
2
y
HI,k
,
(3.23)
.
(3.24)
Estes cálculos têm um custo computacional apreciável quando comparados com o caso
da selecção uniforme. No entanto, dado que a taxa de aceitação passa a ser próxima de 1,
este custo é desprezável quando comparado com o custo de actualização dos campos em
todas as partículas, que acaba sempre por dominar o perfil computacional, principalmente
em sistemas com interacções de longo alcance.
3.1.3 Distribuições de amostragem
Outro grau de liberdade à nossa disposição quando efectuamos simulações de Monte
Carlo é a escolha da distribuição de amostragem, mantendo a mesma conectividade entre
estados, permitindo formas de explorar o espaço de fase que podem ser mais eficientes em
determinadas situações.
Peso de Boltzmann
A escolha mais natural para a distribuição de amostragem consiste em usar uma probabilidade proporcional ao peso de Boltzmann dessa configuração, dado ser essa a importância
da sua contribuição para o integral. Usando a condição de equilíbrio detalhado, isto
equivale a dizer que
pf
P (i → f )
=
= e−β(Ef −Ei ) .
(3.25)
pi
P (f → i)
Neste caso, as médias termodinâmicas são dadas simplesmente pela média aritmética
das medições de cada observável a temperatura constante. Num cenário ideal, esta é a escolha mais eficiente, porque concentra o tempo de simulação precisamente nas regiões do
espaço de fase que mais contribuem para as médias termodinâmicas a uma determinada
temperatura. No entanto, na prática esta escolha pode limitar a convergência da amostragem para esta distribuição, que só está garantida assimptoticamente. Sendo assim, pode
48
3.1 Revisão de métodos de Monte Carlo
compensar modificar a distribuição de amostragem, se isso permitir uma convergência
mais rápida para o regime assimptótico.
Figura 3.5: Representação esquemática do histograma de estados visitados a temperatura
constante.
Peso inversamente proporcional à densidade de estados
É possível também escolher distribuições de amostragem alargadas (Torrie e Valleau [39]),
por comparação com o peso de Boltzmann, no sentido em que a frequência com que uma
determinada configuração ocorre na simulação não corresponde ao peso estatístico no
ensemble canónico. Para ultrapassar problemas de ergodicidade, ou seja, para evitar que a
simulação fique confinada a uma região do espaço de fase, uma escolha comum é usar como
distribuição de amostragem o inverso da densidade de estados n(E) = log S(E) (Lee [40],
Berg [41, 42], Wang e Landau [43]), de tal forma que a probabilidade de encontrar a
simulação numa determinada energia é uniforme. Neste caso, podemos escrever a condição
de equilíbrio detalhado como
P (i → f )
pf
=
= e−(S(Ef )−S(Ei )) .
pi
P (f → i)
(3.26)
Neste formalismo, a escolha do peso de Boltzmann pode ser vista como uma aproximação linear à entropia, S(E) = βE, com β constante, enquanto que a utilização do inverso
da densidade de estados corresponde à construção de uma função β(E) = dS(E)/d(E),
o que está na origem do termo “multicanónico” para designar esta abordagem. Neste
caso, as médias termodinâmicas são acumuladas em função da energia e repesadas para
cada uma das temperaturas desejadas. A aplicação deste método à simulação de redes
quadradas de dipolos pontuais pode ser encontrada na dissertação de Mestrado do autor
[44].
49
3 Métodos numéricos
Figura 3.6: Representação esquemática do histograma de estados visitados numa simulação multicanónica.
3.1.4 Parallel Tempering
O método de Parallel Tempering, com origem em diversas propostas como as de Swendsen
e Wang [45], Marinari e Parisi [46], Lyubartsev et al. [47] e Hukushima e Nemoto [48],
consiste em simular M réplicas independentes do sistema a temperaturas diferentes,
{T1 , T2 , ..., TM } e intercalar as cadeias de Markov de cada simulação com uma cadeia
de Markov adicional que consiste em propor a troca de duas temperaturas adjacentes, Ti
e Ti+1 , proposta esta que é aceite com uma probabilidade
P (Ei , Ti → Ei+1 , Ti+1 ) = min{1, exp(∆B∆E)},
(3.27)
onde ∆B = 1/Ti+1 − 1/Ti e ∆E = Ei+1 − Ei . Esta cadeia de Markov adicional introduz
um passeio aleatório das réplicas no espaço de temperaturas, sem quebrar a condição de
equilíbrio detalhado, o que permite uma exploração mais eficiente do espaço de fase e
consequente diminuição dos tempos de autocorrelação. Naturalmente, simular M réplicas
independentes do sistema custa M vezes mais tempo de computação do que uma única
réplica, mas em geral estamos interessados em obter resultados a temperaturas diferentes,
pelo que já iríamos ter de assumir pelo menos parte desse custo de qualquer forma.
Figura 3.7: Representação esquemática da escolha de temperaturas numa simulação de
Parallel Tempering.
Este esquema de simulação deixa ainda vários parâmetros em aberto que podem ser
50
3.1 Revisão de métodos de Monte Carlo
afinados para optimizar a relação entre a convergência das simulações e o custo computacional, nomeadamente o número de passos de Monte Carlo (MC) entre cada troca de
temperatura, a que vamos chamar passo Parallel Tempering (PT), o número de temperaturas e a sua distribuição.
Em relação ao número de passos MC entre cada passo PT, que vamos designar por
Nlocal , a escolha mais comum é Nlocal = 1, em unidades de passos de Monte Carlo. Bittner
et al. [49] sugeriram recentemente que medir o tempo de autocorrelação de simulações
canónicas τcan (β) a cada uma das temperaturas utilizadas e escolher Nlocal (β) de forma
proporcional a este tempo de autocorrelação permitiria reduzir drasticamente os tempos
de autocorrelação das simulações com Parallel Tempering, chegando mesmo a eliminar o
efeito de abrandamento crítico quando presente nas simulações canónicas. No entanto,
nas nossas implementações verificámos que não é esse o caso - efectivamente, esse aumento
de passos locais corresponde a uma diminuição de frequência de troca de temperaturas,
o que aumenta o tempo de autocorrelação, em vez de o diminuir. A única forma de
reproduzir os resultados desse trabalho foi medir o tempo de autocorrelação sem ter em
conta que foram efectuados Nlocal (β) passos adicionais, o que corresponde a medir o tempo
de autocorrelação em unidades de Nlocal (β), o que explica o aparente desaparecimento do
abrandamento crítico.
No que diz respeito às temperaturas, T1 e TM devem ser naturalmente as temperaturas
máxima e mínima em que estamos interessados, devendo a máxima situar-se acima da
temperatura crítica ou de bloqueio, caso exista. Dado que uma troca de temperaturas é
aceite com a probabilidade da equação 3.27, é importante que ∆β e ∆E sejam pequenos,
o que equivale a dizer que existe uma grande sobreposição entre as energias visitadas
por cada simulação. O esquema utilizado para definir estas temperaturas (Lopes [50]),
foi o de, começando na temperatura mais alta, medir as larguras σ dos histogramas
obtidos à temperatura actual e colocar a temperatura seguinte a βi+1 = βi + f /σ, sendo
o multiplicador f utilizado para controlar a taxa de aceitação de trocas PT. Um valor
de f demasiado elevado leva a que não exista sobreposição entre as distribuições e as
trocas nunca sejam aceites, enquanto que um valor demasiado pequeno tem um custo
computacional demasiado elevado, dado serem necessárias mais temperaturas para cobrir
o mesmo intervalo. Neste trabalho, foi utilizado um valor de f = 1, o que garantiu uma
taxa de aceitação próxima dos 50%, e um número de réplicas da ordem das centenas para
os intervalos considerados (note-se que o número de réplicas necessário cresce também com
√
o tamanho do sistema, sendo proporcional a N ). Katzgraber et al. [51] propuseram
um esquema dinâmico de alocação de temperaturas de forma a maximizar as “viagens de
ida-e-volta” de réplicas entre as temperaturas Ti e Ti+1 , com resultados interessantes a
nível de diminuição dos tempos de autocorrelação. No entanto, este método torna-se algo
51
3 Métodos numéricos
complexo, e não chegou a ser implementado neste trabalho.
3.1.5 Dinâmicas
Para além das probabilidades de selecção e aceitação e da escolha da distribuição de
amostragem, existem ainda diferentes formas de gerar um estado a partir de outro. O
método mais simples é alterar apenas um momento de cada vez, um método conhecido
como single spin flip. Conjugado com o algoritmo de Metropolis e o peso de Boltzmann
como distribuição de amostragem, verifica-se que este método leva a um abrandamento
da convergência das simulações à medida que nos aproximamos de uma temperatura
crítica, um efeito designado por critical slowing down. Para evitar este abrandamento,
foram propostas várias alternativas, como a utilização de algoritmos de cluster, tanto para
interacções de curto alcance (e.g. Swendsen e Wang [52], Wolff [53]), como, dentro de
certas aproximações, para interacções de longo alcance (e.g. [54]).
3.2 Observáveis
A primeira observável a ter em conta ao efectuar simulações de Monte Carlo é a energia
do sistema, dado ser esta que caracteriza o peso estatístico da configuração no cálculo
de médias termodinâmicas. A sua evolução deverá tender para uma distribuição estreita,
aproximadamente gaussiana, à volta de uma energia média que depende da temperatura,
com uma largura dada por
D
E
∆E 2 = E 2 − hEi2 .
(3.28)
As flutuações da energia estão relacionadas com o calor específico do sistema por momento magnético,
D E
dE
(3.29)
= β 2 E 2 − hEi2 /N,
dT
cujo comportamento térmico é utilizado na caracterização do regime crítico dos sistemas.
C=
Nos sistemas constituídos por momentos magnéticos, estaremos também interessados
na magnetização (ou de um parâmetro de ordem relacionado com a magnetização) média
dos sistemas, que iremos estudar em função da temperatura e/ou de um campo externo
aplicado. A magnetização terá uma evolução semelhante à energia, e as suas flutuações
são caracterizadas pela susceptibilidade magnética por momento magnético,
χ=
52
D
E
dM
= β M 2 − hM i2 /N.
dT
(3.30)
3.2 Observáveis
Nos sistemas com uma transição de fase, estas flutuações divergem na temperatura
crítica, mas os efeitos de tamanho finito tornam difícil a tarefa de a localizar, pelo que é
utilizada uma grandeza especificamente concebida para cancelar estes efeitos na temperatura crítica (e.g. Binder e Heermann [55]), o cumulante de Binder
2
hM 2 i
,
U =1−
3 hM 4 i
(3.31)
cujo valor a T = Tc deverá ser aproximadamente independente do tamanho do sistema.
Para cada grandeza A, são então armazenadas médias temporais de hAn i, com n =
1, 2, 4 e 8, para permitir o cálculo destas grandezas e das suas barras de erro. No entanto,
para este último efeito, estimar o erro na medição de cada grandeza, é necessário também medir a correlação, devida à utilização de uma cadeia de Markov, entre valores da
observável em diferentes instantes da evolução.
Para perceber de que forma esta auto-correlação afecta as barras de erro, temos de ter
em conta que aquilo a que temos chamado de valor médio é na verdade um estimador do
valor médio, calculado a partir de uma série temporal finita A1 ...At ...ANt . O estimador
do valor médio é dado por
Nt
1 X
Ā =
Ai
(3.32)
Nt i=1
e o valor médio do estimador é igual ao valor médio da população,
< Ā >=
Nt
1
1 X
< Ai >=
Nt < A >=< A > .
Nt i=1
Nt
(3.33)
No entanto, para calcular a variância do estimador,
var(Ā) = < Ā2 > − < Ā >2
(3.34)
= < Ā2 > − < A >2
é necessário ter em conta a forma como a série foi gerada. Para tal, começamos por notar
que o termo
< Ā2 > =
1 X
< Ati Atj >
Nt2 ti tj
(3.35)
pode ser escrito, separando os casos em que ti = tj , na forma



X
1 X
< Ā2 >= 2
< A2ti > +
< Ati Atj > .

Nt  ti
ti 6=tj
(3.36)
53
3 Métodos numéricos
Se todos os pontos da série fossem descorrelacionados, então para ti 6= tj ,
(3.37)
< Ati Atj >=< Ati >< Atj >=< A >2
e a variância seria dada por
o
1 n
2
2
2
2
2
N
<
A
>
+(N
−
N
)
<
A
>
−N
<
A
>
t
t
t
t
Nt2
o
1 n
=
< A2 > − < A >2 ,
Nt
var(Ā) =
(3.38)
como seria de esperar. No entanto, quando existe correlação entre os valores medidos a
tempos diferentes, a equação 3.37 não se verifica. Se definirmos a função de autocorrelação
temporal da grandeza A como
D
CA (ti , tj ) =
E
D
Ati Atj − hAti i Atj
hA2 i − hAi2
E
,
(3.39)
(3.40)
podemos substituir a equação 3.37 por
< Ati Atj >=< Ati >< Atj > +
D
E
A2 − hAi2 CA (ti , tj ).
Neste caso, a variância do estimador, tendo em conta a autocorrelação, será então dada
por



D E
i
X h
1 
2
2
2
2
2
2
var(Ā) =
N
<
A
>
+
hAi
+
A
−
hAi
C
(t
,
t
)
−
N
<
A
>
t
A i j
t

Nt2 
ti 6=tj


< A2 > − < A >2 X
1 
< A2 > − < A >2 +
CA (ti , tj )
=
Nt
Nt
ti 6=tj

(3.41)

< A2 > − < A >2 
1 X
=
1+
CA (ti , tj ) .
Nt
Nt ti 6=tj
Se assumirmos que a função de autocorrelação só depende da diferença entre os tempos
ti e tj ,
CA (ti , tj ) = CA (|ti − tj |)
(3.42)
e que o tamanho da série, Nt , é suficientemente grande para podermos desprezar efeitos
54
3.3 Ferramentas de computação científica
de fronteira, então podemos ainda simplificar para
Nt
X
< A2 > − < A >2
CA (t) .
1+2
var(Ā) ≈
Nt
t
!
(3.43)
O somatório da função de correlação que aparece nesta expressão é frequentemente
designado por tempo de autocorrelação integrado,
τA =
Nt
X
CA (t),
(3.44)
t
e é esta a grandeza que iremos medir nas simulações para poder estimar o erro na medição
de grandezas termodinâmicas. O fenómeno de critical slowing down, já referido, deve-se
à divergência destes tempos próximo de uma temperatura crítica, e é dependente da
dinâmica de simulação utilizada, como iremos ver.
3.3 Ferramentas de computação científica
Um dos aspectos mais importantes a ter em conta ao programar algoritmos de computação
científica é ter resultados reprodutíveis, o que não é tão fácil de garantir como possa
parecer à partida. Esta secção é uma tentativa de transmitir experiências, métodos e
ferramentas de engenharia de software que possam ajudar o leitor a evitar alguns desses
obstáculos. Tratando-se apenas de um levantamento de alguns aspectos que devem ser
tidos em conta, recomenda-se a consulta de textos mais detalhados (e.g. Katzgraber [56])
e da documentação de cada uma das ferramentas.
Software livre Software livre [57] é, em termos gerais, software que pode ser usado, estudado e modificado sem restrições. A utilização de software livre nas tarefas comuns
do dia a dia já é importante pelas garantias de liberdade de utilização incondicional
e de interoperabilidade com outros sistemas. No desenvolvimento de algoritmos de
computação científica, as liberdades de estudo, modificação e reutilização ganham
especial relevância. Por exemplo, as bibliotecas livres permitem reutilizar código
já existente para as componentes não originais dos algoritmos (e.g. álgebra linear,
transformadas de Fourier, etc.), sem perder a possibilidade de verificar a implementação e, se necessário, alterá-la.
Ambientes Integrados de Desenvolvimento Como é bem sabido, qualquer editor de
texto pode ser usado para escrever um programa de computador, pelo que não é
forçosa a utilização de um ambiente integrado de desenvolvimento (IDE). Frequen-
55
3 Métodos numéricos
Figura 3.8: Ferramentas de computação científica.
56
3.3 Ferramentas de computação científica
temente, estes ambientes são vistos como demasiado complexos, desnecessários e
até como um obstáculo. No entanto, à medida que aumenta a complexidade dos
códigos, a integração de um editor de texto com as ferramentas que iremos ver a
seguir tornam o processo de desenvolvimento mais rápido e ágil, pelo que a sua
utilização é aconselhada desde o início. Entre IDE’s livres que podem ser utilizados
para este fim, salientamos o Eclipse [58], Kdevelop [59] e Anjuta [60], dependendo
das preferências de ambiente de trabalho do programador. Um IDE permite a integração do editor com ferramentas como sistemas de controlo de versões, gestão
de tarefas e projectos, depuradores, profilers, gestão de dependências e geração de
documentação, ferramentas que passamos a descrever a seguir.
Controlo de versões Um sistema de controlo de versões é um conjunto de ferramentas
que permitem gerir, como o nome indica, diferentes versões de um código ou documento. Estas diferentes versões podem referir-se a diferentes momentos da evolução
do projecto, a diferentes cópias distribuídas por pessoas diferentes ou a diferentes
ramos de um tronco comum. Isto é feito através da utilização de cópias de trabalho,
onde as alterações são feitas, de repositórios, onde as diferentes revisões são armazenadas, e de um conjunto de algoritmos que permitem submeter, extrair, comparar
e combinar alterações. Para a reprodutibilidade de simulações, têm o contributo
fundamental de permitirem, a qualquer altura no futuro, extrair do repositório a
revisão exacta do código que gerou um determinado resultado (desde que seja feita
uma submissão ao repositório antes de todos as corridas de produção e registado o
número revisão junto com os resultados.). Entre os sistemas de controlo de versões
mais usados, salientamos o Subversion [61], Git [62] e Mercurial [63].
Gestão de projectos No contexto de desenvolvimento de software, estamos a referir-nos
a sistemas que permitem, tipicamente através de um interface web, a visualização da
informação constante de um sistema de controlo de versões e a sua integração com
tarefas, documentos, etc. Entre estes, salientamos o Redmine [64], pela facilidade de
gestão de múltiplos projectos, mas um outro sistema muito utilizado por projectos
de software livre é o Trac [65].
Depuração Um depurador (debugger) permite inspeccionar a execução de um programa
à procura de erros de programação ou lógica. Para além do conhecido gdb [66],
ferramentas como o Valgrind [67] são particularmente úteis na detecção de fugas de
memória.
Perfis de utilização A medição de perfis de utilização (Profiling) refere-se ao acto de
analisar que componentes do programa consomem mais tempo e/ou recursos, de
57
3 Métodos numéricos
forma a permitir concentrar a optimização nos aspectos que podem ter um impacto
mais relevante. O gprof [68] é a escolha tradicional, mas o valgrind, já referido, pode
ser executado com duas ferramentas adicionais, callgrind e cachegrind, que permitem
gerar essa informação de profiling. A visualização destes dados é particularmente
facilitada pelo kcachegrind [69].
Construção de software Por construção (building) referimo-nos à automatização do processo de compilação e gestão de dependências de um programa. Tipicamente, é
usada uma Makefile feita à mão para automatizar algumas destas tarefas. Ferramentas como as autotools [70] e scons [71] permitem tornar este processo mais
independente do ambiente (sistema operativo, arquitectura, etc.) onde o software
está a ser instalado.
Linguagens
A escolha de linguagem de programação é sempre muito pessoal e alvo até de discussões
acaloradas. No entanto, há distinções objectivas que se podem fazer quanto à tarefa a
desempenhar. Por exemplo, é consensual que para as componentes mais computacionalmente intensivas se deve usar uma linguagem compilada e de baixo nível como Fortran ou
C, e que para tarefas como análise de dados, integração de componentes e automatização
de tarefas, a utilização de linguagens interpretadas e de alto nível como Python ou Perl
poupa muito tempo ao programador.
Uma questão menos consensual poderá ser, por exemplo, a escolha entre C ou C++.
Dado que o autor tem feito um esforço de migração para C++, linguagem que foi usada
para as simulações descritas nesta dissertação, ficam aqui algumas razões que motivaram
a escolha:
• programação orientada por objectos - apesar de obrigar a uma nova forma de pensar
sobre os problemas, o resultado final é mais estruturado e pode facilitar a reutilização
de código para outras simulações. Permite também o encapsulamento de dados e
métodos em objectos, encapsulamento que facilita o desenvolvimento e pode até
constituir uma boa estratégia de paralelização.
• bibliotecas standard de algoritmos e estruturas de dados - não sendo o nosso objecto de estudo a ciência de computadores, e apesar de ser didático implementar os
algoritmos e estruturas de dados mais comuns, na maior parte dos casos é preferível confiar no trabalho de especialistas que se dedicaram a encontrar as melhores
implementações para cada caso. Através do uso de templates, bibliotecas como a
STL e a Boost implementam interfaces genéricos para estas estruturas e algoritmos
58
3.3 Ferramentas de computação científica
de dados com implementações optimizadas para cada caso particular, pelo que o
desempenho não é prejudicado pela generalidade.
• templates - uma característica da linguagem que permite a definição de funções ou
classes genéricas que posteriormente o compilador implementa para cada caso particular usado no programa. Por exemplo, o caso mais simples seria uma função que
pode actuar sobre qualquer tipo de variável (int, float, double, etc.). Através de um
template, a função pode ser definida para um tipo em abstracto e depois o compilador trata de gerar as funções correctas dependendo do tipo de variável efectivamente
usado no programa. Há outros casos que poderiam ser tratados sem templates mas
para os quais a utilização de templates é mais eficiente, por definir um determinado parâmetro em tempo de compilação, como por exemplo a dimensionalidade do
problema.
• referências - em C, é frequente o uso de pointers, apontadores para endereços de
memória, não só para permitir o acesso a variáveis entre funções mas também para
optimizações e gestão de memória, entre outras tarefas. Um “programador perfeito”
não terá qualquer problema com esta abordagem, mas na prática nada impede que
um pointer aponte para um endereço de memória errado ou até não alocado. A
linguagem C++ introduziu o conceito de referência, que, apesar de ser muitas vezes
confundido com um pointer mais seguro, é na verdade um alias, um atalho para
um determinado objecto. O que torna as referências mais seguras é que têm de se
referir sempre a um objecto (um pointer pode ser declarado nulo e só depois lhe
ser atribuído um endereço) e não podem ser alteradas. Para além disso, sendo um
atalho para o objecto, a sintaxe de utilização é a mesma da do objecto, ao contrário
dos pointers que obrigam a uma notação própria. Por estas razões, quando se aceder
a um objecto por referência, podemos simplesmente agir como quando estamos a
aceder directamente ao objecto.
Compiladores A Gnu Compiler Collection [72] é uma colecção de compiladores livres que
inclui compiladores de C, C++, Fortran, Java, etc., e é, na opinião do autor, aquela
que mais garantias oferece de performance e portabilidade entre diferentes arquitecturas e sistemas operativos. No entanto, se a performance numa determinada
arquitectura em particular for um aspecto crítico, existem compiladores optimizados, como o caso da Intel Compiler Collection [73], que apesar de não ser livre, é
gratuita para fins não comerciais.
59
3 Métodos numéricos
Bibliotecas
Dependendo das aplicações, existem inúmeras bibliotecas para cálculo científico que podem ser usadas para diminuir o tempo de programação e aumentar a fiabilidade e performance das simulações.
Uma biblioteca bastante completa, e até didática, escrita em C mas que também pode
ser usada em C++, entre outras, é a Gnu Scientific Library [74]. De especial relevância
para este trabalho foram as funcionalidades de geração de números aleatórios, com diferentes métodos de elevado desempenho que podem ser usados para garantir que os resultados
não são dependentes do gerador. No entanto, a GSL inclui também o cálculo de funções
especiais, algoritmos de minimização, álgebra linear, análise combinatória, transformadas
de Fourier, integração numérica, entre outras.
Para quem optar por C++, já referimos a Standard Template Library, que implementa
estruturas de dados como vectores, listas, conjuntos e mapas e algoritmos como pesquisa
e ordenação nestas estruturas de dados, e a Boost, que é na verdade um conjunto de
bibliotecas tão diversas como o manuseamento de parâmetros e opções de um programa,
arrays multidimensionais, mapas bidireccionais, interacção com o sistema de ficheiros
de forma portável, paralelização e ligação com outras linguagens de programação como
Python.
No que diz respeito a arrays multidimensionais, vale a pena referir ainda uma outra
biblioteca, a Blitz++, que, alegadamente sem perda de performance, permite funcionalidades de alto nível como operações com notação vectorial ou tensorial, extracção de
componentes, entre outras. No entanto, o autor optou pela Boost::multi_array para este
efeito, principalmente pelas garantias de continuidade e eventual integração no standard
que a Boost oferece.
Mais uma vez, existem bibliotecas fornecidas pelos fabricantes de processadores optimizadas para cada arquitectura em particular. É o caso da Intel Math Kernel Library [75] e
da AMD Core Math Library [76] (não livres), que implementam versões optimizadas de
álgebra linear, transformadas de Fourier e geração de números aleatórios.
Paralelização
Apesar de todas as optimizações introduzidas por compiladores, bibliotecas e pelo próprio programador, é frequente encontrar problemas de computação científica para os quais
qualquer computador, sozinho, demoraria demasiado tempo a produzir os resultados esperados. É este tipicamente o caso das simulações de Monte Carlo, mas também de
algoritmos de dinâmica molecular, estrutura electrónica, dinâmica de fluidos, investigação operacional e até algoritmos que podem constituir alguma surpresa, como o cálculo
60
3.3 Ferramentas de computação científica
do preço de derivados financeiros.
Quando um computador não é suficiente, a solução passa por paralelizar o algoritmo,
ou seja, encontrar uma forma de dividir o problema em problemas mais pequenos que
possam ser enviados para diferentes processadores do mesmo ou de outros computadores
e implementar essa divisão. Existem diferentes formas de paralelizar (e.g., a nível das
tarefas ou a nível dos dados), e nalguns algoritmos a paralelização é mais eficiente do que
noutros, devido à quantidade de comunicação necessária em cada etapa e/ou tempos de
espera entre tarefas.
No caso mais simples, será necessário correr o algoritmo para muitos parâmetros independentes, caso em que bastará distribuir o cálculo por diferentes máquinas, o que não é
normalmente considerado uma paralelização. Pode acontecer também que haja um conjunto de tarefas que podem ser realizadas de forma independente, necessitando apenas de
comunicar alguma informação no fim dessas tarefas. Felizmente para nós, as simulações
de Monte Carlo cabem neste domínio e são facilmente paralelizáveis, e ainda mais no
caso de ser necessário fazer médias sobre a desordem, no sentido que se pode acumular
estatísticas durante muito tempo sem qualquer comunicação e pontualmente comunicar
médias da evolução do sistema durante esse tempo.
O caso em que a paralelização pode implicar mais comunicação é o caso da decomposição
dos dados em domínios, como nos algoritmos de álgebra linear, elementos finitos, etc.
Neste caso, a eficácia da paralelização depende em grande medida do tamanho dos dados
(quanto maiores, mais eficiente a paralelização).
Independentemente do tipo de decomposição, a implementação depende da arquitectura
em termos de acesso à memória, ou seja, se todos os processadores acedem directamente à
mesma memória (memória partilhada) ou se é necessário comunicar por rede para aceder
a dados que estão a ser executados noutro processador (memória distribuída). Neste
último caso, terá de ser usada uma biblioteca de passagem de mensagens (MPI [77, 78]),
enquanto que no caso de memória partilhada tanto pode ser usada uma implementação de
MPI como um método mais simples, a nível do compilador, designada por OpenMP [79].
Um desenvolvimento recente nesta área é a utilização de processadores gráficos (GPUs)
que exploraram a paralelização a nível dos dados para obter, para operações particularmente adequadas a este modelo, resultados comparáveis aos de clusters de computadores.
A Nvidia e a AMD / ATI disponibilizam kits de desenvolvimento para as suas plataformas, Cuda [80] e Stream [81], respectivamente, que permitem tirar partido deste tipo de
processadores. Começa também a ser possível encontrar bibliotecas que facilitam este
uso, nalguns casos com interfaces standard, como o da Lapack ou FFTs, de tal forma que
pode não ser necessária nenhuma alteração ao código para que a operação seja executada
nos processadores gráficos. É sem dúvida uma evolução a seguir na área de computação
61
3 Métodos numéricos
de elevado desempenho.
Clusters
Depois de paralelizado o algoritmo, é necessário submetê-lo ao maior número de computadores/processadores possível, se quisermos tirar partido dessa paralelização. Isto pode ser
feito simplesmente ligando vários computadores em rede, mas a complexidade de implementação e gestão deste cluster torna-se rapidamente um obstáculo, gastando os recursos
e tempo cuja optimização era o objectivo da paralelização em primeiro lugar.
Felizmente, existem diversas ferramentas que tornam esta tarefa muito mais simples,
reduzindo essa complexidade a pouco mais do que a de uma única máquina, pelo menos
no que diz respeito ao software. Os problemas que podem surgir a nível de hardware serão
sempre proporcionais ao número de componentes do sistema.
Um cluster pode ser implementado de muitas formas diferentes. As opções que o autor
tem seguido, passam por usar uma distribuição de Linux Enterprise como o CentOS [82],
pelo tempo de vida mais longo que proporcionam, e um sistema de aprovisionamento como
o Perceus [83] (uma alternativa semelhante seria o XCat [84]), que permite a distribuição,
por rede, de imagens stateless do sistema operativo a partir do master para os nodos de
cada vez que estes arrancam. Isto faz com que seja necessário gerir apenas um computador,
o master do cluster, e a imagem que é enviada para os nodos. Um upgrade a todos os
nodos pode ser tão simples quanto mudar a imagem e reiniciar os nodos.
Para submeter trabalhos de computação ao cluster, usa-se um gestor de recursos e
agendamento como o Sun Grid Engine [85] ou o Torque [86] com o Maui [87], que recebem
as submissões e tratam de as prioritizar, verificar quais os nodos disponíveis e enviar os
trabalhos para execução, ao mesmo tempo que permitem a contabilização e relatórios de
utilização.
Por fim, são necessárias ferramentas que permitam monitorizar o estado do cluster, a
nível da ocupação e disponibilidade dos recursos, tarefa esta que pode ser levada a cabo
com o Ganglia [88], o Nagios [89] e/ou o MRTG [90].
Formatos de dados
Um dos problemas com a reprodutibilidade na computação científica é o facto de, frequentemente, os ficheiros de dados resultantes das simulações não conterem informação sobre
as condições em que foram gerados, nem serem suficientemente descritivos para serem
entendidos por outras pessoas ou mesmo pelo seu criador algum tempo mais tarde. Uma
forma simples, e comum, de atenuar este problema é guardar os parâmetros de corrida
62
3.3 Ferramentas de computação científica
junto com os ficheiros de dados, e anexar alguns meta-dados aos ficheiros, o que faz sentido. Algo menos frequente mas que o autor recomenda vigorosamente é que seja incluída
também a revisão do código que gerou os dados (assumindo que se está a usar um sistema
de controlo de versões e que se actualiza o repositório antes de uma corrida de produção),
algo que pode ser feito de forma automática em tempo de compilação ou execução.
Existem também formatos de dados estruturados e descritivos especialmente concebidos
para armazenar e gerir grandes quantidades de dados numéricos (e.g. campos vectoriais),
como o NetCDF [91] ou o HDF [92]. A utilização destes formatos standard, além de tornar
o armazenamento e manuseamento mais eficiente, permite que os dados possam ser lidos
por outras aplicações de visualização ou análise de dados, dado que são auto-descritivos.
Visualização e Análise de Dados
No caso mais simples, o output de uma simulação pode ser simplesmente uma função
unidimensional armazenada na forma de duas colunas num ficheiro de texto, que pode
facilmente ser visualizado e analisado com uma aplicação como o xmgrace [93] ou o gnuplot [94]. Para muitos trabalhos científicos, isto é mais do que suficiente. Nos outros
casos, podem ser usadas bibliotecas e aplicações para a visualização campos escalares
ou vectoriais a duas dimensões (e.g., matplotlib [95]), a três dimensões (e.g. vtk [96],
mayavi [97]), ferramentas estas que produzem gráficos vectoriais que podem ser manipulados pelo inkscape [98]. Estas bibliotecas podem ser integradas com outras bibliotecas
de análise de dados como o SciPy [99][100].
Álgebra Computacional
Uma aplicação de uma utilidade extraordinária na computação científica (e não só) é o
Mathematica [101], da Wolfram, que provavelmente só tem um defeito - não é livre. Uma
alternativa livre poderá ser, para algumas aplicações, o sistema Sage [102].
Documentação
Para além dos comentários que o programador insere no código, ferramentas como o doxygen [103] permitem gerar documentação (html ou latex) com toda a estrutura do código,
incluindo parâmetros de funções, ligações entre classes e ficheiros, perfis de chamada de
funções, etc. Esta documentação pode ser gerada exclusivamente a partir do código,
mas é mais completa e adequada se o programador introduzir comentários num formato
standard que estas ferramentas entendem.
63
3 Métodos numéricos
Este capítulo não ficaria completo sem uma referência à aplicação usada para escrever
artigos e esta dissertação, o Lyx [104]. Trata-se de um frontend para latex cujo objectivo
é permitir aos autores concentrarem-se no conteúdo do que escrever e não na forma.
A recolha e gestão de referências bibliográficas pode ser feita com o Zotero [105], uma
extensão para o browser Firefox que recolhe referências a partir de páginas web com
informação bibliográfica (editoras, revistas, catálogos de biblioteca) e as exporta para o
formato bibtex.
3.4 Conclusões
Neste capítulo, apresentámos formas numéricas de calcular médias termodinâmicas de sistemas constituídos por momentos magnéticos clássicos, nomeadamente através de métodos de Monte Carlo. Foram introduzidas algumas das dificuldades típicas de convergência
das simulações para o regime assimptótico, e os métodos e ferramentas utilizados neste
trabalho para ultrapassar essas dificuldades. Como iremos ver de seguida, a eficácia destes métodos depende bastante do sistema a simular, mas é sempre preferível à utilização
dos métodos canónicos. Nos dois capítulos seguintes, são apresentados os resultados de
simulações para redes tanto ordenadas como desordenadas de dipolos pontuais.
64
4 Redes regulares de dipolos pontuais
Neste capítulo e no seguinte, iremos aplicar os métodos analíticos e numéricos apresentados nos capítulos anteriores para estudar o comportamento termodinâmico de sistemas
constituídos por dipolos pontuais, com três componentes cartesianas (momentos de Heisenberg), cuja interacção é descrita pelo tensor
Uijαβ
3uαij uβij − δαβ
,
= δ ij
3
rij
(4.1)
onde δ ij = (1 − δij ), uij = r ij /rij , e r ij = r i − r j .
Neste capítulo, iremos considerar redes regulares de dipolos pontuais, na sequência dos
trabalhos, tanto analíticos e computacionais, que resumimos de seguida.
Estado fundamental
Um dos resultados mais importantes de sistemas regulares constituídos por dipolos pontuais foi obtido usando o método de Luttinger e Tisza [30] para calcular a energia do
estado fundamental de uma rede cúbica ou quadrada de dipolos pontuais. O principal
resultado de Luttinger e Tisza foi o de que o estado antiferromagnético tem uma energia
mais baixa do que o estado ferromagnético.
Belobrov et al. [106, 107] mostraram que o estado fundamental é na verdade continuamente degenerado, tanto na rede quadrada como na rede cúbica. Brankov e Danchev
[108] generalizaram esta análise a duas dimensões para uma rede rômbica arbitrária, verificando que, dos estados continuamente degenerados na rede quadrada, o levantamento
da simetria da rede selecciona o estado antiferromagnético, mas também que o estado
fundamental pode passar a ferromagnético, como no caso da rede triangular. Zimmerman
et al. [109] consideraram o caso da rede de favo de mel, onde voltaram a encontrar um
estado fundamental continuamente degenerado e sem magnetização total.
As condições para o levantamento da degenerescência foram analisadas por Prakash e
Henley [110], propondo o mecanismo de “ordem devido a desordem”, em que a temperatura e pequenas desordem posicionais ou diluições da rede seleccionam estados particulares
do contínuo de estados degenerados.
65
4 Redes regulares de dipolos pontuais
Transição de Fase
Apesar de se verificar que o estado fundamental de um sistema puramente dipolar pode
ser ordenado, não é claro se é possível existir uma transição de fase a uma temperatura
finita. A uma dimensão, o trabalho de Dyson [111] excluiu essa hipótese, e a duas dimensões o teorema de Mermin-Wagner [112] diz-nos que, para sistemas de Heisenberg
com interacções isotrópicas e de curto alcance, não pode existir magnetização espontânea
a temperatura finita. No entanto, como vimos, a interacção dipolar não satisfaz estes
critérios, além de que o estado fundamental pode ter magnetização nula.
Bruno [113] generalizou estes resultados, excluindo ordem magnética a temperatura
finita para sistemas Heisenberg ou XY com interacções de longo alcance a decair como
r−α desde que o expoente α seja superior a 2d, sendo d a dimensionalidade do sistema.
Para a interacção dipolar, temos α = 3, pelo que este resultado exclui uma transição
de fase a temperatura finita em uma dimensão mas não em duas. Lax [114] estudou
uma rede de dipolos usando o modelo esférico, ou seja, substituindo a condição de todos
os momentos terem magnitude fixa pela condição de fixar a soma dos quadrados dos
momentos, e verificou que neste modelo não existe transição nem a uma nem a duas
dimensões, existindo uma transição para um estado antiferromagnético a três dimensões.
Simulações de Monte Carlo
S. Romano efectuou simulações de Monte Carlo de sistemas puramente dipolares [115,
116], que verificou as previsões teóricas, seguido de Rastelli et al. [117, 118, 119, 120,
121, 122, 123, 124, 125] e De ’Bell et al. [126, 127], que estudaram diversas variantes
de sistemas dipolares em redes regulares. No entanto, estas simulações foram feitas com
momentos XY, que podem não descrever correctamente as excitações a temperatura finita,
e usando o algoritmo de Metropolis, cujos problemas de convergência são particularmente
limitadores em sistemas com interacção dipolar, como iremos ver de seguida.
Muito recentemente, Tomita [128] efectuou simulações em redes planares de dipolos
pontuais do tipo Heisenberg, usando uma combinação de Parallel Tempering com um
novo algoritmo de Monte Carlo para interacções de longo alcance proposto por Fukui e
Todo [54]. Seria interessante comparar os tempos de autocorrelação obtidos com este método e os resultados apresentados de seguida, mas a publicação não fornece esses dados. A
análise do comportamento crítico da rede quadrada de dipolos pontuais será directamente
comparada com esse trabalho.
66
4.1 Rede quadrada
4.1 Rede quadrada
Como referido em cima, o estado fundamental de uma rede quadrada (infinita) de dipolos
pontuais é continuamente degenerado, em que a orientação de um momento no ponto
(nx , ny ) da rede pode ser descrito por
mnx ,ny = ((−1)ny cos θ, (−1)nx sin θ) ,
(4.2)
onde θ é um ângulo arbitrário. A forma deste estado, conhecido como “estado de microvórtice”, está esquematizada na figura 4.1.
Figura 4.1: Estado de micro-vórtice e os casos particulares de Luttinger e Tisza.
Quando θ é um múltiplo de π/2, encontramos um estado antiferromagnético por colunas, com os momentos magnéticos orientados ao longo dos eixos da rede, embora qualquer
ângulo θ produza um estado com a mesma energia.
A magnetização total deste estado é nula, apesar do estado fundamental ser ordenado,
pelo que é necessário encontrar um outro parâmetro que descreva o estabelecimento desta
ordem. Por analogia com um antiferromagnete comum, é possível construir uma “magnetização alternada” (staggered), da seguinte forma:
1
mst =
N
v
u
u X
u
t
nx ,ny
2
2

(−1)ny mxnx ,ny  + 
X
(−1)nx mynx ,ny  ,
(4.3)
nx ,ny
onde mαnx ,ny é a componente cartesiana α do momento magnético na posição (nx , ny ) da
rede. Este parâmetro de ordem aplica-se quer a magnetização tenha apenas duas componentes cartesianas, no plano da rede (ditos rotores planares ou momentos clássicos XY)
como três componentes (momentos clássicos de Heisenberg). No que se segue consideramos sempre momentos com três componentes cartesianas.
67
4 Redes regulares de dipolos pontuais
4.1.1 Condições de fronteira
Este estado de micro-vórtice é o estado fundamental de uma rede quadrada de dipolos
pontuais apenas no limite de amostra infinita. Se considerarmos uma amostra finita,
quadrada, com condições de fronteira abertas, a degenerescência é levantada. Em vez
de um ângulo de micro-vórtice arbitrário mas uniforme em toda a amostra encontramos
uma variação contínua deste ângulo, desde π/4 no centro da amostra a múltiplos de π/2
na periferia, de forma a alinhar os momentos com a fronteira da amostra (Costa [44],
Matsushita [129] e Sugano [130]).
Para eliminar este efeito, é necessário introduzir condições de fronteira periódicas, através de réplicas periódicas da amostra de simulação que preenchem todo o espaço. Em
sistemas com interacções de curto alcance, considera-se apenas a imagem mais próxima,
o que é equivalente a “dobrar” a amostra num toro a d + 1 dimensões.
Figura 4.2: Representação esquemática das condições de fronteira periódicas.
No nosso caso, com interacções de longo alcance, pode ser necessário considerar a
interacção com réplicas mais afastadas, até a distância tornar a interacção desprezável.
Na prática, como a configuração posicional e magnética das réplicas é igual à da amostra
de simulação, basta substituir a matriz de interacção da equação 4.1 por uma matriz de
interacção efectiva somada sobre todas as réplicas (figura 4.2), ou seja,
Uijαβ
=
X
n
68
3ũαij ũβij − δαβ
δ ij
,
3
r̃ij
(4.4)
4.1 Rede quadrada
onde n é a posição de cada réplica, ũij = r̃ ij /r̃ij , e r̃ ij = r ij + n.
4.1.2 Simulações Metropolis
Utilizando simulações de Metropolis convencionais, a diferentes temperaturas, para estudar o comportamento deste parâmetro em redes de N = L × L partículas, com condições
de fronteira periódicas, verificamos que de facto parece existir uma transição de fase a
temperatura finita para este estado ordenado. As figuras 4.3 e 4.4 mostram o comportamento, em função da temperatura e do tamanho L, das médias canónicas do parâmetro
de ordem mst , calor específico C, susceptibilidade χst e cumulante de Binder Ust , como
definidos na secção 3.2.
1
b)
a)
L=8
L=16
L=24
L=8
L=16
L=24
2
0.6
C
< mst >
0.8
1
0.4
0.5
0.2
0
0.125
0.25
0.5
T
1
2
0.25
0.5
1
2
T
Figura 4.3: Parâmetro de ordem (a) e calor específico (b) em função da temperatura,
obtido com simulações Metropolis convencionais e condições de fronteira periódicas e 107
passos de Monte Carlo.
Estes resultados são coerentes com a existência de uma transição de fase a uma temperatura próxima de T = 0.6 µ2 /(kb a3 ). No entanto, para caracterizar o comportamento
crítico do sistema, é necessário quantificar a precisão das simulações e simular sistemas
maiores. Como vimos no capítulo 3, a precisão com que estimamos uma determinada
grandeza termodinâmica está relacionada com o tempo de autocorrelação integrado dessa
grandeza na série temporal gerada pela cadeia de Markov, definido pela equação 3.44.
A medição deste tempo de autocorrelação em corridas finitas de simulações de Monte
Carlo levanta alguns problemas técnicos. O comportamento assimptótico da função CA (t),
para tempos muito longos, deverá ser um decaimento exponencial, mas, por um lado, este
comportamento poderá aparecer muito tarde e, por outro, o erro estatístico na medição
da correlação faz com que esta nunca seja de facto nula, mas sim que flutue à volta de
69
4 Redes regulares de dipolos pontuais
100
0.7
a)
b)
L=8
L=16
L=24
L=8
L=16
L=24
0.6
χst
Ust
10
0.5
1
0.4
0.1
0.25
0.5
1
2
0.3
0.125
0.25
T
0.5
1
2
T
Figura 4.4: Susceptibilidade e cumulante de Binder do parâmetro de ordem em função da
temperatura, obtida com Simulações Metropolis convencionais e condições de fronteira
periódicas.
zero com uma variância da ordem de 1/Nt , em que Nt é o número de medições para cada
t, pelo que o somatório não converge como seria de esperar (e.g. Madras e Sokal [131]).
O método que adoptámos para medir o tempo de autocorrelação foi calcular o somatório
e o desvio padrão do somatório em função do cut-off tmax ,
tX
max
CA (t)
(4.5)
< τA2 (tmax ) > − < τA (tmax ) >2
(4.6)
τA (tmax ) =
t=0
e
δτA (tmax ) =
q
e fixar tmax quando o erro relativo passa a ser superior a 1% (figura 4.5).
Os tempos de autocorrelação da energia e do parâmetro de ordem, estimados com este
método, para simulações de Monte Carlo convencionais da rede quadrada dipolos pontuais,
são apresentados na figura 4.6. Note-se que, mesmo para estes tamanhos relativamente
pequenos, os tempos de autocorrelação chegam, na região crítica, a valores da ordem
de 100 passos de Monte Carlo na energia e da ordem de 1000 passos de Monte Carlo
no parâmetro de ordem, e que nessa região os tempos aumentam com o tamanho do
sistema. Por outro lado, seria de esperar que esse comportamento desaparecesse longe da
região crítica, mas tal não acontece a temperaturas baixas, voltando mesmo a crescer e
eventualmente divergir à medida que a temperatura tende para zero, de uma forma que
parece independente do tamanho do sistema.
70
4.1 Rede quadrada
a)
b)
8
τ(tmax)
100*δτ(tmax)
τ
(MCSteps)
(MCSteps)
100
4
2
1
1
100
10000
1
1
1e+06
100
tmax (MCSteps)
10000
1e+06
tmax (MCSteps)
Figura 4.5: Exemplo de cálculo do tempo de autocorrelação em duas situações distintas.
No painel da esquerda, o tempo de autocorrelação é baixo mas é importante truncar o
somatório antes de o erro ser demasiado grande. No painel da direita, o regime exponencial
só aparece para tempos muito longos.
1000
10000
a)
L=8
L=16
L=24
L=8
L=16
L=24
1000
< τE >
< τM >
100
b)
100
10
10
1
0.25
0.5
T
1
2
1
0.25
0.5
1
2
T
Figura 4.6: Tempos de autocorrelação da energia (à esquerda) e parâmetro de ordem
(à direita) em função da temperatura e do tamanho do sistema, usando o algoritmo de
Metropolis convencional.
Se analisarmos a taxa de aceitação em função da temperatura, verificamos que ela
tende rapidamente para zero à medida que baixamos a temperatura, independentemente
do tamanho do sistema, pelo que deverá estar relacionada com o comportamento a baixas
temperaturas do tempo de autocorrelação. Na secção que se segue, vamos estudar o efeito
da aplicação do método de selecção optimizada apresentado no capítulo 3 nos tempos de
autocorrelação das simulações.
71
4 Redes regulares de dipolos pontuais
4.1.3 Taxa de aceitação optimizada
A figura 4.7 mostra uma secção da evolução da energia, para uma temperatura baixa,
obtidas usando o método de selecção uniforme e o método de selecção optimizada apresentado na secção 3.1.2, expondo a diferença entre as duas dinâmicas. Para o mesmo
intervalo de tempo (passos de Monte Carlo), as flutuações com a selecção optimizada
são muito superiores, permitindo uma melhor amostragem do espaço de fase. Este facto
só por si deverá contribuir para uma diminuição dos tempos de autocorrelação, quando
comparados com os medidos na secção anterior.
-2.3
-2.3
E/N
b)
E/N
a)
-2.4
-2.4
-2.5
100250
100500
100750
101000
-2.5
100250
t (MCSteps)
100500
100750
101000
t (MCSteps)
Figura 4.7: Comparação entre dinâmicas na rede quadrada - à esquerda, selecção uniforme, e à direita, selecção optimizada.
10000
a)
0.8
b)
S. Optimizada
L=8
L=16
L=24
S. Uniforme
L=8
L=16
L=24
1000
< τM >
<A>
0.6
0.4
100
10
0.2
S. Optimizada
S. Uniforme
0
0.125
0.25
0.5
T
1
2
1
0.125
0.25
0.5
1
2
T
Figura 4.8: Comparação de taxa de aceitação e tempo de autocorrelação do parâmetro
de ordem em função da temperatura.
72
4.1 Rede quadrada
Na figura 4.8a, é apresentada a melhoria em termos de taxa de aceitação devida à
utilização do método de selecção optimizada. Para além de aumentar nominalmente a taxa
de aceitação, esta deixa de diminuir com a temperatura, tendendo para um valor constante
superior a 80%. Note-se que tanto no caso de selecção uniforme como optimizada, a taxa
de aceitação praticamente não depende do tamanho do sistema. A figura 4.8b exibe o
efeito desta optimização nos tempos de autocorrelação do parâmetro de ordem, que são
uma ordem de grandeza mais baixos na região de baixas temperaturas.
Ainda assim, na região crítica o efeito desta optimização não é tão relevante, pelo que
são necessários outros métodos para diminuir os tempos de autocorrelação nesta região.
4.1.4 Parallel Tempering
O método de Parallel Tempering, descrito na secção 3.1.4, diminui ainda mais os tempos
de autocorrelação, ao permitir a mistura de cadeias de Markov a temperaturas diferentes.
A figura 4.9 mostra os tempos de autocorrelação na energia e no parâmetro de ordem com
as diferentes variantes descritas.
Metropolis
S. Optimizada (SO)
Parallel Tempering (PT)
PT + SO
1000
a)
10000
b)
1000
< τE >
< τM >
100
100
10
10
1
0.125
0.25
0.5
T
1
2
1
0.125
0.25
0.5
1
2
T
Figura 4.9: Comparação entre as diferentes dinâmicas, para L=24.
Se compararmos o método de selecção optimizada (linhas vermelhas na figura 4.9), com
o método de Parallel Tempering (linhas verdes), verificamos que este último leva a tempos
de autocorrelação mais baixos na região crítica, mas à medida que a temperatura diminui
é a selecção optimizada que produz melhores resultados. A combinação dos dois métodos
(linhas azuis) leva aos melhores resultados a qualquer temperatura, com diminuições dos
tempos de autocorrelação que podem chegar a duas ordens de grandeza, pelo que será
este o método a utilizar para estudar o comportamento crítico do sistema.
73
4 Redes regulares de dipolos pontuais
4.1.5 Expoentes críticos
Para estimar os expoentes críticos, é necessário, em primeiro lugar, determinar a temperatura crítica, Tc . Para tal, usaremos o cumulante de Binder
2
hm2st i
,
U =1−
3 hm4st i
(4.7)
cujo valor em T = Tc deverá ser aproximadamente independente do tamanho do sistema,
como vimos no capítulo anterior. Na figura 4.10, podemos encontrar o cumulante de
Binder em função da temperatura para os tamanhos L = 16 a 64, e, na magnificação,
desprezando os tamanhos mais pequenos, verifica-se que a Tc é aproximadamente 0.58.
L=16
L=24
L=32
L=40
L=48
L=56
L=64
0.7
0.6
0.68
U
0.66
0.5
Tc ~ 0.58
0.64
0.62
0.4
0.6
0.58
0.55 0.56 0.57 0.58 0.59 0.6 0.61
0.3
0.3
T
0.4
0.6
0.5
0.7
T
Figura 4.10: Cumulante de Binder em função da temperatura para redes quadradas de
dipolos pontuais com L = 16 a 64.
Se definirmos uma temperatura reduzida como uma medida da distância à temperatura
crítica
T − Tc
(4.8)
t=
Tc
então podemos definir os expoentes críticos ν, γ, α e β através do comportamento, na
região crítica e no limite termodinâmico, do comprimento de correlação
ξ ∼ |t|−ν ,
(4.9)
susceptibilidade
γ
χ ∼ |t|−γ = ξ ν ,
74
(4.10)
4.1 Rede quadrada
calor específico
C ∼ |t|−α = ξ ν
α
(4.11)
β
(4.12)
e magnetização
m ∼ |t|β = ξ − ν .
Num sistema finito, o comprimento de correlação nunca chega de facto a divergir, não
podendo ser maior do que o tamanho L do sistema, pelo que o comportamento destas
grandezas na região em que ξ seria maior do que L depende do tamanho do sistema (e.g.
Fisher e Barber [132]). Isto pode ser expresso à custa de funções de scaling, tal que, para
um função termodinâmica F (t), a sua dependência com o tamanho do sistema é dada por
ρ
1
F (t, L) = L ν F̃ (L ν t),
(4.13)
onde a função F̃ , adimensional e de uma única variável, não tem qualquer outra dependência em L para além da que é exibida explicitamente nesta equação, e ρ é o expoente
crítico de F (t). Se tentarmos medir esta função através de simulações de Monte Carlo
com sistemas de tamanhos diferentes, o resultado deverá ser independente do tamanho do
sistema, desde que usemos os valores correctos dos expoentes críticos e da temperatura
crítica, pelo que estes podem ser estimados a partir da tentativa de colapso das curvas
medidas para diferentes tamanhos do sistema. Estes colapsos, bem como os valores de Tc
e dos expoentes que lhes deram origem, são exibidos nas figuras 4.11.
ν = 1.2; Tc = 0.58
0.7
ν = 1.2; β = 0.125; Tc = 0.58
a)
c)
b)
24
32
40
48
56
64
24
32
40
48
56
64
1
L
χst
−γ/ν
L
β/ν
0.5
24
32
40
48
56
64
0.1
< mst >
0.6
Ust
ν = 1.2; γ = 2; Tc = 0.58
0.5
0.05
0.4
0.3
-4
-2
0
1/ν
L t
2
4
0
-4
-2
0
1/ν
L t
2
4
0
-4
-2
0
1/ν
2
4
L t
Figura 4.11: Colapsos das curvas de a) cumulante de Binder, b) magnetização e c) susceptibilidade para L = 24 a 64.
75
4 Redes regulares de dipolos pontuais
Os colapsos foram obtidos com os expoentes ν = 2, β = 0.125, γ = 2 e com a temperatura crítica Tc = 0.58. Estes valores estão razoavelmente de acordo com os apresentados
por Tomita [128], tendo em conta os limites desta forma de análise de tamanho finito e os
tamanhos considerados. Os expoentes encontrados nesse trabalho são ν = 1.19, β = 0.104
e Tc = 0.581, e não foi apresentado nenhum valor para γ. No entanto, o parâmetro de
ordem considerado é diferente - enquanto que o usado neste trabalho não distingue entre
o contínuo de estados do tipo micro-vórtice com diferentes ângulos θ, o parâmetro considerado por Tomita pretende descrever a selecção dos estados alinhados com os eixos da
rede. Para além isso, nesse trabalho é utilizada uma função de scaling “não convencional”
proposta por Campbell et al. [133], em alternativa à da eq. 4.13, que consiste em usar
como temperatura reduzida, em vez de t = (T − Tc )/Tc , o parâmetro t0 = (T − Tc )/T , e
a função de scaling
√ ρ/ν √ 1/ν G t0 L T
.
(4.14)
G(T, L) = L T
A motivação apresentada é que esta forma de scaling, derivada a partir da forma estrutural de expansões de alta temperatura, permitirá o colapso das curvas numa gama
mais alargada de temperaturas, não só próximo da região crítica, resultando nos mesmos
expoentes críticos com um erro inferior. Para permitir a comparação, efectuámos também
a análise com esta forma não convencional, e dentro da precisão deste procedimento não
encontrámos diferenças significativas.
4.2 Cadeia unidimensional
Como vimos no início deste capítulo, a uma dimensão o decaimento da interacção dipolar
é demasiado forte para que seja possível o estabelecimento de ordem de longo alcance a
temperatura finita [113]. As medições do parâmetro de ordem (neste caso, a magnetização
ao longo da direcção da cadeia) por simulação de Monte Carlo exibem um comportamento
semelhante ao parâmetro de ordem da rede quadrada (figura 4.12a), mas a análise do
cumulante de Binder revela a ausência de um ponto crítico, dado que a temperatura a
que as curvas se cruzam diminui sucessivamente com o tamanho do sistema (figura 4.12b),
sugerindo que a temperatura crítica tende para zero, como seria de esperar.
No entanto, a medição dos tempos de autocorrelação na energia e na magnetização
revela dificuldades de convergência que nem mesmo a combinação de selecção optimizada com o método de Parallel Tempering consegue ultrapassar (figura 4.13). Com estes
tempos de autocorrelação, seriam necessárias simulações bastante mais longas do que as
efectuadas para tornar possível uma análise semelhante à que fizemos para a rede quadrada.
76
4.2 Cadeia unidimensional
1
a)
b)
L=64
L=256
L=576
0.8
L=64
L=256
L=576
0.6
0.6
0.4
U
<m>
0.5
0.3
0.4
0.2
0.2
0.1
0
0.125
0.25
0.5
1
0
2
0.125
0.25
T
0.5
1
2
T
Figura 4.12: Magnetização e cumulante de Binder para uma cadeia unidimensional de
dipolos pontuais, usando a combinação de selecção optimizada com o método de Parallel
Tempering.
1000
1e+05
a)
b)
L=64
L=256
L=576
L=64
L=256
L=576
10000
100
< τE >
< τM >
1000
100
10
10
1
1
0.25
0.5
T
1
2
0.25
0.5
1
2
T
Figura 4.13: Tempos de autocorrelação na energia e na magnetização para uma cadeia
unidimensional de dipolos pontuais. Mesmo com a combinação de selecção optimizada
e Parallel Tempering, os tempos de autocorrelação tornam-se demasiado grandes para
medir com o tempo de simulação utilizado (107 passos de Monte Carlo).
Quanto à origem destes tempos de autocorrelação elevados, a taxa de aceitação, com as
optimizações utilizadas, continua a ser alta (da ordem de 85%), pelo que não será essa a
causa das dificuldades de convergência. O comportamento dos tempos de autocorrelação
sugere que a baixas temperaturas, os estados visitados não diferem muito na energia e
magnetização, podendo corresponder simplesmente ao movimento de paredes de domínio,
pelo que este comportamento deverá ser devido à utilização de uma dinâmica de single
spin flip.
77
4 Redes regulares de dipolos pontuais
4.3 Conclusões
Neste capítulo, aplicámos os métodos de simulação expostos no capítulo anterior a dois
sistemas ordenados de dipolos pontuais. Na rede quadrada, as dificuldades de convergência foram ultrapassadas com uma combinação de selecção optimizada com o método de
Parallel Tempering, com reduções de uma a duas ordens de grandeza nos tempos de autocorrelação. Estas optimizações tornaram possível a caracterização do comportamento
crítico do sistema, embora seja necessário simular sistemas maiores. Já na cadeia unidimensional, foi identificado um comportamento patológico a baixas temperaturas que
sugere a necessidade de utilização de algoritmos de cluster para permitir uma exploração
eficiente do espaço de fase.
78
5 Plano desordenado de dipolos
pontuais
Para além das redes ordenadas consideradas no capítulo anterior, é relevante considerar
também o caso de sistemas com desordem posicional, muito comuns entre os sistemas
magnéticos experimentais. São conhecidos relatos de “ferromagnetismo dipolar” (ordem
ferromagnética devida às interacções dipolares e não às interacções de troca) em diversas
condições (e.g. Zhang et al. [134], Kechrakos et al. [135], Takzei et al. [136], Denisov et
al. [137, 138, 139], Held et al. [140], Yamamoto et al. [141]). No entanto, esta ordem está
relacionada com realizações específicas da desordem, existência de anisotropia uniaxial e
elevados tempos de relaxação de estados meta-estáveis. Apesar de alguns destes trabalhos
incluírem análises de campo médio que parecem indicar a possibilidade de ordem de longo
alcance a temperatura finita, se fizermos simulações com médias sobre a desordem de um
sistema de dipolos pontuais uniformemente distribuídos e com acoplamento puramente
dipolar, essa ordem não se verifica.
Mesmo sem a existência de ordem de longo alcance, em sistemas desordenados usados como sensores graças ao efeito de magneto-resistência gigante (GMR) ou magnetoresistência de túnel (TMR), as correlações magnéticas devidas às interacções entre as
partículas têm um efeito substancial na resistência eléctrica da amostra em que estão
inseridas (e.g. Lopes et al. [142]). Neste capítulo, apresentamos resultados da influência
da concentração de partículas e da temperatura nessas correlações.
5.1 Modelo
Para medir as correlações entre momentos magnéticos num sistema desordenado em função do campo aplicado, foram geradas pequenas amostras circulares (figura 5.1) de partículas (esferas uniformemente magnetizadas), com uma partícula no centro e as restantes
partículas colocadas aleatoriamente à volta da primeira. As partículas não se podem
sobrepor, ou seja, para cada par de partículas i e j
(xi − xj )2 + (yi − yj )2 > D2
(5.1)
79
5 Plano desordenado de dipolos pontuais
onde D é o diâmetro das partículas. O raio da amostra é dado por
s
Ra =
N D2
4C
(5.2)
onde N é o número de partículas e C é a concentração planar de partículas, que foi variada
entre 10 e 30%.
Figura 5.1: Representação esquemática de cada amostra desordenada.
O facto de estarmos a lidar com um sistema desordenado obriga-nos a fazer médias
sobre a desordem, o que volta a ser bastante exigente do ponto de vista computacional.
A figura 5.2 mostra a média sobre 10000 amostras de 50 partículas cada, com uma concentração de 30%, da densidade de centros de partículas em função da posição relativa à
partícula central.
Figura 5.2: Exemplo da distribuição de partículas - média sobre 10000 amostras de 50
partículas cada, com uma concentração de 30%. Em todas as amostras existe uma partícula no centro, que não é exibida no gráfico, e o centro das outras partículas não pode
estar a menos de um diâmetro D de distância do centro.
80
5.2 Efeitos de tamanho finito
A energia do sistema é dada pela energia de interacção entre cada uma das partículas
de cada amostra mais o termo de Zeeman devido ao campo aplicado H, ou seja,
E=−
X χ
µk H χ −
k
1 X χλ χ λ
U µ µ ,
2 k,l kl k l
(5.3)
onde, mais uma vez, a matriz de interacção é dada por
Uklχλ = δ kl
3uχkl uλkl − δχλ
.
3
rkl
(5.4)
Dado que o sistema é desordenado, não vamos usar o método de réplicas periódicas
descrito no capítulo anterior, pelo que é necessário tratar de outra forma a questão do
tamanho finito da amostra.
5.2 Efeitos de tamanho finito
Como podemos ver na a distribuição espacial da magnetização (figura 5.3), o tamanho
finito da amostra provoca uma deformação da magnetização média no interior da amostra.
Junto à fronteira da amostra, a magnetização decresce significativamente.
Figura 5.3: Exemplo da distribuição da magnetização com condições de fronteira abertas.
Para ultrapassar este problema, foi gerado um plano com Np N partículas, que
inclui a amostra de simulação, com N partículas, no centro (figura 5.4). A diferença
entre as partículas que estão dentro e fora da amostra de simulação é que as partículas
que estão fora vão ter todas a mesma magnetização, calculada de forma auto-consistente.
Isto permite calcular uma matriz de interacção, Lχλ
k , que, para um determinado valor
da magnetização externa mλL , produz um campo local, em cada uma das partículas no
81
5 Plano desordenado de dipolos pontuais
interior da amostra, dado por
χ
λ
HL,k
= −Lχλ
k mL .
(5.5)
Este campo local é simplesmente adicionado ao campo externo, e a energia será então
dada por
X χ
1 X χλ χ λ
χ
(5.6)
U µ µ .
E=−
µk (H χ + HL,k
)−
2 k,l kl k l
k
Figura 5.4: Representação esquemática de cada amostra desordenada inserida num plano
de partículas uniformemente magnetizadas.
O número de partículas no exterior da amostra deve ser tal que o campo criado por
esta distribuição de partículas uniformemente magnetizadas, HL , deixe de depender desse
número. A figura 5.5 mostra a variação da componente x desde campo, ao longo do eixo
x e em função da razão Np /N , onde se verifica que é necessário chegar a um factor da
ordem de 100 para eliminar a dependência de HL com essa razão.
0.8
Np = 2 N
Np = 10 N
Np = 50 N
Np = 100 N
HL
0.6
0.4
0.2
0
-6
-4
-2
0
2
4
6
x/D
Figura 5.5: Variação do campo criado pela distribuição de partículas uniformemente
magnetizadas no exterior da amostra com a razão entre o número de partículas total e o
tamanho da amostra.
82
5.3 Simulações
A magnetização das partículas exteriores é definida de forma auto-consistente. Inicialmente, a magnetização das partículas que se encontram no exterior da amostra é colocada
a zero, e é calculada a magnetização média, por Monte Carlo, do interior da amostra. A
magnetização das partículas exteriores é então igualada à magnetização média do interior, e o procedimento é repetido até a diferença entre as magnetizações no interior e no
exterior ser inferior a 1%.
Figura 5.6: Exemplo da distribuição da magnetização paralela ao campo com magnetização constante e não nula no exterior da amostra.
Com este esquema, a magnetização da amostra pode aumentar, devido à magnetização
no exterior da amostra, e o perfil de magnetização no interior da amostra passa a ser mais
uniforme (figura 5.6), diminuindo dessa forma os efeitos de tamanho finito da amostra.
Para cada uma das amostras, foi calculada a aproximação de altas temperaturas, até
segunda ordem, da magnetização média de cada partícula e da correlação média de cada
partícula com a partícula central. Essas grandezas são também calculadas através de
simulações de Monte Carlo usando o algoritmo de Metropolis e a optimização de taxas
de aceitação referida anteriormente. No final, foram feitas médias destas grandezas sobre
todas as amostras em função da posição de cada partícula na amostra.
5.3 Simulações
A existência de um campo aplicado facilita a convergência das simulações, pelo que não
é necessário o uso do método de Parallel Tempering. No entanto, a taxa de aceitação a
campos altos pode ser demasiado baixa, se a proposta for uniforme, o que nos levou a
utilizar uma optimização semelhante à usada no capítulo anterior. A figura 5.7a mostra
o decaimento da taxa de aceitação, quando a proposta é uniforme, à medida que o campo
aumenta, devido ao custo energético de propostas que fariam os momentos desviar-se
83
5 Plano desordenado de dipolos pontuais
significativamente da direcção do campo.
1
1
a)
b)
S. Optimizada
S. Optimizada
(com campo externo e dipolar)
0.9
(com campo externo)
0.6
<A>
<A>
0.8
0.4
0.8
S. Optimizada
(com campo externo)
0.7
0.2
0.6
S. Uniforme
0
0.01
S. Uniforme
β=1
β=1
0.1
1
10
βH
100
0.5
0.01
0.1
1
βH
10
100
Figura 5.7: Exemplos de taxa de aceitação quando a proposta é uniforme, optimizada com
inclusão do campo externo e optimizada com inclusão do campo externo e da componente
planar do campo dipolar.
A introdução do campo externo na probabilidade de selecção, como descrito na secção 3.1.2, permite evitar este problema. Como podemos ver na figura 5.7a, esta optimização, incluindo o termo proporcional ao campo externo em
Ef − Ei = −H I,k · (mk,f − mk,i ) − Hext (cos θk,f − cos θk,i ).
(5.7)
na probabilidade de selecção,
Sf i ∝ eβHext cos θk,f e Sif ∝ eβHext cos θk,i
faz com que a taxa de aceitação tenda para 100% a campos altos, enquanto que anteriormente tendia para zero. Ainda assim, à medida que baixamos a temperatura, a taxa de
aceitação baixa também, devido ao termo dipolar que deixámos de fora.
Para melhorar a taxa de aceitação a campos fracos e baixas temperaturas, podemos
voltar a introduzir a componente planar dos campos locais na probabilidade de selecção,
k
k
Sf i ∝ eβHk cos θk,f e Sif ∝ eβHk cos θ̃k,i
e
84
5.4 Magnetização
k
onde θe é o ângulo que o momento faz com a componente planar do campo local H k =
k
H ext + H I,k , como descrito na secção 3.1.2. Como podemos ver na figura 5.7b, esta
alteração melhora significativamente a taxa de aceitação em toda a gama de parâmetros
em que estamos interessados. Os resultados que se seguem foram obtidos utilizando este
método.
5.4 Magnetização
Para partículas livres em que o momento magnético tem três componentes cartesianas
(apêndice A), a magnetização total na direcção do campo aplicado é descrita pela função
de Langevin,
1 2π
D
µk
E
0
u eβHu dφdu
1
= coth(βH) −
≡ L(βH),
= −11 0 2π
βHu dφdu
βH
e
−1 0
(5.8)
o que corresponde à aproximação de ordem zero da expansão de altas temperaturas. A
temperaturas elevadas espera-se que o comportamento da magnetização seja bem descrito
por esta função, para qualquer valor do campo aplicado, e que comecem a ser necessárias
correcções à medida que baixamos a temperatura, como se pode ver nas figuras 5.8 e5.9.
C = 30%
β = 0.1
β = 0.5
β=1
β = 1.5
1
Ordem 0
Ordem 1
Ordem 2
Sim.
Mx - L(βH)
0.8
0.6
0.4
0.2
0
0.01
0.1
1
βH
10
0.01
100
0.1
1
βH
10
0.01
100 0.1
1
βH
10
0.01
100
0.1
1
βH
10
100
Figura 5.8: Comparação da magnetização na direcção do campo aplicado, em função da
magnitude do mesmo. A concentração foi fixada em 30% e a temperatura diminui da
esquerda para a direita. O resultado da ordem zero da expansão de altas temperaturas é
exibido a preto, até primeira ordem a vermelho, até segunda ordem a verde e os resultados
da simulação a azul.
85
5 Plano desordenado de dipolos pontuais
C = 30%
β = 0.1
β = 0.5
β=1
0.14
Ordem 1
Ordem 2
Sim.
0.12
Mx - L(βH)
β = 1.5
0.1
0.08
0.06
0.04
0.02
0
0.01
0.1
1
βH
10
0.01
0.1
1
βH
10
0.01
0.1
1
βH
10
0.01
0.1
1
βH
10
100
Figura 5.9: Comparação da diferença entre a magnetização e a função de Langevin, em
função do campo aplicado. A concentração foi fixada em 30% e a temperatura diminui
da esquerda para a direita. O resultado da expansão de altas temperaturas até primeira
ordem é exibido a vermelho, até segunda ordem a verde e os resultados da simulação a
azul.
H = 1/β
C = 0.2
C = 0.1
C = 0.3
0.49
0.48
Ordem 0
Ordem 1
Ordem 2
Sim.
0.47
Mx
0.46
0.45
0.44
0.43
0.42
0
0.5
β
1
1.5 0
0.5
β
1
1.5 0
0.5
β
1
1.5
Figura 5.10: Magnetização paralela ao campo para diferentes temperaturas e H = 1/β, e
concentração C = 10% a 30%, com as diferentes ordens da expansão (até primeira ordem
a vermelho, até segunda ordem a verde) e os resultados da simulação (a azul).
86
5.5 Correlações magnéticas
Repare-se que a aproximação de altas temperaturas não se aproxima do resultado da
simulação de forma monótona. Comparando a aproximação de ordem zero (linha preta)
com a de primeira ordem (linha vermelha) e o resultado da simulação (linha azul), verificamos que correcção de primeira ordem estima a magnetização por excesso, o que depois
é corrigido pelo termo de segunda ordem (linha verde).
A concentração afecta também a magnitude das correcções, ao tornar as interacções
mais relevantes. A figura 5.10 mostra a magnetização paralela ao campo em função da
temperatura, para H = 1/β, e concentrações C = 10% a 30%.
À medida que baixamos a temperatura (aumentamos β) ou aumentamos a concentração
C, as interacções tornam-se mais relevantes e aumentam a magnetização ao longo do
campo. Note-se, mais uma vez, que a expansão estima esta correcção por excesso, e que
o desvio aumenta com o aumento da concentração e a diminuição da temperatura.
5.5 Correlações magnéticas
Para partículas livres, a correlação entre os momentos magnéticos de partículas diferentes,
definida como
D
E
D E
Cijαβ = µαi µβj − hµαi i µβj
(5.9)
é sempre nula, como vimos no capítulo 2, pelo que só iremos apresentar as correcções de
primeira e segunda ordem, bem como o resultado das simulações. Para cada partícula,
foi medida a correlação do seu momento magnético com o da partícula central, e esta
correlação foi mapeada em função da posição relativa ao centro (figura 5.11). A correlação
entre as diferentes componentes dos momentos magnéticos depende essencialmente da
posição relativa das partículas, como se vê na figura.
Figura 5.11: Exemplo da distribuição da correlação com a partícula central das componentes paralela e perpendicular da magnetização. Note-se a importância da relação entre
as componentes do momento magnético e a posição relativa das partículas.
87
5 Plano desordenado de dipolos pontuais
C = 30%
β = 0.1
β = 0.5
β=1
β = 1.5
0.005
Ordem 1
Ordem 2
Sim.
0.004
Ci0
xx
0.003
0.002
0.001
0
0.01
0.1
1
βH
10
0.01
0.1
1
βH
10
0.01
0.1
1
βH
10
0.01
0.1
1
βH
10
100
Figura 5.12: Comparação da correlação na direcção do campo, em função do campo, para
diferentes temperaturas e concentração C = 30%, com as diferentes ordens da expansão
(primeira ordem a vermelho, até segunda ordem a verde) e os resultados da simulação (a
azul).
H = 1/β
C = 0.1
C = 0.2
C = 0.3
0.002
Ordem 1
Ordem 2
Sim.
xx
Ci0
0.0015
0.001
0.0005
0
0
0.5
β
1
1.5 0
0.5
β
1
1.5 0
0.5
1
β
1.5
2
Figura 5.13: Magnetização paralela ao campo para diferentes temperaturas e H = 1/β, e
concentração C = 10% a 30%, com as diferentes ordens da expansão (primeira ordem a
vermelho, até segunda ordem a verde) e os resultados da simulação (a azul).
88
5.6 Conclusões
A figura 5.12 mostra a correlação na direcção paralela ao campo, integrada sobre toda a
amostra, à semelhança do que fizemos com a magnetização. Esta análise permite comparar
os valores numéricos obtidos com as diferentes aproximações. A correlação integrada
também aumenta com β e com C, como podemos ver, para H = 1/β, na figura 5.13.
No entanto, verificamos que para temperaturas baixas e/ou concentrações altas, o termo
de segunda ordem da expansão já não é suficiente para descrever o comportamento do
sistema.
5.6 Conclusões
Foi desenvolvida uma aplicação de simulações de Monte Carlo de um plano desordenado
de dipolos pontuais com médias sobre a desordem e com o cálculo, para cada amostra,
da expansão de altas temperaturas até segunda ordem da magnetização e da correlação
média com a partícula central. Isto permite verificar a validade das diferentes ordens
da expansão de altas temperaturas, em função da temperatura e da concentração. Em
particular, verificamos que a aproximação de primeira ordem começa a descrever o efeito
das interacções, quando comparada com a aproximação de ordem zero, mas estima por
excesso a magnetização e a correlação magnética, quando comparadas com os resultados
das simulações. O termo de segunda ordem tipicamente melhora a aproximação, desde
que a concentração não seja demasiado alta (inferior a 30%). Estes resultados para as
propriedades magnéticas estão a ser introduzidos em modelos de transporte para obter as
consequências em termos de magneto-resistência de sistemas constituídos por partículas
magnéticas em substratos metálicos e isoladores. Esta tarefa permitirá aplicações práticas
dos resultados apresentados.
89
5 Plano desordenado de dipolos pontuais
90
6 Partículas de tamanho finito
Nos capítulos anteriores, utilizámos a aproximação de dipolos pontuais para descrever as
partículas que constituem os sistemas nano-estruturados e as interacções entre elas. Para
ter em conta efeitos de tamanho e forma das partículas, é necessário, em primeiro lugar,
saber calcular, no regime clássico, a energia magnetostática de uma partícula magnética
(e.g. Aharoni [143]). Esta energia não é mais do que a interacção da sua magnetização
M (r) com o campo magnético H(r) criado por essa magnetização, ou seja
µ0
Em = −
2
M (r) · H(r)dr
(6.1)
V
(onde µ0 = 4π × 10−7 N/A2 é a permeabilidade magnética do vácuo e o integral é feito
sobre o volume da partícula). Para avaliar a energia de um sistema de partículas, é
normalmente usado o teorema de reciprocidade (e.g. Jackson [28]), que diz que a energia
de interacção de uma distribuição de magnetização M 1 (r) com o campo magnético H 2 (r)
criado por uma distribuição de magnetização M 2 (r) é igual à energia de interacção desta
com o campo magnético H 1 (r), ou seja
M 1 (r) · H 2 (r)dr =
V
M 2 (r) · H 1 (r)dr.
(6.2)
V
Figura 6.1: Representação esquemática do teorema da reciprocidade.
Para calcular esta energia de interacção, é necessário calcular o campo magnético produzido por uma determinada distribuição de magnetização, como veremos de seguida.
91
6 Partículas de tamanho finito
6.1 Função de forma
Genericamente, o cálculo de energias de energias magnetostáticas directamente a partir
da definição não é tão simples como possa parecer. No que se segue, iremos usar um
formalismo no espaço de Fourier, referido por Aharoni [143] e mais recentemente sistematizado por Beleggia et al. [144, 145, 146, 147], que simplifica em grande medida esse
cálculo. Para tal, vamos começar por definir uma função de forma
D(r) =


1
dentro da partícula

0
fora da partícula
(6.3)
para extender os integrais sobre o volume a integrais sobre todo o espaço
f (r)dr =
D(r)f (r)dr.
(6.4)
V
0
1
Figura 6.2: Representação esquemática da função de forma.
Isto permite transferir a questão da forma das partículas para esta função de forma,
nomeadamente para o cálculo da transformada de Fourier da função de forma, designada
por amplitude de forma,
D(k) =
D(r)e
−ik·r
dr =
e−ik·r dr.
(6.5)
V
Este último integral é o único ponto neste formalismo onde a informação sobre a forma
é usada nos limites de integração (no apêndice B são calculadas as amplitudes de forma
para os casos particulares de formas esféricas e cilíndricas). Se M (r) incluir o factor de
92
6.1 Função de forma
forma, e.g.
M (r) = M (r)D(r)m̂(r),
(6.6)
onde M (r) é a magnitude da magnetização e m̂(r) a sua direcção, a energia magnetostática fica então
µ0
M (r) · H(r)dr,
(6.7)
Em = −
2
com o integral extendido a todo o espaço em vez de confinado ao volume da partícula.
Figura 6.3: Representação esquemática de um corpo com magnetização não uniforme.
Para calcular H(r), começamos por usar o conhecido facto (e.g. [28]) de que o potencial
vector magnético devido à magnetização da partícula é dado por
µ0
A(r) =
4π
M (r 0 ) ×
r − r0
dr 0 ,
|r − r 0 |3
(6.8)
e pode ser escrito, com o auxílio do teorema da convolução e a linearidade do produto vectorial (e representando por F o operador transformada de Fourier, F[f (r)] =
f (r)e−ik·r dr), como
A(k) = F [A(r)]
r
µ0
=
F [M (r)] × F 3
4π
r
iµ0
= − 2 M (k) × k,
k
(6.9)
o que reduz o cálculo do potencial vector a um produto vectorial, se for possível calcular
a transformada de Fourier da magnetização. O potencial vector permite-nos calcular a
93
6 Partículas de tamanho finito
indução magnética B(r) = ∇ × A(r) através da sua transformada de Fourier
B(k) = −ik × A(k)
µ0
= 2 k × M (k) × k
k
(6.10)
o que, usando a identidade vectorial k ×M (k)×k = M (k)k 2 −k [k · M (k)], nos permite
escrever
)
(
k [k · M (k)]
.
(6.11)
B(k) = µ0 M (k) −
k2
Relacionando o campo magnético H com a indução B através da relação comum B =
µ0 (M + H), identificamos a grandeza
H(k) = −
k [k · M (k)]
k2
(6.12)
com a transformada de Fourier do campo desmagnetizante H(r), que pode então ser
escrito, usando a transformada inversa
F
−1
1
[f (k)] = 3
8π
f (k)eik·r dk = f (r),
como
1
H(r) =
8π 3
(6.13)
H(k)eik·r dk.
Escrevendo a equação 6.7 com recurso a esta transformada inversa, obtemos a expressão
para a energia magnetostática em termos da transformada de Fourier acima definida
Em
µ0
= −
2
1
M (r) ·
8π 3
(
)
H(k)e
ik·r
dk dr,
(6.14)
o que, substituindo H(k) pela equação 6.12 e reconhecendo no integral em dr a transformada de Fourier M (−k), pode ser simplificado para
µ0
Em =
16π 3
[k · M (−k)] [k · M (k)]
dk.
k2
(6.15)
Esta é a expressão mais geral para a energia magnetostática de um sistema magnético.
Na secção 6.2, vamos introduzir a simplificação de considerar a magnetização uniforme
em todas as regiões do espaço onde está definida.
94
6.2 Magnetização uniforme
6.2 Magnetização uniforme
Para o caso mais simples de um corpo uniformemente magnetizado (figura 6.4), a transformada de Fourier da magnetização torna-se simplesmente
M (k) = M D(k)m̂,
(6.16)
onde m̂ é o versor unitário na direcção da magnetização, e a energia passa a ser dada por
µ0 M 2
Em =
16π 3
|D(k)|2
(m̂ · k)2 dk.
2
k
(6.17)
Figura 6.4: Representação esquemática de um corpo com magnetização uniforme.
Para exemplificar o uso deste formalismo, são considerados no Apêndice B algumas
formas particulares como esferas (Osborn [148]) e cilindros (Joseph e Schlomann [149])
uniformemente magnetizados, neste último caso em duas configurações distintas.
6.2.1 Tensor desmagnetizante local
No caso de corpos elipsoidais uniformemente magnetizados, o campo desmagnetizante é
também uniforme, e a energia magnetostática pode ser expressa através de um tensor
desmagnetizante global. Se a forma for diferente, isto já não se verifica, mas podemos definir, por analogia, um tensor desmagnetizante local N αβ (r), que relaciona a magnetização
M (r) com o campo H(r) produzido pela mesma,
H α (r) = −
X
N αβ (r)M β (r).
(6.18)
β
95
6 Partículas de tamanho finito
Comparando com a expressão 6.12, podemos identificar o termo
D(k)
kαkβ
k2
(6.19)
com a transformada de Fourier do tensor desmagnetizante,
kαkβ
D(k) 2
k
"
N
αβ
(r) = F
−1
#
(6.20)
que, usando o teorema da convolução, pode ser escrito como
"
N
αβ
(r) = D(r) ⊗ F
onde a função
−1
kαkβ
,
k2
#
(6.21)
f (r) ⊗ g(r) ≡
f (r 0 )g(r 0 − r)dr 0
(6.22)
é a convolução das funções f e g. O segundo termo desta convolução,
"
F
−1
kαkβ
k2
#
(6.23)
é o tensor de interacção entre dipolos pontuais, que iremos designar aqui simplesmente
por tensor dipolar, Dαβ (r). Isto permite olhar para o tensor desmagnetizante como a
convolução entre a função de forma e o tensor dipolar
N αβ (r) = D(r) ⊗ Dαβ (r),
(6.24)
isto é, o tensor desmagnetizante local é obtido através do efeito do tensor dipolar a partir
de todos os pontos da partícula. No caso da partícula ser de facto um dipolo pontual, a
função de forma é uma função delta e o tensor desmagnetizante é igual ao tensor dipolar.
É importante notar que, para um sistema magnético uniformemente magnetizado, toda a
informação está contida nesta expressão. Em todo o caso, no final deste capítulo vamos
abordar também o caso de magnetização não uniforme, para o qual isto deixa de se
verificar.
6.2.2 Interacção entre partículas
Até agora, temos tratado apenas o caso de uma única partícula, mas não há nada que nos
impeça de considerar que a magnetização está distribuída por diferentes regiões do espaço
96
6.2 Magnetização uniforme
(figura 6.5). Vamos considerar então duas partículas magnetizadas uniformemente,
M (r) = Mi Di (r − r i )m̂i + Mj Dj (r − r j )m̂j .
(6.25)
Figura 6.5: Representação esquemática da magnetização distribuída por diferentes regiões
do espaço.
A transformada de Fourier desta magnetização é dada por
M (k) = Mi Di (k)e−ik·ri m̂i + Mj Dj (k)e−ik·rj m̂j ,
(6.26)
que pode ser substituída na equação 6.15 resultando nos dois termos das energias magnetostáticas de cada partícula mais um termo de interacção
Mi Mj µ0
EI =
8π 3
Di (k)Dj (−k)eik·rij
(m̂i · k) (m̂j · k)
dk,
k2
(6.27)
onde r ij = r j − r i é a posição relativa das duas partículas e onde contribui apenas a parte
real do integrando. Este integral corresponde a uma transformada de Fourier inversa,
relativamente a r ij , pelo que a energia de interacção pode ser escrita como
"
EI =
µ0 mαi mβj Fr−1
ij
kαkβ
Di (k)Dj (−k) 2
k
#
(6.28)
e, usando mais uma vez o teorema da convolução,
"
EI =
µ0 mαi mβj Fr−1
ij
[Di (k)Dj (−k)] ⊗ Fr−1
ij
kαkβ
.
k2
#
(6.29)
O segundo termo desta convolução é, como já vimos, o tensor dipolar Dαβ (r ij ), enquanto
que o primeiro termo corresponde a uma correlação cruzada
C(r) = Di (r) ⊗ Dj (−r),
(6.30)
97
6 Partículas de tamanho finito
ou seja
EI = µ0 mα1 mβ2 C(r ij ) ⊗ Dαβ (r ij ).
(6.31)
O aspecto mais interessante desta expressão é que fornece uma forma geral de calcular
o tensor de interacção
Uijαβ = C(r ij ) ⊗ Dαβ (r ij )
entre duas partículas de forma arbitrária. Dado que os métodos analíticos e numéricos
apresentados no capítulo anterior não dependem da forma particular do tensor de interacção, a introdução deste tensor, uma vez calculado, é a única alteração necessária para
tratar sistemas constituídos por partículas de tamanho finito uniformemente magnetizadas.
6.2.3 Sistemas periódicos
No caso de um sistema periódico de elementos idênticos (com a mesma forma, tamanho
e distribuição de magnetização),
M (r) =
X
m(r − r i ),
(6.32)
i
onde r i é a posição de cada ponto da rede e m(r −r i ) é a magnetização de cada partícula,
incluindo o seu factor de forma D(r − r i ). A transformada de Fourier desta magnetização
é dada por
M (k) =
=
M (r)e−ik·r dr
X
(6.33)
m(r − r i )e−ik·r dr.
i
Dada a periodicidade de m(r), fazemos a mudança de variável r 0 = r − r i para escrever
M (k) =
X
−ik·r 0
m(r )e
0
i
=
X
i
98
m(k)e−ik·ri ,
!
dr
0
e−ik·ri
(6.34)
6.3 Rede de discos
onde m(k) = m(r)e−ik·r dr é a transformada de Fourier da magnetização de uma única
partícula. A energia de interacção por partícula é então dada por
Em =
µ0
16π 3
h
[k · m(−k)] k ·
−ik·r i
i m(k)e
P
k2
i
dk,
(6.35)
onde o somatório sobre a rede só não será nulo se k = g, sendo g um vector da rede
recíproca, resultado este que iremos usar de seguida para o caso particular de uma rede
de discos (cilindros em configuração planar).
6.3 Rede de discos
Recentemente, Kakazei et al. [150] encontraram experimentalmente uma anisotropia de
quarta ordem no espectro de ressonância ferromagnética (FMR) de um sistema constituído
por uma rede quadrada de discos (figura 6.6).
Figura 6.6: Representação esquemática de uma rede quadrada, com período a, de discos
com diâmetro D e espessura t.
Na experiência, a amostra é sujeita a um campo magnético suficientemente forte para
manter a magnetização próxima a saturação, no plano da rede e das faces dos discos, e
a anisotropia é encontrada à medida que é variada a direcção do campo aplicado relativamente aos eixos da rede. O efeito é mais forte quando os discos estão próximos do
contacto entre si e desaparece à medida que se aumenta o período da rede relativamente
ao diâmetro dos discos, pelo que deverá ser devido às interacções entre discos. Esta manifestação explícita dos efeitos da interacção magnetostática é interessante tanto do ponto
de vista fundamental como das aplicações, pelo que a iremos estudar com recurso aos
métodos anteriormente expostos neste capítulo.
99
6 Partículas de tamanho finito
6.3.1 Resultados experimentais
A amostra experimental é constituída por discos de Permalloy com D = 1 µm de diâmetro
e t = 50 nm de espessura, fabricados por litografia de feixe de electrões e técnicas de liftoff, dispostos em redes quadradas com espaçamento a entre 1.1 e 2.5 µm. Estas dimensões
foram confirmadas por imagens microscópicas de microscopia de força atómica e de efeito
de túnel.
Anisotropia de segunda ordem
Anisotropia de quarta ordem
95
95
Hr (kA/m)
a=1.1µm
87,17 + 0.35 cos(2φH)
87,17 + 0.35 cos(2φH) + 1.94 cos(4φH)
90
90
85
85
a=2.5µm
89,88 + 0.38 cos(2φH)
0
π/2
π
φH
3π/2
0
π/2
π
φH
3π/2
2π
Figura 6.7: Dependência experimental do campo de ressonância Hr em função do ângulo φH , para parâmetro de rede a = 2.5 µm (a preto, à esquerda), ajustada por uma
anisotropia uniaxial e a = 1.1 µm (a vermelho, à direita), ajustada por uma sobreposição de anisotropia de quarta ordem e a anisotropia uniaxial, exibida separadamente em
pontilhado.
A experiência de ressonância ferromagnética (FMR) foi realizada à temperatura ambiente aplicando um campo magnético constante H, bastante forte, no plano da amostra
e fazendo um ângulo φH com um eixo da rede, e um (fraco) campo alternado de alta
frequência (9.8 GHz). O campo de ressonância Hr é a magnitude do campo H que maximiza a perda de energia na cavidade, ou seja, a que faz com que a frequência própria de
precessão dos momentos da amostra seja igual à frequência do campo alternado. O facto
100
6.3 Rede de discos
desta frequência ser bastante elevada faz com que o campo de ressonância seja superior
ao campo de saturação dos discos, evitando assim a formação de domínios ou vórtices.
91
4
a)
b)
90
89
H4 (kA/m)
Hr,av (kA/m)
3
88
2
1
87
86
0
0.2
0.4
(D/a)
0.6
0.8
1
0
0
0.1
0.2
2
0.3
0.4
(D/a)
0.5
0.6
4
Figura 6.8: Dependência experimental de a) Hr,av e b) H4 com o parâmetro de rede a.
Hr,av é inversamente proporcional a a2 e H4 é inversamente proporcional a a4 .
As dependências do campo de ressonância Hr com o ângulo φH para os casos de maior
e menor parâmetro da rede são apresentadas na figura 6.7. No caso do maior parâmetro
de rede, a = 2.5 µm, verifica-se a presença de uma anisotropia uniaxial descrita por
Hr = Hr,av + H2 cos 2φH .
(6.36)
A pequena anisotropia uniaxial H2 é aproximadamente constante em todas as amostras,
e o alinhamento com um dos eixos da rede sugere que seja devida a um pequeno desvio
da simetria da rede em relação a uma rede quadrada.
Para a amostra com parâmetro de rede a = 1.1 µm verificamos uma diminuição do valor
médio Hr,av e, para além da anisotropia uniaxial já referida (linha pontilhada vermelha
na figura 6.7), encontramos uma anisotropia de quarta ordem que se reflecte nos máximos
de Hr na direcção dos eixos da rede (geometria [10]) e mínimos nas direcções diagonais
(geometria [11]). Neste caso, Hr (φH ) é bem descrito por
Hr = Hr,av + H2 cos 2φH + H4 cos 4φH .
(6.37)
A variação de Hr,av e H4 com o parâmetro de rede a é exibida na figura 6.8.
O valor médio do campo de ressonância, Hr,av é inversamente proporcional a a2 , en-
101
6 Partículas de tamanho finito
quanto que a anisotropia de quarta ordem H4 é inversamente proporcional a a4 e anula-se
no mesmo limite.
6.3.2 Modelo teórico
Se considerarmos então esta rede quadrada de lado a no plano xy, formada por discos
de raio R e espessura t, e que a magnetização M , para além de não ter componente
segundo a direcção z, também não depende desta coordenada [151], podemos escrever a
transformada de Fourier da magnetização de um disco como
m(k) =
0
=
m(r 0 )e−ik·r dr0
(6.38)
m(ρ)e−iκ·ρ e−iqz dρdz,
onde, mais uma vez, mudámos para coordenadas cilíndricas tanto no espaço real, r =
ρ + z = (ρ cos θ, ρ sin θ, z), como no espaço de Fourier, k = κ + q = (κ cos φ, κ sin φ, q).
A integração em z pode ser feita como no apêndice B, de onde resulta
m(k) =
2 sin (qt/2)
m(κ).
q
(6.39)
Substituindo na equação 6.35, temos a energia magnetostática por partícula e por unidade
de espessura t
µ0
Em =
16π 3 t
2 sin (qt/2)
q
!2
[κ · m(−κ)] κ · i m(κ)e−iκ·ri
dκdq
κ2 + q 2
P
(6.40)
e usando o resultado da integração em q (eq. B.11), chegamos a
Em
µ0
=
8π 2
[κ · m(−κ)] κ ·
f (κt)
κ2
P
i
m(κ)e−iκ·ri
dκ,
(6.41)
onde f (x) = (e−x − 1 + x) /x. O somatório sobre os pontos da rede ri só não é nulo se κ
for um vector da rede recíproca, κ = g,
X
i
e−iκ·ri f n(κ)dκ =
2π
a
2 X
f n(g)
g
pelo que a energia pode ser escrita como um somatório sobre a rede recíproca,
102
(6.42)
6.3 Rede de discos
Em =
[g · m(−g)] [g · m(g)]
µ0 1 X
f
(gt)
.
2 a2 g
g2
(6.43)
De seguida, vamos começar por assumir uma magnetização uniforme na direcção do
campo.
Magnetização uniforme
Se considerarmos que a magnetização é uniforme e paralela ao campo aplicado, que vamos
manter ao longo do eixo x, então
R
2π
m (g) = M
ei g r cos θ rdrdθ = 2πRM
x
0
0
J1 (gR)
g
(6.44)
onde Jα (x) é a função cilíndrica de Bessel, e my,z (g) = 0. A energia magnetostática é
então dada por
Em
2π 2 R2 M 2 µ0 X
gx
=
f
(gt)
a2
g
g
!2
J1 (gR)
g
!2
.
(6.45)
A variação com o ângulo φH entre a direcção do campo e os eixos da rede pode ser
obtida numericamente mantendo fixa a direcção do campo e rodando a rede recíproca, ou
seja,
2π
(n1 cos φH − n2 sin φH , n1 sin φH + n2 cos φH ) ,
(6.46)
g=
a
verificando-se que a energia de uma rede quadrada de discos uniformemente magnetizados
no plano da rede não depende da direcção da magnetização, pelo que não é suficiente
para explicar os resultados experimentais. No entanto, a curva de absorção medida na
experiência está relacionada com o campo magnético total sentido por cada momento dos
discos. Este campo é constituído não só pelos campos aplicados mas também pelo campo
desmagnetizante devido à forma da partícula e devido à presença de outras partículas. O
primeiro é independente da direcção do campo em relação aos eixos da rede e da distância
entre discos, pelo que será necessário considerar a contribuição para o campo local devida
à presença dos outros discos da rede.
103
6 Partículas de tamanho finito
Campo magnetostático
O campo magnetostático no interior de um dos discos da rede pode ser obtido combinando
as equações 6.13 e 6.34,
1
H(r) = − 3
8π
h
k k·
−ik·r i
i m(k)e
P
k2
i
eik·r dk
(6.47)
e integrando em coordenadas cilíndricas no espaço de Fourier. As componentes do campo
local, integrando em coordenadas cilíndricas e fazendo a média sobre a espessura do disco
(apêndice B), serão dadas por, para α, β = x, y,
gαgβ β
1 X
f
(gt)
m (g) cos(g · ρ),
a2 β,g
g2
(6.48)
1 X
[1 − f (gt)] mz (g) cos(g · ρ).
a2 g
(6.49)
H α (ρ) = −
H z (ρ) = −
Supondo uma magnetização próxima da saturação na direcção do campo aplicado,
admitindo na equação 6.44 as amplitudes, infinitesimais e uniformes, de precessão µy e
µz , ou seja, m(ρ) = (M, µy , µz ), e calculando as transformadas de Fourier de cada uma das
componentes da magnetização, podemos identificar os factores de proporcionalidade entre
as componentes da magnetização e do campo com os elementos do tensor desmagnetizante
(secção 6.2.1),
N
αβ
g α g β J1 (gR)
2πR X
(ρ) = 2
f (gt) 2
cos(g · ρ)
a
g
g
g
(6.50)
2πR X
J1 (gR)
[1 − f (gt)]
cos(g · ρ),
2
a
g
g
(6.51)
N zz (ρ) =
sendo nulos os elementos N zx = N xz = N zy = N yz .
Segundo a teoria clássica de ressonância ferromagnética proposta por Kittel [27]1 , a
condição de ressonância é dada por
s
Hr (ρ) =
1
H02 + M 2
[N zz (ρ) − N yy (ρ)]2
N zz (ρ) + N yy (ρ) − 2N xx (ρ)
−M
4
2
(6.52)
A teoria de Kittel assume que o tensor desmagnetizante é diagonal e que a variação da componente
da magnetização paralela ao campo interno é desprezável. Neste caso, as equações de movimento
têm solução analítica, resultando numa divergência na susceptibilidade relativa ao campo alternado
quando a frequência própria de precessão, que depende tanto do campo aplicado como do campo
desmagnetizante, é igual à frequência do campo alternado.
104
6.3 Rede de discos
onde H0 = ω/γ, sendo ω a frequência do campo alternado aplicado e γ é o factor giromagnético do material. O campo de ressonância medido experimentalmente corresponde
ao máximo da absorção em toda a amostra, que iremos aproximar por
δ (Hext − Hr (ρ)) dρ,
I(Hext ) ∝
(6.53)
ρ<R
onde Hext é o campo aplicado. As curvas de absorção calculadas numericamente, usando
as equações 6.53, 6.52, 6.51 e 6.50, são apresentadas na figura 6.9.
a = 1.1 µm
a)
a = 2.5 µm
b)
[10]
[11]
u. a.
u. a.
[10]
[11]
60
80
H (kA/m)
100
120
60
80
100
120
H (kA/m)
Figura 6.9: Intensidade de absorção em função do campo aplicado Hext para a) espaçamento mínimo (a = 1.1 µm) e b) espaçamento máximo (a = 2.5 µm) e geometrias [10]
(φH = 0, a preto) e [11] (φH = π/4, a vermelho). O integral (eq. 6.53) foi aproximado pela
soma sobre 500x500 elementos do disco e a soma na rede recíproca até g = 500 × (2π/a).
Estes resultados, obtidos usando a aproximação de discos uniformemente magnetizados,
exibem já uma anisotropia em relação à direcção do campo aplicado. Na geometria [10]
a curva é mais larga e deslocada para valores mais elevados do campo aplicado do que
na geometria [11]. No entanto, a forma das curvas exibe comportamentos peculiares que
não são observados na experiência e a magnitude da anisotropia é mais pequena do que
a observada experimentalmente.
Estes aspectos devem ser devidos à aproximação de considerar magnetização uniforme
em todo o disco e apenas os termos diagonais do tensor desmagnetizante. O termo N yx
é responsável por um campo H y = −N yx M x , que é perpendicular ao campo aplicado
mas proporcional à magnetização na direcção do campo, o que induz uma deformação na
distribuição da magnetização dentro do disco. De seguida é apresentado um procedimento
variacional que permite estimar essa deformação.
105
6 Partículas de tamanho finito
Magnetização não uniforme
Mantendo a condição de magnetização planar e com magnitude constante em todo o disco,
m(ρ) = M (cos θ(ρ), sin θ(ρ), 0),
(6.54)
a deformação pode ser completamente descrita pelo ângulo θ(ρ), tal que a sua variação
arbitrária e infinitesimal produz uma correspondente variação de magnetização
δm(ρ) = ẑ × m(ρ)δθ(ρ)
(6.55)
onde ẑ é o vector unitário perpendicular ao plano, ou pela sua transformada de Fourier,
na forma da convolução
δm(g) = ẑ ×
X
m(g − g 0 )δθ(g 0 ).
(6.56)
g0
A condição de equilíbrio para esta deformação vai resultar da competição entre a energia
de Zeeman, devido à presença de um campo externo, e a energia magnetostática dentro
de cada partícula e entre partículas, ou seja, a energia total é dada por
Em = −µ0 H ext · m0 +
[g · m(g)] [g · m(−g)]
µ0 X
f (gt)
2
2a g
g2
(6.57)
onde m0 = m(g = 0) é a magnetização total do sistema e H ext é o campo aplicado.
Usando um tratamento variacional, esta energia pode ser minimizada substituindo m(g)
por m(g) + δm(g) e escrevendo os coeficientes de cada δθ(g 0 ) na condição δE = 0, o que,
desprezando os termos de segunda ordem em 1/Hext , leva ao sistema de equações
h
[g · m(g)] g · ẑ ×
1 X
H ext · (ẑ × m(−g 0 )) = 2
f (gt)
a g
g2
P
g0
i
m(g − g 0 )
(6.58)
ou, usando as propriedades do produto triplo,
(H ext × m(−g 0 )) · ẑ =
[g · m(g)] [(g × m(g − g 0 )) · ẑ]
1 X
f
(gt)
.
a2 g
g2
(6.59)
Este cálculo é simplificado se supusermos que o campo está sempre aplicado na direcção
do eixo x, e rodando a rede para estudar a anisotropia em função do ângulo que o campo
faz com os eixos da rede. Com o campo aplicado na direcção x, é esta também a maior
componente da magnetização, e a deformação é descrita pela componente y
106
6.3 Rede de discos
my (g) =
1
X
Hext a2
g 0 6=0
f (g 0 t)
[m(g 0 ) · g 0 ] [(m(g − g 0 ) × g0 ) · ẑ]
g 02
(6.60)
(onde trocámos os índices g e g 0 ). Esta expressão pode ser resolvida de forma iterativa,
calculando a iteração n a partir da iteração n − 1,
m(n)
y (g) =
h
1
Hext
X
a2
f (g 0 t)
m(n−1) (g 0 ) · g 0
i h
m(n−1) (g − g 0 ) × g0 · ẑ
i
g 02
g0
.
(6.61)
Começamos por notar que, para campos suficientemente fortes para saturar a magnetização, esta seria dada, como já vimos, por
m(0)
x (g) = 2πM R
J1 (gR)
e m(0)
y (g) = 0,
g
(6.62)
onde o índice zero corresponde à condição inicial da iteração auto-consistente. Sendo
assim, a primeira iteração da equação resulta
m(1)
y (g)
0 0
0
0
4π 2 M 2 R2 X
0 gx gy J1 (g R) J1 ((g − g )R)
=−
f
(g
t)
Hext a2 g0
g 02
g0
g − g0
(6.63)
onde já podemos encontrar uma fonte de anisotropia no termo gx0 gy0 , que não é invariante face a uma rotação da rede. A magnetização no espaço real é obtida através da
transformada de Fourier inversa, resultando em
m(1)
y (ρ) = −Θ (R − ρ)
gx gy
2πM 2 R2 X
f
(gt)
J1 (gR) cos (g · ρ)
Hext a2 g
g3
(6.64)
e a componente x pode ser calculada através da condição de normalização mx (ρ) =
q
M 2 − m2y (ρ). Somando sobre os vectores da rede recíproca
g = (2π/a)(n1 cos φH − n2 sin φH , n1 sin φH + n2 cos φH )
(6.65)
para φH = 0 e φH = π/4, é possível verificar que a deformação é substancialmente menor
no caso da orientação φH = π/4 ([11]) do que no caso φH = 0 [10] (figura 6.10), o que
está qualitativamente de acordo com os resultados experimentais, nomeadamente o facto
de a largura da linha de absorção ser maior no caso [10] do que no caso [11].
Para além disso, no caso [11], a deformação é convexa perto do centro do disco, e côncava
perto da periferia, enquanto que no caso [10] é sempre côncava, como no caso de um disco
isolado. Este efeito deverá ser devido a uma competição entre o campo desmagnetizante
107
6 Partículas de tamanho finito
0.1
[10]
[11]
0.2
0
0.1
my/M
my/M
-0.1
0
-0.2
-0.1
r = 0.3
r = 0.45
-0.3
-0.4
0
-0.2
0.1
0.2
r
0.3
0.4
0.5
0
π/4
π/2
θ
3π/4
π
Figura 6.10: Variação radial (a), ao longo da direcção de deformação máxima, e angular
(b), para r = 0.3 (linha sólida) e r = 0.45 (linha tracejada), da componente y da magnetização (eq. 6.64), para as geometrias [10] (a preto) e [11] (a vermelho), com espaçamento
entre discos a = 1.1 µm. Note-se que a magnitude da deformação é sempre maior no caso
[10] do que [11], tendo sinais opostos na zona central e coincidentes junto à periferia do
disco.
devido à forma da partícula, que produz sempre uma deformação dirigida para o interior
do disco, e o campo devido aos outros discos da rede, que depende da orientação do campo
externo relativamente aos eixos da rede. Este comportamento foi também verificado nas
simulações micro-magnéticas apresentadas por Kakazei et al. [150], e não pode ser descrito
por modelos de um único parâmetro como o que foi apresentado por Metlov [152].
6.4 Conclusões
Neste capítulo, foi abandonada a aproximação de dipolos pontuais e apresentado um
procedimento geral de introdução da forma das partículas no cálculo de campos e energias magnetostáticas. Foram verificados alguns resultados conhecidos para partículas
com magnetização uniforme e o procedimento foi generalizado para o caso de sistemas
periódicos. No caso particular de uma rede de discos uniformemente magnetizados, foi
calculada a distribuição do campo magnético no interior de cada disco e desenvolvido um
tratamento de dinâmicas de precessão magnética quase uniforme para esta distribuição
dos campos internos. A diferença entre as distribuições de campos de ressonância locais
em duas geometrias características permite explicar a anisotropia verificada em experiências de ressonância ferromagnética. Foi ainda introduzido um procedimento variacional
iterativo para calcular deformações na magnetização de partículas com forma de discos
devido tanto a esta forma como à interacção com outras partículas idênticas, deforma-
108
6.4 Conclusões
ção esta que deverá contribuir para o aumento da anisotropia relativamente à verificada
quando considerámos a magnetização uniforme. A introdução desta deformação torna
necessário o abandono da aproximação de Kittel à condição de ressonância, pelo que deverá ser necessária a integração numérica das equações de movimento para reproduzir
quantitativamente os resultados experimentais.
109
6 Partículas de tamanho finito
110
7 Considerações finais
Esta dissertação reflecte alguns dos aspectos mais relevantes do trabalho científico levado
a cabo pelo autor durante o período de Doutoramento. Segue-se um resumo dos principais
resultados.
• Foi apresentado um método geral de expansão de altas temperaturas que permite
a validação nesse regime de simulações de Monte Carlo para sistemas de momentos
magnéticos com interacção arbitrária, e foi proposta uma forma diagramática que
pode permitir levar esta expansão a ordens superiores.
• Dados os problemas de convergência de simulações de Monte Carlo em sistemas
dipolares, foram aplicados diversos métodos de simulação com vista a ultrapassar
estes problemas. Nesta dissertação, optou-se por uma combinação do método de Parallel Tempering com uma optimização específica da taxa de aceitação que permitiu
uma redução significativa dos tempos de autocorrelação e consequente melhoria da
precisão das simulações.
• Estes métodos permitiram caracterizar o comportamento crítico de redes regulares
de dipolos pontuais e as correlações magnéticas em sistemas desordenados, contribuindo para um conhecimento mais completo dos efeitos colectivos com aplicações
em diversas áreas da nanotecnologia.
• No sentido de descrever melhor os sistemas experimentais, foi levantada a aproximação de dipolos pontuais. Combinando um método geral de introdução da forma
das partículas no cálculo dos campos magnetostáticos, foi calculada a distribuição
de campos no interior de discos dispostos numa rede quadrada, em função do espaçamento entre discos e da magnitude e orientação do campo aplicado, assumindo
discos uniformemente magnetizados, o que permitiu explicar a anisotropia encontrada em experiências de ressonância ferromagnética. Utilizando um procedimento
variacional, foi ainda calculada a deformação da magnetização no interior dos discos, o que mostrou a existência de uma competição entre o campo desmagnetizante
devido à forma das partículas e o campo de interacção com os outros discos da rede.
111
7 Considerações finais
Naturalmente, estes resultados levantaram outras questões que ficaram em aberto. Algumas formas de continuar este trabalho seriam:
• Automatização do cálculo da expansão de altas temperaturas para ordens superiores
e para outras grandezas, de forma a poder validar, no regime de altas temperaturas,
mais resultados das simulações de Monte Carlo.
• Optimização das simulações de Parallel Tempering, através da maximização do número de viagens de ida-e-volta nas temperaturas e/ou através da utilização de algoritmos de cluster, de forma a permitir a simulação de sistemas maiores.
• Aplicação destes métodos ao estudo de diferentes redes ordenadas de dipolos pontuais, e introdução gradual de desordem posicional nessas redes.
• Modificação da interacção entre partículas para ter em conta as correcções de tamanho finito das partículas.
• Tratamento completo da dinâmica magnética na rede quadrada de discos, abandonando a aproximação de Kittel e incluindo tanto os elementos não-diagonais do
tensor desmagnetizante como a deformação da magnetização no interior dos discos.
Para além dos resultados científicos, resultaram também destes trabalhos contribuições
técnicas para a capacidade computacional do Centro de Física do Porto, como a implementação e administração de clusters para computação científica, e o aprofundamento de
técnicas de programação, optimização e paralelização e algoritmos de computação científica. Também nestes aspectos mais técnicos, muitos aspectos podem ser melhorados,
entre eles:
• Reformulação do código, tendo em conta metodologias de programação orientada
por objectos como os padrões de design, com vista a criar uma framework de simulação que permita experimentar de forma mais ágil diferentes algoritmos de simulação.
• Utilização de formatos de dados mais apropriados à representação de séries temporais e campos escalares e vectoriais, incluindo mais meta-dados relativos à simulação.
• Paralelização híbrida (tanto em memória partilhada como distribuída) das simulações, com vista à simulação de maiores tamanhos de rede.
• Aceleração por hardware, com recurso a processadores gráficos, da rotina de actualização de campos (a mais intensiva computacionalmente), com vista à simulação
de maiores tamanhos de rede.
112
Apêndices
113
7 Considerações finais
114
A Médias livres
Para calcular médias num sistema com interacções à custa de médias num sistema sem
interacções, é necessário calcular exactamente as médias livres
dµ O e−βHZ
hOi0 = −βH
dµ e
Z
(A.1)
e as suas derivadas em ordem aos campos locais. Para tal vamos escrever o campo em cada
partícula como uma componente constante, o campo externo, mais pequenas variações
em torno desse campo,
H i = H + hi
Figura A.1: Ângulo entre o momento e o campo externo
Se definirmos o ângulo θHµ como o ângulo que o momento µi faz com o campo externo
H (figura A.1), o produto escalar no Hamiltoniano de Zeeman pode escrever-se, em
qualquer sistema de coordenadas, como
µi · H i = µi H cos θHµ + µi · hi
e a exponencial
eβµi ·H i = eβµi H cos θHµ eβµi ·hi
115
A Médias livres
pode ser expandida, para hαi H, e escrever-se
eβµi H cos θHµ (1 + βµi · hi + ...)
Substituindo em A.1
hOi i0,H+h '
dµ Oi eβµi H cos θHµ + β
dµ eβµi H cos θHµ + β
dµ Oi µi eβµi H cos θHµ · hi
dµ µi eβµi H cos θHµ · hi
.
Dado que
a bc − ad
a+b·h
' +
·h
c+d·h
c
c2
podemos escrever, em primeira ordem,
hOi i0,H+h ' hOi i0 + βhαi {hOi µαi i0 − hOi i0 hµαi i0 } .
A derivada em ordem ao campo, no limite em que as variações tendem para zero, é
dada então por
∂
hAi i0,h→0 ' β {hAi µαi i0 − hAi i0 hµαi i0 } .
α
∂hi
As derivadas que encontrámos na expansão da energia livre até segunda ordem podem
então ser escritas como
∂
κ α
hµκk i0 ≡ δik ∂α mκi = δik β (mκα
i − mi mi )
∂hαi
∂ D κ λE
κλ
κλα
κλ α
m
µ
µ
≡
δ
∂
m
=
δ
β
m
−
m
ik
α
ik
i
i
i
i .
k k 0
∂hαi
Estes integrais vão depender do grau de liberdade dos momentos, isto é, se estão confinados a um eixo (Ising), a um plano (XY) ou livres para rodar no espaço (Heisenberg).
Ao longo deste trabalho, considerámos sempre momentos livres para rodar no espaço.
Escolhendo coordenadas esféricas com inclinação θ relativamente à direcção do campo e
azimute φ no plano perpendicular a essa direcção, a média livre de qualquer observável
de uma única partícula (sabendo que podemos sempre separar médias livres em médias
de uma única partícula) é sempre dada por
116
π 2π
hOi i0 =
Oi eβhi cos θ sin θdθdφ
π 0 2π
0
0
0
eβhi cos θ sin θdθdφ
(A.2)
1 2π
Oi eβhi u dφdu
−1 0
1 2π
eβhi u dφdu
−1 0
=
onde fizemos a mudança de variável u = cos θ.
Quando, como nos casos que nos interessam, a observável O depende apenas de combinações de componentes dos momentos, que se escrevem à custa de senos e cossenos de θ e
φ, estes integrais têm sempre solução analítica. Em particular, até ao termo de segunda
ordem é necessário calcular
1 2π
k
µi
0
D
E
u eβhi u dφdu
1
≡ L(βhi )
= −11 0 2π
= coth(βhi ) −
βhi u dφdu
βh
i
e
−1 0
1 2π
D
E
µ⊥
i 0
=
−1
0
cos φeβhi u dφdu
1 2π
−1
1 2π
=
−1
0
0
sin φeβhi u dφdu
1 2π
−1
eβhi u dφdu
0
eβhi u dφdu
= 0
k 2
µi
1 2π
0
=
=
u2 eβhi u dφdu
−1 0
1 2π
−1 0
eβhi u dφdu
L(βhi )
1−2
βhi
!
117
A Médias livres
2
µ⊥
i
1 2π
=
0
−1
10
−1
2π
0
1 2π
=
cos2 φ eβhi u dφdu
−1
sin2 φ eβhi u dφdu
10 2π
−1
eβhi u dφdu
0
eβhi u dφdu
L(βhi )
βhi
=
e assim sucessivamente. As expressões para as médias livres envolvendo produtos de
uma, duas e três componentes podem ser resumidas da seguinte forma:
Uma componente
• Média livre da componente paralela ao campo
k
D
µi
E
0
= L(βH)
• Média livre da componente perpendicular ao campo
D
µ⊥
i
E
=0
0
Produto de duas componentes
• Média livre do quadrado da componente paralela ao campo
k 2
µi
0
=1−2
L(βH)
βH
• Média livre do quadrado da componente perpendicular ao campo
µ⊥
i
2 0
=
L(βH)
βH
• Termos cruzados
D
k
µi µ⊥
i
E
0
=0
Produto de três componentes
• Média livre do cubo da componente paralela ao campo
118
k 3
µi
0
=
−2βH + (6 + H 2 β 2 ) L(βH)
β 2H 2
• Média livre do produto da componente paralela com o quadrado de uma componente
perpendicular
k
µi
2
µ⊥
i
0
=
Hβ − 3L(βH)
β 2H 2
• Média livre do produto do quadrado da componente paralela com uma componente
perpendicular
k 2
µi
µ⊥
i
0
=0
• Em todos os restantes casos, em que não existe nenhuma componente paralela ao
plano, o resultado é nulo.
Genericamente, estes integrais podem escrever-se como
2π
=
0
k n
µi
m
1
µ⊥
i
(cos φ)m (sin φ)l dφ
2π
1
−1
l
2
µ⊥
i
0
=
m+l
un (1 − u2 )
eβHu du
2 sinh (βH) / (βH)
1
2
onde µ⊥
e µ⊥
são duas componentes ortogonais no plano perpendicular ao campo aplii
i
cado.
119
A Médias livres
120
B Amplitudes de forma
Para exemplificar o formalismo apresentado no capítulo 6, calculamos aqui as amplitudes
de forma de uma esfera e de um cilindro, bem como as respectivas energias magnetostáticas
no caso de corpos uniformemente magnetizados.
Esfera
Figura B.1: Representação esquemática de uma esfera uniformemente magnetizada
A amplitude de forma de uma esfera de raio R (figura B.1) é obtida integrando
D(k) =
e−ik·r dr
(B.1)
V
em coordenadas esféricas, resultando em
D(k) =
4πR2
j1 (kR),
k
(B.2)
onde jn (x) é a função esférica de Bessel, relacionada com a mais comum função cilíndrica
de Bessel Jn (x) por
r
π
jn (x) =
Jn+1/2 (x)
(B.3)
2x
121
B Amplitudes de forma
A energia magnetostática de uma esfera uniformemente magnetizada resulta de substituir a amplitude de forma correspondente na equação 6.17 e integrar mais uma vez em
coordenadas esféricas,
2π
µ0 M 2 R 3 .
(B.4)
Em =
9
O mesmo formalismo permite mostrar que o campo no exterior de uma esfera uniformemente magnetizada corresponde exactamente ao campo produzido por um dipolo pontual
[147].
Cilindro
Já se a partícula for de forma cilíndrica de altura t e raio R , a sua amplitude de forma é
dada por
D(k) =
e−ik·r dr
V
! t/2
=
e−iqz dz
!
e−ig·ρ dρ
(B.5)
−t/2
= 4πR
sin
q
t
q
2
J1 (gR)
,
g
onde mudámos para coordenadas cilíndricas tanto no espaço real, r = (ρ cos θ, ρ sin θ, z),
como no espaço de Fourier, k = (g cos φ, g sin φ, q).
a)
b)
Figura B.2: Representação esquemática de um cilindro uniformemente magnetizado em
a) configuração axial e b) configuração planar.
122
Neste caso, a energia magnetostática
µ0 M 2
Em =
16π 3
|D(k)|2
(m̂ · k)2 dk
2
k
(B.6)
depende da orientação da magnetização, através do termo (m̂ · k)2 , pelo que vamos considerar dois casos em separado: magnetização paralela ou perpendicular ao eixo do cilindro
(figura B.2).
Configuração axial
No caso de magnetização paralela ao eixo do cilindro (figura B.2a), a que iremos chamar
configuração axial, o termo (m̂ · k)2 é dado simplesmente por q 2 e o integral em q escrevese
!
∞
πt 1 − e−gt
sin2 (qt/2)
.
(B.7)
dq =
2
2
2
gt
−∞ (g + q )
Neste caso, a energia magnetostática é dada por
µ0 M 2 2
R (8R + 3πt) + ,
6 #
"
"
2 #!
2
R
R
2
3
2
3
− 4R t + t K −4
,
+ −4R t + t E −4
t
t
axial
Em
=
(B.8)
onde K é o integral elíptico completo de primeira ordem,
∞
(2n − 1)!!
πX
K[x] =
2 n=0
2n!!
!2
(B.9)
x2n
e E o integral elíptico completo de segunda ordem,


!2
∞
X
(2n − 1)!!
x2n 
π
E[x] =
1−
.
2
2n!!
2n − 1 
n=1
(B.10)
Configuração planar
No caso de magnetização perpendicular ao eixo do cilindro (figura B.2b), a que iremos
chamar configuração planar (no plano das faces do cilindro), o integral em q resulta
∞
−∞
sin2 (qt/2)
πt
1 − e−gt
dq
=
1
−
q 2 (g 2 + q 2 )
2g 2
gt
!
(B.11)
123
B Amplitudes de forma
e a energia magnetostática é dada por
µ0 M 2 3
8R +
# 12
planar
Em
=−
+ −4R t + t
2
3
"
R
E −4
t
2
− 4R t + t
2
3
(B.12)
"
K −4
R
t
2 #!
o que permite verificar, numericamente, que para t . 1.81R, a configuração planar tem
uma energia menor à da configuração axial.
124
C Lista de publicações
O resultado visível, em termos de publicações, do trabalho científico levado a cabo durante
o período de Doutoramento, incluindo colaborações em trabalhos para além dos que estão
descritos nesta dissertação, pode ser encontrado nos seguintes artigos1 :
Métodos de Monte Carlo
Optimized multicanonical simulations: A proposal based on classical
fluctuation theory
• Autores: Lopes, JV; Costa, MD; Lopes dos Santos, JMB; Toral, R
• Publicação: Physical Review E 74, 046702 (2006)
Abstract: “We propose a recursive procedure to estimate the microcanonical density of
states in multicanonical Monte Carlo simulations which relies only on measurements of
moments of the energy distribution, avoiding entirely the need for energy histograms.
This method yields directly a piecewise analytical approximation to the microcanonical
inverse temperature β(E) and allows improved control over the statistics and efficiency of
the simulations. We demonstrate its utility in connection with recently proposed schemes
for improving the efficiency of multicanonical sampling, either with adjustment of the
asymptotic energy distribution or with the replacement of single spin flip dynamics with
collective updates.“
Analytical study of tunneling times in flat histogram Monte Carlo
• Autores: Costa, MD; Lopes, JV; Lopes dos Santos, JMB
• Publicação: Europhysics Letters 72, 5, 802-808 (2005)
1
A lista actualizada de publicações pode ser encontrada aqui: http://www.researcherid.com/rid/
A-6968-2008
125
C Lista de publicações
Abstract: “We present a model for the dynamics in energy space of multicanonical simulation methods that lends itself to a rather complete analytic characterization. The
dynamics is completely determined by the density of states. In the ±J 2D spin glass the
transitions between the ground state level and the first excited one control the long time
dynamics. We are able to calculate the distribution of tunneling times and relate it to the
equilibration time of a starting probability distribution. In this model, and possibly in
any model in which entering and exiting regions with low density of states are the slowest
processes in the simulations, tunneling time can be much larger (by a factor of O(N))
than the equilibration time of the probability distribution. We find that these features
also hold for the energy projection of single spin flip dynamics.“
Modelos com Interacção Dipolar
Apparent phase transition in a chain of magnetic dipoles?
• Autores: Pogorelov, YG; Costa, MD
• Publicação: Journal Of Magnetism And Magnetic Materials 272, E979-E980 (2004)
Abstract: “Considering chains of classical spins coupled by long-range interactions ∝ r−σ ,
we show a qualitative difference between the critical behaviour of systems with isotropic
and anisotropic interactions, as seen from the analysis of the Binder cumulant by Monte
Carlo simulations. Unlike the known absence of long range order in isotropic case for
sv>2, we found indications of finite critical temperatures in the anisotropic case even for
sv as high as 3, related to dipolar coupling.”
Continuum model for dipolar coupled planar lattices
• Autores: Costa, MD; Pogorelov, YG
• Publicação: Journal Of Magnetism And Magnetic Materials 258, 32-34 (2003)
Abstract: “In an effective continuum approach alike the phenomenological Landau theory,
we study low energy excitations in a square lattice of dipolar coupled magnetic moments
µ, over continuously degenerate microvortex (MV) ground states defined by an arbitrary
angle 0 < θ < π/2. We consider two vector order parameters: the MV vector v =
µ(cos θ, sin θ) and the ferromagnetic (FM) vector f = 21 (∂y vx , −∂x vy ). The excitation
energy density ∼ f 2 leads to a non-linear Euler equation. It allows, besides common linear
waves of small amplitude, also non-linear excitations with unlimited (but slow) variation of
126
θ(r). For plane wave excitations θ(r) = θ(n · r) propagating along n = (cos φ, sin φ), exact
integrals of Euler equation are found. The density of excitation states turns anisotropic
in θ, conforming to the enhanced occurrence of MV-like states with θ close to 0 or π/2 in
our Monte Carlo simulations of this system at low excitation energies.“
Criteria for long-range magnetic order in planar lattices of dipolar
coupled magnetic nanoparticles
• Autores: Costa, MD; Pogorelov, YG
• Publicação: Physica Status Solidi A-Applied Research 189, 3, 923-927 (2002)
Abstract: “The conditions are sought for the dipolar interactions in a planar lattice of
separate magnetic nanogranules to bring this system from the initially frustrated ground
state to a long-range ferromagnetic (FM) order. This is promoted by lowering the symmetry: either of the lattice (e.g., from square to rectangular (a > b)) for point-like or
spherical granules, or of a finite-size granule (e.g., from spherical to rod-like or to disk-like)
on square lattice. The FM order turns possible after some critical values of asymmetry
parameters, though the system energy may behave non-analytically in these parameters.
The comparison is done with the available experimental data.”
Interpretação de Resultados Experimentais
Collective dynamics and ferromagnetic order in random planar arrays
of magnetic granules
• Autores: Pogorelov, YG; Kakazei, GN; Costa, MD; Sousa, JB
• Publicação: Journal Of Applied Physics 103, 07B723 (2008)
Abstract: “A dynamical study is done on existence and stability of ferromagnetically ordered ground state in a positionally disordered planar array of magnetic moments coupled
only by dipolar forces. Starting from almost aligned ground state under a strong enough
applied field, the excitation energy spectrum and related eigenmodes are found, permitting to develop a mean-field analysis of the static magnetization in function of magnetic
field and temperature. In the limit of zero applied field, the stability conditions are obtained for the onset of in-plane spontaneous magnetization, manifesting a specific ’order
from disorder’ mechanism.”
127
C Lista de publicações
Probing arrays of circular magnetic microdots by ferromagnetic
resonance
• Autores: Kakazei, GN; Mewes, T; Wigen, PE; Hammel, PC; Slavin, AN; Pogorelov,
YG; Costa, MD; Golub, VO; Guslienko, KYu; Novosad, V
• Publicação: Journal Of Nanoscience And Nanotechnology 8, 6, 2811-2826 (2008)
Abstract: “X-band ferromagnetic resonance (FMR) was used to characterize in-plane
magnetic anisotropies in rectangular and square arrays of circular nickel and Permalloy
microdots. In the case of a rectangular lattice, as interdot distances in one direction
decrease, the in-plane uniaxial anisotropy field increases, in good agreement with a simple theory of magnetostatically interacting uniformly magnetized dots. In the case of a
square lattice a four-fold anisotropy of the in-plane FMR field H(r) was found when the
interdot distance a gets comparable to the dot diameter D. This anisotropy, not expected
in the case of uniformly magnetized dots, was explained by a non-uniform magnetization
m(r) in a dot in response to dipolar forces in the patterned magnetic structure. It is
well described by an iterative solution of a continuous variation procedure. In the case
of perpendicular magnetization multiple sharp resonance peaks were observed below the
main FMR peak in all the samples, and the relative positions of these peaks were independent of the interdot separations. Quantitative description of the observed multiresonance
FMR spectra was given using the dipole-exchange spin wave dispersion equation for a
perpendicularly magnetized film where in-plane wave vector is quantized due to the finite
dot radius, and the inhomogenetiy of the intradot static demagnetization field in the nonellipsoidal dot is taken into account. It was demonstrated that ferromagnetic resonance
force microscopy (FMRFM) can be used to determine both local and global properties of
patterned submicron ferromagnetic samples. Local spectroscopy together with the possibility to vary the tip-sample spacing enables the separation of those two contributions to
a FMRFM spectrum. The global FMR properties of circular submicron dots determined
using magnetic resonance force microscopy are in a good agreement with results obtained
using conventional FMR and with theoretical descriptions.“
Origin of fourfold anisotropy in square lattices of circular
ferromagnetic dots
• Autores: Kakazei, GN; Pogorelov, YG; Costa, MD; Mewes, T; Wigen, PE; Hammel,
PC; Golub, VO; Okuno, T; Novosad, V
• Publicação: Physical Review B 74 060406(R) (2006)
128
Abstract: “We discuss the fourfold anisotropy of the in-plane ferromagnetic resonance
field Hr, found in a square lattice of circular Permalloy dots when the interdot distance
a becomes comparable to the dot diameter d. The minimum Hr along the lattice h11i
axes and the maximum along the h10i axes differ by ~50 Oe at a/d=1.1. This anisotropy,
not expected in uniformly magnetized dots, is explained by a mechanism of nonuniform
magnetization m(r) in a dot in response to dipolar forces in the patterned magnetic
structure under strong enough applied field. It is well described by an iterative solution
of a continuous variational procedure.”
Interlayer dipolar interactions in multilayered granular films
• Autores: Kakazei, GN; Pogorelov, YG; Costa, MD; Golub, VO; Sousa, JB; Freitas,
PP; Cardoso, S; Wigen, PE
• Publicação: Journal Of Applied Physics 97, 10A723-10A723-3 (2005)
Abstract: “Using the ferromagnetic resonance (FMR) technique, the interlayer coupling
is studied in granular multilayered system (CoFe/Al2O3)n with varying number of layers
n=1,...,8. The main difference between the cases of n=1 and n>=3 is a notable splitting
of FMR line under external field normal to layers. This is explained by an interlayer
dipolar coupling, only possible for discontinuous layers. The relative magnitudes of two
absorption peaks and distance between them are described by a simple model of planar
square lattices of magnetic dipoles.“
129
C Lista de publicações
130
Bibliografia
[1] R. P. Cowburn. Property variation with shape in magnetic nanoelements. Journal
of Physics D: Applied Physics, 33(1):R1–R16, 2000.
[2] K. Yu. Guslienko. Magnetostatic interdot coupling in two-dimensional magnetic dot
arrays. Applied Physics Letters, 75(3):394–396, July 1999.
[3] Philip Moriarty.
64(3):297, 2001.
Nanostructured materials.
Reports on Progress in Physics,
[4] S.D. Bader. Magnetism in low dimensionality. Surface Science, 500(1-3):172–188,
March 2002.
[5] S. D. Bader. Colloquium: Opportunities in nanomagnetism. Reviews of Modern
Physics, 78(1):1–15, 2006.
[6] R. Skomski. Nanomagnetics. Journal of Physics: Condensed Matter, 15(20):R841,
2003.
[7] Gary A. Prinz. Magnetoelectronics. Science, 282(5394):1660–1663, November 1998.
[8] Ca Ross. Patterned magnetic recording media. Annual Review of Materials Research, 31(1):203–235, 2001.
[9] Andreas Moser, Kentaro Takano, David T Margulies, Manfred Albrecht, Yoshiaki
Sonobe, Yoshihiro Ikeda, Shouheng Sun, and Eric E Fullerton. Magnetic recording:
advancing into the future. Journal of Physics D: Applied Physics, 35(19):R157,
2002.
[10] A O Adeyeye and N. Singh. Large area patterned magnetic nanostructures. Journal
of Physics D: Applied Physics, 41(15):153001, 2008.
[11] C.L. Chien, F.Q. Zhu, and J.G. Zhu. Patterned nanomagnets. Physics Today,
60(6):40–45, June 2007.
131
Bibliografia
[12] S. Tehrani, E. Chen, M. Durlam, M. DeHerrera, J. M. Slaughter, J. Shi, and
G. Kerszykowski. High density submicron magnetoresistive random access memory
(invited). In J. Appl. Phys., volume 85, pages 5822–5827. AIP, April 1999.
[13] S. Tehrani, J.M. Slaughter, M. Deherrera, B.N. Engel, N.D. Rizzo, J. Salter, M. Durlam, R.W. Dave, J. Janesky, B. Butcher, K. Smith, and G. Grynkewich. Magnetoresistive random access memory using magnetic tunnel junctions. Proceedings of
the IEEE, 91(5):703–714, 2003.
[14] R. P. Cowburn and M. E. Welland. Room temperature magnetic quantum cellular
automata. Science, 287(5457):1466–1468, February 2000.
[15] G. Csaba, A. Imre, G.H. Bernstein, W. Porod, and V. Metlushko. Nanocomputing
by field-coupled nanomagnets. Nanotechnology, IEEE Transactions on, 1(4):209–
213, 2002.
[16] D. A. Allwood, Gang Xiong, M. D. Cooke, C. C. Faulkner, D. Atkinson, N. Vernier,
and R. P. Cowburn. Submicrometer ferromagnetic NOT gate and shift register.
Science, 296(5575):2003–2006, June 2002.
[17] A. Imre, G. Csaba, L. Ji, A. Orlov, G. H. Bernstein, and W. Porod. Majority logic
gate for magnetic Quantum-Dot cellular automata. Science, 311(5758):205–208,
2006.
[18] Briefing
2:
Nanotechnology
and
http://www.foresight.org/updates/Briefing2.html.
enabling
technologies.
[19] A. E. Berkowitz, J. R. Mitchell, M. J. Carey, A. P. Young, S. Zhang, F. E. Spada,
F. T. Parker, A. Hutten, and G. Thomas. Giant magnetoresistance in heterogeneous
Cu-Co alloys. Physical Review Letters, 68(25):3745, June 1992.
[20] M. N. Baibich, J. M. Broto, A. Fert, F. Nguyen Van Dau, F. Petroff, P. Etienne, G. Creuzet, A. Friederich, and J. Chazelas. Giant magnetoresistance of
(001)Fe/(001)Cr magnetic superlattices. Physical Review Letters, 61(21):2472, November 1988. Copyright (C) 2009 The American Physical Society; Please report
any problems to [email protected].
[21] Richard P. Feynman. There’s plenty of room at the bottom. Engineering and
Science, 23(5), 1960.
132
Bibliografia
[22] K. Eric Drexler. Molecular engineering: An approach to the development of general
capabilities for molecular manipulation. Proceedings of the National Academy of
Sciences of the United States of America, 78(9):5275–5278, 1981.
[23] G. Binnig, H. Rohrer, Ch. Gerber, and E. Weibel. Surface studies by scanning
tunneling microscopy. Physical Review Letters, 49(1):57, July 1982.
[24] G. Binnig, C. F. Quate, and Ch. Gerber. Atomic force microscope. Physical Review
Letters, 56(9):930, March 1986.
[25] Y. Martin and H. K. Wickramasinghe. Magnetic imaging by “force microscopy” with
1000 [A-ring] resolution. Applied Physics Letters, 50(20):1455–1457, May 1987.
[26] J. H. E. Griffiths. Anomalous high-frequency resistance of ferromagnetic metals.
Nature, 158(4019):670–671, 1946.
[27] Charles Kittel. On the theory of ferromagnetic resonance absorption. Physical
Review, 73(2):155, 1948.
[28] John David Jackson. Classical Electrodynamics Third Edition. Wiley, 3 edition,
August 1998.
[29] Ralph Skomski. Simple models of magnetism. Oxford University Press, 2008.
[30] J. M. Luttinger and L. Tisza. Theory of dipole interaction in crystals. Physical
Review, 70(11-12):954, December 1946.
[31] E. J. S Lage. Física Estatística. Manuais universitários. Fundação Calouste Gulbenkian, Lisboa, 1995.
[32] Ryogo Kubo. Generalized cumulant expansion method. Journal of the Physical
Society of Japan, 17:1100–1120, 1962.
[33] R. A. Fisher and J. Wishart. The derivation of the pattern formulae of Two-Way
partitions from those of simpler patterns. Proceedings of the London Mathematical
Society, 2(1):195, 1932.
[34] Donald L. Kreher and Douglas Robert Stinson. Combinatorial algorithms. CRC
Press, 1999.
[35] Noboru Fukushima. A new method of the high temperature series expansion. Journal of Statistical Physics, 111(5):1049–1090, June 2003.
133
Bibliografia
[36] Faà di Bruno. Sullo sviluppo delle funzioni. Annali di Scienze Matematiche e
Fisiche, 6:479–480, 1855.
[37] Mark E. J. Newman and G. T. Barkema. Monte Carlo methods in statistical physics.
Oxford University Press, 1999.
[38] Nicholas Metropolis, Arianna W. Rosenbluth, Marshall N. Rosenbluth, Augusta H.
Teller, and Edward Teller. Equation of state calculations by fast computing machines. The Journal of Chemical Physics, 21(6):1087–1092, June 1953.
[39] G. M. Torrie and J. P. Valleau. Nonphysical sampling distributions in monte carlo
free-energy estimation: Umbrella sampling. Journal of Computational Physics,
23(2):187–199, February 1977.
[40] Jooyoung Lee. New monte carlo algorithm: Entropic sampling. Physical Review
Letters, 71(2):211, July 1993.
[41] Bernd A. Berg and Tarik Celik. New approach to spin-glass simulations. Physical
Review Letters, 69(15):2292, October 1992.
[42] Bernd A. Berg and Thomas Neuhaus. Multicanonical ensemble: A new approach
to simulate first-order phase transitions. Physical Review Letters, 68(1):9, 1992.
[43] Fugao Wang and D. P. Landau. Efficient, Multiple-Range random walk algorithm
to calculate the density of states. Physical Review Letters, 86(10):2050, March 2001.
[44] Miguel Dias Costa. Estados magnéticos de sistemas com interacção dipolar baseados
em redes planares. PhD thesis, Departamento de Física da FCUP, 2001. Supervisor:
Yuri G. Pogorelov.
[45] Robert H. Swendsen and Jian-Sheng Wang. Replica monte carlo simulation of SpinGlasses. Physical Review Letters, 57(21):2607, November 1986.
[46] E. Marinari and G. Parisi. Simulated tempering: a new monte carlo scheme. Europhys. Lett, 19(6):451–458, 1992.
[47] A. P. Lyubartsev, A. A. Martsinovski, S. V. Shevkunov, and P. N. VorontsovVelyaminov. New approach to monte carlo calculation of the free energy: Method of
expanded ensembles. The Journal of Chemical Physics, 96(3):1776–1783, February
1992.
134
Bibliografia
[48] Koji Hukushima and Koji Nemoto. Exchange monte carlo method and application
to spin glass simulations. Journal of the Physical Society of Japan, 65:1604–1608,
1996.
[49] Elmar Bittner, Andreas Nussbaumer, and Wolfhard Janke. Make life simple: Unleash the full power of the parallel tempering algorithm. Physical Review Letters,
101(13):130603–4, 2008.
[50] Joao Manuel Viana Parente Lopes. Métodos de Monte Carlo para sistemas com
desordem e interacção de longo alcance. Tese de Doutoramento. 2006.
[51] Helmut G Katzgraber, Simon Trebst, David A Huse, and Matthias Troyer.
Feedback-optimized parallel tempering monte carlo. cond-mat/0602085, February
2006. J. Stat. Mech. P03018 (2006).
[52] Robert H. Swendsen and Jian-Sheng Wang. Nonuniversal critical dynamics in monte
carlo simulations. Physical Review Letters, 58(2):86, 1987.
[53] Ulli Wolff. Collective monte carlo updating for spin systems. Physical Review
Letters, 62(4):361, 1989.
[54] Kouki Fukui and Synge Todo. Order-N cluster monte carlo method for spin systems
with long-range interactions. Journal of Computational Physics, 228(7):2629–2642,
April 2009.
[55] Kurt Binder and Dieter W. Heermann. Monte Carlo Simulation in Statistical Physics. Springer, 4th, enlarged ed. edition, August 2002.
[56] H. G Katzgraber. Scientific software engineering in a nutshell. Imprint, 2009.
[57] Free software foundation. http://www.fsf.org/.
[58] Eclipse.org. http://www.eclipse.org/.
[59] KDevelop. http://www.kdevelop.org/.
[60] Anjuta. http://projects.gnome.org/anjuta/index.shtml.
[61] Subversion. http://subversion.tigris.org/.
[62] Git. http://git-scm.com/.
[63] Mercurial. http://mercurial.selenic.com/wiki/.
135
Bibliografia
[64] Redmine. http://www.redmine.org/.
[65] The trac project. http://trac.edgewall.org/.
[66] GDB: the GNU project debugger. http://www.gnu.org/software/gdb/.
[67] Valgrind home. http://valgrind.org/.
[68] GNU gprof. http://www.cs.utah.edu/dept/old/texinfo/as/gprof_toc.html.
[69] KCachegrind. http://kcachegrind.sourceforge.net/html/Home.html.
[70] Autoconf. http://www.gnu.org/software/autoconf/.
[71] SCons: a software construction tool. http://www.scons.org/.
[72] GCC, the GNU compiler collection. http://gcc.gnu.org/.
[73] Intel® compilers. http://software.intel.com/en-us/intel-compilers/.
[74] GSL - GNU scientific library. http://www.gnu.org/software/gsl/.
[75] Intel® MKL. http://software.intel.com/en-us/intel-mkl/.
[76] AMD core math library. http://developer.amd.com/cpu/Libraries/acml/Pages/default.aspx.
[77] Message passing interface. http://www.mcs.anl.gov/research/projects/mpi/.
[78] Open MPI: open source high performance computing. http://www.open-mpi.org/.
[79] OpenMP. http://openmp.org/wp/.
[80] CUDA zone. http://www.nvidia.com/object/cuda_home.html.
[81] Advanced micro devices, AMD. http://ati.amd.com/technology/streamcomputing/.
[82] The community ENTerprise operating system. http://www.centos.org/.
[83] Cluster software portal. http://www.perceus.org/portal/.
[84] xCAT - extreme cluster administration toolkit. http://xcat.sourceforge.net/.
[85] Sun grid engine. http://gridengine.sunsource.net/.
[86] Cluster resources. http://www.clusterresources.com/products/torque/.
[87] Maui cluster scheduler. http://www.clusterresources.com/products/maui-clusterscheduler.php/.
136
Bibliografia
[88] Ganglia monitoring system. http://ganglia.sourceforge.net/.
[89] Nagios. http://www.nagios.org/.
[90] The multi router traffic grapher. http://oss.oetiker.ch/mrtg/.
[91] NetCDF (network common data form). http://www.unidata.ucar.edu/software/netcdf/.
[92] HDF group - HDF5. http://www.hdfgroup.org/HDF5/.
[93] Grace. http://plasma-gate.weizmann.ac.il/Grace/.
[94] Gnuplot. http://www.gnuplot.info/.
[95] matplotlib. http://matplotlib.sourceforge.net/.
[96] VTK - the visualization toolkit. http://www.vtk.org/.
[97] The MayaVi data visualizer. http://mayavi.sourceforge.net/.
[98] Inkscape. http://www.inkscape.org/.
[99] SciPy. http://www.scipy.org/.
[100] Enthought, inc. :: Scientific computing solutions. http://www.enthought.com/.
[101] Wolfram research:
Mathematica,
http://www.wolfram.com/.
technical
and
scientific
software.
[102] Sage: Open source mathematics software. http://www.sagemath.org/.
[103] Doxygen. http://www.stack.nl/~dimitri/doxygen/.
[104] LyX - the document processor. http://www.lyx.org/.
[105] Zotero. http://www.zotero.org/.
[106] PI Belobrov, RS Gehkt, and VA Ignatchenko. Ground-state on systems with dipole
interactions. Zhurnal Eksperimentalnoi I Teoreticheskoi Fiziki, 84(3):1097–1108,
1983.
[107] PI Belobrov, VA Voevodin, and VA Ignatchenko. Ground-state of a dipole system in a planar rhombic lattice. Zhurnal Eksperimentalnoi I Teoreticheskoi Fiziki,
88(3):889–893, 1985.
137
Bibliografia
[108] J. G. Brankov and D. M. Danchev. Ground state of an infinite two-dimensional
system of dipoles on a lattice with arbitrary rhombicity angle. Physica A: Statistical
and Theoretical Physics, 144(1):128–139, July 1987.
[109] George O. Zimmerman, A. K. Ibrahim, and F. Y. Wu. Planar classical dipolar
system on a honeycomb lattice. Physical Review B, 37(4):2059, February 1988.
[110] Sona Prakash and Christopher L. Henley. Ordering due to disorder in dipolar
magnets on two-dimensional lattices. Physical Review B, 42(10):6574, October 1990.
[111] Freeman J. Dyson. Existence of a phase-transition in a one-dimensional ising ferromagnet. Communications in Mathematical Physics, 12(2):91–107, June 1969.
[112] N. D. Mermin and H. Wagner. Absence of ferromagnetism or antiferromagnetism
in one- or Two-Dimensional isotropic heisenberg models. Physical Review Letters,
17(22):1133, November 1966.
[113] P. Bruno. Absence of spontaneous magnetic order at nonzero temperature in oneand Two-Dimensional heisenberg and XY systems with Long-Range interactions.
Physical Review Letters, 87(13):137203, 2001.
[114] Melvin Lax. Dipoles on a lattice : the spherical model. The Journal of Chemical
Physics, 20(9):1351–1359, 1952.
[115] S. Romano. Computer simulation study of a two-dimensional dipolar lattice. Il
Nuovo Cimento D, 9(4):409–430, April 1987.
[116] S. Romano. Computer-simulation study of a three-dimensional lattice-spin model
with dipolar-type interactions. Physical Review B, 49(17):12287, May 1994.
[117] E. Rastelli and A. Tassi. Absence of long-range order in three-dimensional heisenberg models. Physical Review B, 40(7):5282, 1989.
[118] E. Rastelli, A. Tassi, A. Pimpinelli, and S. Sedazzari. Triangular planar antiferromagnet in an external magnetic field. Physical Review B, 45(14):7936, April 1992.
[119] E. Rastelli, A. Carbognani, S. Regina, and A. Tassi. Pseudodipolar anisotropy in
the square rotator model. In J. Appl. Phys., volume 87, pages 5902–5904. AIP, May
2000.
[120] E. Rastelli, S. Regina, and A. Tassi. Planar triangular model with long-range interactions. Physical Review B, 66(5):054431, 2002.
138
Bibliografia
[121] E. Rastelli, S. Regina, and A. Tassi. Phase diagram of the square planar model
with long range isotropic antiferromagnetic and dipolar interactions. physica status
solidi (a), 189(2):555–558, 2002.
[122] E. Rastelli, S. Regina, and A. Tassi. Short-range exchange and long-range dipole
interactions in a triangular planar model. Physical Review B, 67(9):094429, March
2003.
[123] E. Rastelli, S. Regina, and A. Tassi. Monte carlo simulation for square planar model
with a small fourfold symmetry-breaking field. Physical Review B, 70(17):174447,
November 2004.
[124] E. Rastelli, S. Regina, and A. Tassi. Monte carlo simulations on a triangular ising
antiferromagnet with nearest and next-nearest interactions. Physical Review B,
71(17):174406, May 2005.
[125] E. Rastelli, S. Regina, and A. Tassi. Phase transitions in a square ising model
with exchange and dipole interactions. Physical Review B (Condensed Matter and
Materials Physics), 73(14):144418–11, April 2006.
[126] K. De’Bell, A. B. MacIsaac, I. N. Booth, and J. P. Whitehead. Dipolar-induced
planar anisotropy in ultrathin magnetic films. Physical Review B, 55(22):15108,
June 1997. Copyright (C) 2009 The American Physical Society; Please report any
problems to [email protected].
[127] K. De Bell, A. B. MacIsaac, and J. P. Whitehead. Dipolar effects in magnetic thin
films and quasi-two-dimensional systems. Reviews of Modern Physics, 72(1):225,
2000.
[128] Yusuke Tomita. Monte carlo study of Two-Dimensional heisenberg dipolar lattices.
Journal of the Physical Society of Japan, 78:114004, 2009.
[129] Katsuyoshi Matsushita, Ryoko Sugano, Akiyoshi Kuroda, Yusuke Tomita, and Hajime Takayama. Peculiar ‘from-Edge-to-Interior’ spin freezing in a magnetic dipolar
cube. Journal of the Physical Society of Japan, 74:2651–2654, 2005.
[130] Ryoko Sugano, Katsuyoshi Matsushita, Akiyoshi Kuroda, Yusuke Tomita, and Hajime Takayama. Magnetic properties of Two-Dimensional dipolar squares: Boundary geometry dependence. Journal of the Physical Society of Japan, 76:044705,
2007.
139
Bibliografia
[131] Neal Madras and Alan D. Sokal. The pivot algorithm: A highly efficient monte
carlo method for the self-avoiding walk. Journal of Statistical Physics, 50(1):109–
186, 1988.
[132] Michael E. Fisher and Michael N. Barber. Scaling theory for Finite-Size effects in
the critical region. Physical Review Letters, 28(23):1516, June 1972.
[133] I. A. Campbell, K. Hukushima, and H. Takayama. Extended scaling scheme for
critically divergent quantities in ferromagnets and spin glasses. Physical Review
Letters, 97(11):117202, 2006.
[134] H. Zhang and M. Widom. Spontaneous magnetic order in random dipolar solids.
Physical Review B, 51(14):8951, April 1995.
[135] D. Kechrakos and K. N. Trohidou. Magnetic properties of dipolar interacting singledomain particles. Physical Review B, 58(18):12169, November 1998.
[136] G. Takzei, L. Gun’ko, I. Sych, A. Surzhenko, S. Cherepov, Yu. Troshchenkov, and
I. Mirebeau. Onset of long-range ferromagnetic order in a system of ferromagnetic
particles with giant magnetic moments. Journal of Experimental and Theoretical
Physics, 87(5):1003–1008, November 1998.
[137] S. Denisov. Long-range order and magnetic relaxation in a system of single-domain
particles. Physics of the Solid State, 41(10):1672–1677, October 1999.
[138] S. I. Denisov, V. F. Nefedchenko, and K. N. Trohidou. Dipolar ferromagnetism
in nanoparticle ensembles. Journal of Physics: Condensed Matter, 12:7111–7115,
2000.
[139] S. I. Denisov and K. N. Trohidou. Fluctuation theory of magnetic relaxation for
two-dimensional ensembles of dipolar interacting nanoparticles. Physical Review B,
64(18):184433, October 2001.
[140] G. A. Held, G. Grinstein, H. Doyle, Shouheng Sun, and C. B. Murray. Competing
interactions in dispersions of superparamagnetic nanoparticles. Physical Review B,
64(1):012408, June 2001.
[141] Kazuo Yamamoto, Sara A. Majetich, Martha R. McCartney, Madhur Sachan, Saeki
Yamamuro, and Tsukasa Hirayama. Direct visualization of dipolar ferromagnetic
domain structures in co nanoparticle monolayers by electron holography. Applied
Physics Letters, 93(8):082502–3, 2008.
140
Bibliografia
[142] J. Viana Lopes, J. M. B. Lopes dos Santos, and Yu. G. Pogorelov. Dipolar interactions and anisotropic magnetoresistance in metallic granular systems. Physical
Review B, 66(6):064416, 2002.
[143] A. Aharoni. Magnetostatic energy calculations. Magnetics, IEEE Transactions on,
27(4):3539–3547, 1991.
[144] M. Beleggia, S. Tandon, Y. Zhu, and M. De Graef. On the magnetostatic interactions
between nanoparticles of arbitrary shape. Journal of Magnetism and Magnetic
Materials, 278(1-2):270–284, July 2004.
[145] M. Beleggia and M. De Graef. General magnetostatic shape-shape interactions.
Journal of Magnetism and Magnetic Materials, 285(1-2):L1–L10, 2005.
[146] Marco Beleggia, Shakul Tandon, Yimei Zhu, and Marc De Graef. On the computation of the demagnetization tensor for particles of arbitrary shape. Journal
of Magnetism and Magnetic Materials, 272-276(Supplement 1):E1197–E1199, May
2004.
[147] M. Beleggia and M. De Graef. On the computation of the demagnetization tensor
field for an arbitrary particle shape using a fourier space approach. Journal of
Magnetism and Magnetic Materials, 263(1-2):L1–L9, July 2003.
[148] J. A. Osborn. Demagnetizing factors of the general ellipsoid. Physical Review,
67(11-12):351, June 1945.
[149] R. I. Joseph and E. Schlomann. Demagnetizing field in nonellipsoidal bodies. Journal of Applied Physics, 36(5):1579–1593, May 1965.
[150] G. N. Kakazei, Yu. G. Pogorelov, M. D. Costa, T. Mewes, P. E. Wigen, P. C.
Hammel, V. O. Golub, T. Okuno, and V. Novosad. Origin of fourfold anisotropy in
square lattices of circular ferromagnetic dots. Physical Review B, 74(6), 2006.
[151] K. Yu. Guslienko. Magnetic anisotropy in two-dimensional dot arrays induced by
magnetostatic interdot coupling. Physics Letters A, 278(5):293–298, 2001.
[152] Konstantin L. Metlov. Quasiuniform In-Plane magnetization state of thin cylindrical dots in a square array and related anisotropy. Physical Review Letters,
97(12):127205–4, 2006.
141