Costa, Lidiane da Silva Araújo
Comportamento crítico de filmes Ising � S=1
em estruturas cúbicas simples / Lidiane da Silva
Araújo Costa. � Recife : O Autor, 2008.
xiv, 61 folhas: il., fig., tab.
Dissertação (mestrado) - Universidade Federal de
Pernambuco. CCEN. Física, 2008.
Inclui bibliografia.
1. Física estatística. 2. Monte Carlo. 3. Histograma
múltiplo. 4. Comportamento crítico. 5. Filmes finos
magnéticos I. Título.
530.13 CDD (22.ed.) FQ2008-036
À minha querida filha Ennya.
Agradecimentos
Sou grata ao pessoal do Laboratório de Computação Científica, André, Aranildo, Cesar
Sampaio e Paulo Gustavo. Obrigada pelas participações efetivas, contribuições e sugestões,
durante nossos seminários de grupo. Agradeço o apoio e amizade de todos.
Agradeço ao meu orientador Brady, pela grande amizade. Paciente, sempre receptivo e
disposto a responder qualquer tipo de pergunta. É um privilégio ter tido a oportunidade de
ser orientanda de uma pessoa que sabe ouvir e ensinar desta maneira singular. Que consegue
cuidar da formação básica do estudante e, ao mesmo tempo, incentivá-lo à adquirir autonomia.
Fui muito beneficiada por suas perguntas e observações pertinentes nos seminários, pela sua
sinceridade nas críticas e elogios. Pelas cobranças e pressões (discretas) na fase de escrita da
dissertação. Pela paciência durante as criteriosas correções de cada um dos capítulos desta
dissertação, sempre preocupado com a clareza. Pelo aprendizado...
Sou muita grata ao meu Coorientador, Adauto, com quem aprendimuito desde a iniciação
científica. Pelos vários puxões de orelha (coisas que agentesó agradece depois que vê o re-
sultado... É quando passamos a entender que é necessário). Durante todo este tempo, foram
muitas lições de coragem para enfrentar adversidades. Obrigada pelos conselhos acadêmicos
e pelas longas discussões sobre este trabalho. Agradeço muito pelas grandes contribuições à
minha formação básica. Pelo apoio e amizade.
Sou uma pessoa de sorte, pelo privilégio de ser orientada porBrady e Adauto.
Agradeço ao professor Jairo Rolim, pela paciência e clareza nas explicações durante os
seminários do nosso grupo.
Agradeço ao Cesar Sampaio, pelo companheirismo durante o Mestrado e pelas gentilezas
v
AGRADECIMENTOS vi
que me fez na fase de escrita desta dissertação.
Obrigada Paulo Gustavo, pelas boas conversas e preocupaçãocom o meu bem estar e segu-
rança.
Meus agradecimentos ao Aranildo, pela manutenção de software e ao André pela manu-
tenção de hardware dos recursos do Laboratório de ComputaçãoCientífica. Agradeço também
pela torcida e interesse que eles mostraram a respeito do andamento deste trabalho.
Obrigada meus amigos e companheiros de jornada, Manoel Vieira, Marco Polo, Gilvânia
Lúcia e Rebeca Cabral. Obrigada pelas longas discussões sobreas “listas” e elas foram impor-
tantes para minha formação. Agradeço à Rebeca Cabral pela amizade.
Sou grata ao amigo Francisco Vieira, pelo apoio moral em um momento tão crucial.
Agradeço aos professores Alexandre Ricalde e Lúcio Acioli, pelas contribuições à minha
formação básica complementar na graduação, sem as quais me seria impossível ingressar e
concluir o Mestrado.
Aos professores Giovani Vasconcelos, Rita Zorzenon e SandraSampaio, agradeço pela
formação básica no Mestrado.
Sou grata ao professor Jairo Rocha, da UFRPE, pelo apoio com os livros-texto.
Agradeço à secretária da pós-graduação, Sara Fonsêca, pelaeficiência e boa vontade sempre
que precisei resolver assuntos burocráticos, e pelas boas conversas.
Ana Paula e Flávia me ajudaram muito durante várias fases desta jornada e eu sou imensa-
mente grata às duas.
Agradeço ao meu marido Gilmar, pela paciência e por ter sido um pai totalmente presente
e ter cuidado da nossa filha durante as minhas ausências.
Sou muito grata aos meus irmãos Aninha, Luizinho e Lili, peloincentivo e estima.
Finalmente, os meus profundos agradecimentos aos meus pais, Luiz Carlos e Maria Cé-
lia, pelo incentivo para prosseguir nos estudos, pelo apoioem vários aspectos da vida, pela
credibilidade e pelo carinho.
somewhere i have never travelled,gladly beyond any experience,your eyes
have their silence: in your most frail gesture are things which enclose me,
or which i cannot touch because they are too near
your slightest look easily will unclose me though i have closed myself as
fingers, you open always petal by petal myself as Spring opens (touching
skillfully,mysteriously)her first rose
or if your wish be to close me,i and my life will shut very
beautifully,suddenly, as when the heart of this flower imagines the snow
carefully everywhere descending;
nothing which we are to perceive in this world equals the power of your
intense fragility:whose texture compels me with the colour ofits countries,
rendering death and forever with each breathing
(i do not know what is is about you that closes and opens; only something
in me understands the voice of your eyes is deeper than all roses)
nobody,not even the rain,has such small hands
—E. E. CUMMINGS (somewhere i have never travelled,gladly beyond,
1931)
Resumo
Usamos simulações Monte Carlo e teoria de escala de tamanho finito para estudar o com-
portamento magnético de filmes finos, em função da espessuraLZ dos filmes. Consideramos
sistemas de spinS= 1, com interações entre primeiros vizinhos em estruturas cúbicas simples
na geometria de filmes,L×L×LZ, com condições de contorno periódicas ao longo do plano
do filme e abertas na direção perpendicular ao filme. Obtivemos a energia, a magnetização, a
susceptibilidade e o calor específico em função da temperatura, para diversos valores de L e
espessuraLZ no intervalo[1 : 20]. Para estudar o comportamento do sistema na região crítica,
são geradas longas séries de dados para um conjunto de valores da temperatura e, através da
técnica de histogramas múltiplos, são obtidas as grandezasrelevantes em uma faixa contíua
de temperaturas. Em particular obtivemos a dependência da temperatura críticaTC(LZ) e dos
expoentes críticos, com a espessura do filme. Concluímos queTC(LZ) se aproxima do valor da
temperatura crítica do sistema volumétrico(3D), segundo uma lei de potência com expoente
λ = 1.40(3). Quanto aos expoentes críticosβ , γ e ν não observamos o crossover de 2D para
3D, e os valores obtidos são compatíveis com os expoentes exatos do sistema bidimensional
(2D).
Palavras-chave: Monte Carlo, histograma múltiplo, comportamento crítico, filmes finos mag-
néticos.
viii
Abstract
We use Monte Carlo simulations and finite size scaling theory to study the magnetic properties
of Ising thin films, as a function of temperature and film thicknesses. We consider films with
S= 1 and nearest-neighbor ferromagnetic interactions on simple cubic structures of sizeL×
L×LZ. We calculate the magnetization, the magnetic internal energy, susceptibility and specific
heat as a function of temperature, for several layer sizesL and film thicknessesLZ in the interval
[1 : 20]. In order to investigate the critical behavior of these systems we perform Monte Carlo
simulations at some values of temperature within the critical region, and extrapolate the results
over a wide range of temperatures through the multiple histogram technique. In particular, we
obtain how the critical temperatureTC(LZ) and critical exponents depend on thickness. We
conclude thatTC(LZ) varies from 2D to 3D values with increasing film thickness, and it shows
a convergence to the bulk critical temperature according toa power law with shift exponent
λ = 1.40(3). However, we do not observe the dimensional crossover of thecritical exponents
β , γ andν ; our results are consistent with the exact 2D exponents.
Keywords: Monte Carlo simulation, multiple histogram, critical behavior, magnetic thin
films.
ix
Sumário
1 Introdução 1
2 O Método 8
2.1 O Método de Monte Carlo 8
2.1.1 Conceitos fundamentais e simulações Monte Carlo 8
2.1.2 O Ensemble Canônico 11
2.1.3 O Algoritmo de Metropolis 12
2.2 O Método do Histograma 16
2.2.1 O método do histograma simples 16
2.2.2 O método do histograma múltiplo 20
3 Resultados e Discussões 30
3.1 Procedimentos computacionais 30
3.1.1 As grandezas termodinâmicas 33
3.1.2 As funções respostas 34
3.2 Comportamento termodinâmico 35
3.3 Comportamento crítico 40
3.3.1 Escalamento de tamanho finito 40
3.3.2 O cálculo da temperatura crítica e dos expoentes críticos 43
4 Conclusões 56
x
Lista de Figuras
2.1 Energia como uma função do tempo, para um filme com extensão L = 100
e espessuraLZ = 10. As simulações foram realizadas nas temperaturasT =
3.025, 3.069 e 3.15, sendoTC = 3.069 a temperatura crítica deste sistema. O
valor deτ∗ foi estimado com base na relaxação obtida emT = TC. 15
2.2 Magnetização como uma função do tempo, para o mesmo sistema apresentado
na figura (2.1) e simulações realizadas nas mesmas temperaturas. O tempo
característico,τ∗, estimado neste caso é obtido com base na relaxação da mag-
netização emT < TC. 15
2.3 Energia em função da temperatura para um filme com espessura LZ = 5 e ta-
manhoL = 80. Os pontos são resultados de simulações Monte Carlo usuais.
A linha contínua é a realização da extrapolação do resultadoda simulação em
T0 = 2.885, representado pelo símbolo cheio, via método do histograma sim-
ples. A linha tracejada, que resulta da extrapolação dos resultados obtidos de
todas as simulações, foi obtida através do método do histograma múltiplo. 19
2.4 Histogramas simples resultantes de simulações realizadas nas temperaturasT =
2.80, 2.84 e 2.88 para um filme de tamanhoL = 50 com espessuraLZ = 5. 21
2.5 Histograma múltiplo resultante da combinação dos resultados de simulações
para um filme comL = 50 eLZ = 5, em várias temperaturas correspondentes
aos histogramas simples especificados na figura. 23
xi
LISTA DE FIGURAS xii
2.6 Histogramas múltiplos gerados a partir da combinação dos resultados de si-
mulações em duas e três temperaturas, para um filme com espessuraLZ = 4 e
estensãoL = 50. 24
3.1 Filmes magnéticos formados por planos atômicos, de ladoL e espessuraLZ. O
valor deLZ indica o número de planos atômicos da camada magnética. 30
3.2 Estrutura cúbica simples em geometria de filme composta por duas camadas
atômicas. Os pontos representam as partículas de spinS= 1 que residem nos
sítios da rede. 31
3.3 Energia média em função da temperatura para vários tamanhos de filmes com
espessurasLZ = 5 eLZ = 18. 36
3.4 Calor específico em função da temperatura para vários tamanhosL de filmes
com espessurasLZ = 5 eLZ = 18. 37
3.5 Magnetização como uma função contínua deT para vários tamanhos de filmes
com espessurasLZ = 5 eLZ = 18. 39
3.6 Susceptibilidade magnéticaχ em função da temperatura para vários valores de
L para filmes com espessurasLZ = 5 eLZ = 18. 39
3.7 Cumulante de BinderuL em função da temperatura comLZ = 5 e vários valores
deL. O cruzamento das curvas resulta emTC = 2.881(3). Em (b) mostramos
a região crítica em uma escala menor. O erro apresentado na estimativa da
temperatura crítica indica que o valor deTC é a média aritmética dos valores de
temperaturas correspondentes ao primeiro e último cruzamentos. 45
3.8 Cumulante de BinderuL em função deT para vários valores deL e espessuras
do filmeLZ = 18. O cruzamento das curvas forneceTC = 3.144(3). 45
3.9 Temperatura crítica como uma função da espessura dos filmes. Indicamos os
valores deTC em duas e três dimensões obtidos neste trabalho. Estes valores
concordam com aqueles publicados nas referências [42] e [43], respectivamente. 46
LISTA DE FIGURAS xiii
3.10 Deslocamento da temperatura crítica do filme em direçãoao valor tridimensio-
nal. A inclinação da reta dá o expoente de deslocamentoλ . 47
3.11 Derivada do cumulante de Binder, emT = TC, como uma função deL para vá-
rias espessurasLZ. As inclinações das retas fornecem os valores de 1/ν mos-
trados na Tabela 3.1. 48
3.12 Derivada do cumulante de Binder, emT = T∗u (L), como uma função deL para
várias espessurasLZ. As inclinações das retas fornecem outra estimativa para
o expoente 1/ν (Tabela 3.1). 49
3.13 Magnetização em função deL, calculado na temperatura críticaTC. As inclina-
ções das retas fornecem os valores deβ/ν mostrados na Tabela 3.2. 51
3.14 Derivada da magnetização em função deL, calculado emTC. As inclinações
das retas fornecem os valores de(1−β )/ν mostrados na Tabela 3.2. 52
3.15 Susceptibilidade em função deL, calculado na temperatura críticaTC. As in-
clinações das retas fornecem os valores deγ/ν mostrados na Tabela 3.2. 52
3.16 Derivada da susceptibilidade em função deL, calculado emTC. As inclinações
das retas fornecem os valores de(1+ γ)/ν mostrados na Tabela 3.2. 53
3.17 Calor específico em função do logaritmo deL, calculado na temperatura pseudo-
críticaT∗c (L) e na temperatura críticaTC. Estes resultados indicam que a diver-
gência do calor específico é logarítimica, o que implica emα = 0. 54
Lista de Tabelas
3.1 Temperatura críticaTC e expoente 1/ν calculados a partir do cumulante de
Binder. Os valores de 1/ν foram obtidos dos ajustes lineares (a) da Eq. (3.27)
e em (b) da Eq. (3.29). Na última coluna mostramos o valor do cumulante de
Binder emT = TC. Em destaque, os valores conhecidos dessas grandezas para
duas (2D) [39, 47, 44] e três (3D) dimensões [37]. 50
3.2 Expoentesβ/ν e γ/ν obtidos dos ajustes lineares (a) das expressões (3.18) e
(3.20) e em (b) das suas respectivas derivadas (3.21) e (3.22) emT = TC. Em
destaque, os valores conhecidos desses expoentes para duas(2D) e três (3D)
dimensões [39, 46]. 50
3.3 Temperatura críticaTC e expoente crítico 1/ν obtidos da relação de escala
(3.28), na temperatura pseudo-críticaT∗χ associada à susceptibilidade. 55
xiv
CAPÍTULO 1
Introdução
O estudo de sistemas magnéticos em geometria de filme é importante tanto do ponto de
vista da ciência básica quanto das aplicações tecnológicas. Do ponto de vista das aplicações,
tais sistemas são peças fundamentais na construção de sensores [1, 2], dispositivos de leitura
e gravação em discos rígidos de computadores [3] e, mais recentemente, memória magnética
de acesso randômico [4]. Além dessas e outras aplicações tecnológicas, os filmes finos são
interessantes porque, com eles, se pode explorar e testar concepções fundamentais como a
hipótese de universalidade em fenômenos críticos [5, 6]. Defato, estudos de filmes magnéticos
ultrafinos revelaram uma infinidade de novos fenômenos que seriam inimagináveis em sistemas
estritamente bidimensionais ou muito menos em estruturas tridimensionais [7].
Esses fatos despertam o interesse de vários pesquisadores eum grande esforço tem sido
despendido para a caracterização do comportamento das propriedades magnéticas de heteroes-
truturas formadas pelo empilhamento de diferentes materiais magnéticos. E, graças ao desen-
volvimento das técnicas de fabricação de materiais, houve um grande impulso na investigação
das propriedades físicas de filmes finos nos anos recentes. Emparticular, investiga-se como
as propriedades magnéticas variam com a espessura das camadas. Um caso especialmente in-
teressante ocorre em filmes finos magnéticos homogêneos paraos quais a dimensionalidade
espacial não é bem estabelecida. Neste caso, observa-se umamudança de um comportamento
bidimensional para um comportamento tridimensional. A geometria de filme é caracterizada
pelo fato de uma das três dimensões do sistema ser muito menorque as demais.
Experimentalmente observa-se que o comportamento dos filmes muda quando a sua espes-
sura é variada. A temperatura crítica e os expoentes críticos de filmes finos ferromagnéticos
1
CAPÍTULO 1 INTRODUÇÃO 2
têm sido exaustivamente investigados experimentalmente [8, 9, 10]. É bem conhecido que o
valor da temperatura de Curie de filmes magnéticos muda de um valor tipicamente bidimensi-
onal para um valor correspondente a um sistema tridimensional à medida que a espessura dos
filmes aumenta. Todavia a classe de universalidade de um filmemagnético deve corresponder
a duas dimensões [11, 12]. O deslocamento da temperatura crítica do filme em direção ao valor
tridimensional à medida que a espessura cresce é descrita pelo expoente de deslocamentoλ .
Argumentos de grupo de renormalização sugerem queλ = 1/ν , ondeν é o expoente crítico
em três dimensões associado à divergência do comprimento decorrelação [11, 13, 12].
Um filme magnético formado por vários planos atômicos é uma estrutura tridimensional.
No limite de um único plano atômico o sistema é essencialmente bidimensional eTC e os cor-
respondentes expoentes críticos assumem valores bidimensionais. À medida que a espessura
do filme aumenta e a temperatura é mantida bem longe deTC, o comprimento de correlação é
muito menor que a espessura do filme e o seu comportamento termodinâmico é tipicamente tri-
dimensional. Para temperaturas suficientemente próximas deTC o comprimento de correlação,
ξ , tornar-se maior do que a espessura do filme e seu comportamento muda de três para duas
dimensões. Ou seja, um cruzamento dimensional ocorre quando o comprimento de correlação
começa a divergir e ultrapassa a espessura do filme em alguma temperatura próxima deTC.
Resultados de simulações Monte Carlo e análise de campo médio para filmes finos mode-
lados por sistemas Ising spin 1/2, comprovam o cruzamento dimensional para a temperatura
crítica,TC, bem como o caráter bidimensional dos expoentes críticos [14].
Não obstante, foi observado experimentalmente uma drástica mudança no valor do expoente
crítico associado à magnetização em filmes finos de níquel [9]. Os autores creditam tal variação
do expoenteβ a um cruzamento de um caráter bidimensional para tridimensional do filme
quando a espessura atinge entre 5 e 7 camadas atômicas. Já as propriedades críticas dos filmes
mais finos mostraram um caráter basicamente bidimensional.Uma possível explicação para
tal fenômeno pode ser atribuida ao cruzamento de um caráter Ising para um comportamento
CAPÍTULO 1 INTRODUÇÃO 3
Heisenberg, induzido pela geometria (anisotropia de forma). Ou seja, o comportamento crítico
de filmes ultrafinos é bidimensional, enquanto que os filmes espessos podem apresentar outras
características críticas. Dito de outra forma, os graus de liberdade magnéticos de filmes finos
podem ser bem descritos por spins de Ising enquanto que os de filmes mais espessos por spins
de outra simetria. Um estudo conduzido na referência [14] emfilmes magnéticos de Ising
com spin 1/2 em estruturas cúbicas revelou que o caráter crítico bidimensional é mantido,
mesmo para filmes com até 20 camadas atômicas, apesar de um leve desvio dos valores típicos
bidimensionais ter sido observado para os filmes mais espessos.
Em contraste, mais recentemente [15], simulações Monte Carlo similares às apresentadas
na referência [14] detectaram uma mudança no comportamentodos expoentes críticos com a
espessura de filmes, descrito por um modelo de Ising(S= 1/2) com interações entre primeiros
vizinhos. Os autores calcularam a magnetização e a susceptibilidade magnética associadas
tanto às superfícies dos filmes quanto ao seu interior, separadamente. Os expontes extraídos
com as propriedades apenas de superfície, apenas do volume (interior) ou total, concordaram
entre si e todos indicam um cruzamento de duas para três dimensões. Em suas análises, eles
definiram uma dimensão efetiva,de f, através da relação de hiperescala,de f = γ/ν + 2β/ν e
concluíram quede f aumenta monotonicamente de um valor aproximadamente 2, para filmes
formados por poucas monocamadas atômicas, em direção à 3 quando a espessura cresce.
Muitas das investigações teóricas sobre filmes magnéticos de Ising foram realizadas para
S= 1/2, onde os spins assumem apenas os valoresσ = ±1. Porém, o caso de spins de Ising
comS= 1, onde os spins assumem os valoresσ = ±1 ou zero, é pouco estudado no limite em
que a anisotropia de íon único é nula. Não obstante, esse casolimite fornece um excelente teste
para a hipótese de universalidade para os expoentes críticos. Além disso, o estadoσ = 0 pode
ser interpretado como a possibilidade do parâmetro de ordemapontar numa direção perpen-
dicular ao eixo de fácil magnetização. Estudamos nesta dissertação o comportamento crítico
de sistemas magnéticos com simetria de filme. Os filmes são modelados por spins de Ising,
CAPÍTULO 1 INTRODUÇÃO 4
que podem assumir os estados−1, 0 e+1, e que interagem apenas com os seus vizinhos mais
próximos. Os filmes são modelados por estruturas cúbicas simples, com condições de con-
torno apropriadas como descrito no Capítulo 3. Estudamos este sistema através de simulações
Monte Carlo, onde empregamos o algoritmo de Metropolis [16].Detalhes do método Monte
Carlo (MMC) que empregamos será apresentado no Capítulo 2. Aquifaremos uma breve re-
visão de sua aplicação à mecânica estatística e discutiremos brevemente alguns dos algoritmos
empregados em simulações de modelos de spins na rede.
O método Monte Carlo é uma poderosa e versátil técnica numérica para resolver problemas
matemáticos. Basicamente, o método consiste em simular o comportameto de variáveis aleató-
rias. Embora seja uma técnica estocástica, o método é aplicável tanto à problemas de natureza
intrinsicamente estocástica quanto à problemas determinísticos.
A invenção do MMC é atribuída ao matemático Stanislaw Ulam que, durante um repouso
forçado em 1946, iniciou uma reflexão sobre a probabilidade de alguém vencer o jogo de cartas
conhecido como paciência ou solitário. Depois de algumas tentativas de resolver o problema
analiticamente através de técnicas combinatórias, ele percebeu que uma boa solução poderia
ser obtida por um meio mais simples e intuitivo. Ulam sugere jogar algo como uma centena de
partidas e contar o número de partidas ganhas. A razão entre onúmero de partidas vitoriosas e o
número de jogadas é uma boa aproximação para a solução do problema proposto. Na verdade
Ulam não inventou o método Monte Carlo, a base matemática do método já era muito bem
estabelicida bem antes de 1946. Porém, ele foi o primeiro a reconhecer a potencialidade das
novas máquinas computadoras na aplicação da técnica em problemas realísticos.
Simulações Monte Carlo tem sido empregadas com sucesso no estudo de sistemas mode-
los da Mecânica Estatística tanto nas suas versões clássicas quanto as correspondentes versões
quânticas [17, 18, 19]. A área de transições de fases e fenômenos críticos alcançou apreciável
progresso com a disseminação do método e o aumento do poder computacional dos computa-
dores digitais.
CAPÍTULO 1 INTRODUÇÃO 5
Outra razão para o sucesso da técnica em descrever com precisão o comportamento de
sistemas físicos complexos, deve-se ao desenvolvimento dealgoritmos rápidos e eficientes que
reduzem efeitos indesejáveis intrínsecos ao método. Por exemplo, o conhecido fenômeno de
amortecimento crítico que ocorre próximo a uma transição defases, devido às fortes correlações
temporarais, entre as configurações geradas, induzidas pelas grandes flutuações que geralmente
estão presentes nesses casos [20].
Originalmente, a formulação do método MC foi no ensemble canônico [16] e, desde então,
um enorme número de estudos foram conduzidos nos mais diversos sistemas empregando a
prescrição de Metropolis. No contexto de uma simulação MC à la Metropolis, o parâmetro de
controle é a temperatura e o sistema é considerado estar em contato com um reservatório de
calor à temperatura especificada. Além disso, como será detalhado no Capítulo 2, a dinâmica
de Metropolis consiste em tentar atualizar o estado de uma única partícula de cada vez. Sob
certas circunstâncias, esta dinâmica pode se tornar bastante ineficiente, uma vez que o sistema
pode ficar armadilhado em uma configuração para certos valores do parâmetro de controle.
Problemas dessa natureza motivaram a invenção de novos algoritmos. O amortecimento crí-
tico foi praticamente superado com a introdução dos algoritmos de ilhas por Swendsen e Wang
[21]. A idéia básica subjacente à proposta consiste em identificar a escala das correlações espa-
ciais e atualizar coletivamente todas as partículas dentrode uma certa região. Especificamente,
aglomerados ou ilhas de partículas são identificadas introduzindo-se uma ligação entre partí-
culas interagentes que se encontram no mesmo estado, com umaprobabilidade dependente da
temperatura. As partículas que não estão no mesmo estado sãoconsideradas desacopladas.
Um aglomerado é formado pelo conjunto de partículas conectadas entre si. Note que todas as
partículas em um dado aglomerado estão no mesmo estado, por construção. Note ainda que
nem sempre duas partículas interagentes pertencem ao mesmoaglomerado, devido à maneira
estocástica de sua construção. Assim, o sistema de partículas é decomposto em um conjunto
de ilhas desacopladas entre si. Após as ilhas terem sido formadas, é escolhido um novo estado
CAPÍTULO 1 INTRODUÇÃO 6
aleatório para cada uma delas. Todas as partículas de uma determinada ilha vão para o estado
selecionado para ela. Observe que, como o novo estado é escolhido com igual probabilidade
entre os possíveis estados de uma partícula, há uma probabilidade finita de qualquer ilha per-
manecer no mesmo estado. Maiores detalhes podem ser encontrados nas referências [18] e
[17].
Um outro algoritmo de ilhas foi proposto por Wolff [22] e baseia-se no algoritmo de Swend-
sen e Wang. Neste, as ilhas são estocasticamente formadas e depois atualizadas também proba-
bilisticamente. Assim, ilhas de todos os tamanhos (número de partículas) podem ser formadas
e qualquer delas será atualizada ou não com a mesma probabilidade. Um procedimento com-
putacionalmente mais eficiente seria, por exemplo, atualizar apenas a maior ilha. Isto garantiria
uma maior diferença entre a nova e a velha configurações. Na verdade, não precisamos cons-
truir todas as ilhas. De fato, se escolhermos uma partícula do sistema ao acaso a probabilidade
dela pertencer a uma dada ilha é proporcional ao tamanho da ilha. Wolff então propõe sortear
uma partícula qualquer do sistema e considerá-la a semente da ilha a ser atualizada (muito pro-
vavelmente, a maior das possíveis ilhas). A semente escolhida é imediatamente atualizada e
novas partículas são incorporadas ao aglomerado de forma similar ao procedimento de Swend-
sem e Wang. Cada nova partícula acrescentada à ilha é atualizada e tornar-se uma semente
também. A conexão entre uma semente e uma partícula não incorporada é testada uma única
vez. A ilha para de crescer quando deixar de existir conexõesque ainda não foram testadas.
Na descrição dos algoritmos anteriores, assumimos tacitamente que os possíveis estados das
partículas são discretos. A generalização para o caso contínuo pode ser encontrado no artigo
original do Wolff [22].
Uma outra classe importante de algoritmos Monte Carlo consiste na tentativa de se estimar
diretamente a densidade de estados do sistema de interesse.O broad histogram method[23]
é uma das primeiras tentativas bem sucedida de um tal algoritmo. Um avanço em relação ao
broad histogram methodfoi introduzido por Wang e Landau [24]. Por fim, citamos o método
CAPÍTULO 1 INTRODUÇÃO 7
de Hüller e Pleimling que consiste em estimar diretamente a entropia microcanônica [25], que
obteve algum sucesso em descrever o comportamento de sistemas discretos. Uma generalização
para o caso contínuo pode ser encontrado na referência [26]
Nessa dissertação escolhemos empregar o método Monte Carlo segundo o algoritmo de
Metropolis. Dado a sua simplicidade de implementação, ele constitui-se no procedimento mais
geral e menos susceptível a erros de implementação.
Qualquer método Monte Carlo fornece um meio de se calcular valores médios dos obser-
váveis físicos apenas para sistemas de tamanhos finitos. Um meio de se extrair informações no
limite termodinâmico é através da teoria de escalamento como tamanho finito [27, 28, 29, 30]
que será tratada na secção 3 do Capítulo 3.
Por último, mais não menos importante, chegamos ao ponto de como analisar os dados
obtidos numa simulação Monte Carlo. De um certo ponto de vista, os dados gerados numa
simulação são analisados mais ou menos como um físico experimental analisaria os dados
obtidos em seu laboratório. Ou seja, calcula-se os valores médios e tenta-se estimar os erros
estatísticos e os erros sistemáticos envolvidos no procedimento de medição. Porém, no caso
de uma simulação Monte Carlo existem alguns procedimentos que são peculiares ao método.
No Capítulo 2, apresentamos o método do histograma múltiplo [31] que permite extrapolar os
resultados de simulações realizadas em um conjunto de temperaturas, para todo um contínuo
de temperaturas em uma ampla faixa de valores.
CAPÍTULO 2
O Método
2.1 O Método de Monte Carlo
2.1.1 Conceitos fundamentais e simulações Monte Carlo
O Método de Monte Carlo (MMC) é um método de solução numérica quesimula o com-
portamento de variáveis aleatórias, e consiste essencialmente em gerar números aleatórios.
Embora introduzido [32] por Nicolas Metropolis e StanislawUlam em 1949, o princípio deste
método é conhecido há muito tempo e, não fosse o volume considerável de cálculo associado
à simulação de variáveis aleatórias, ele teria sido amplamente aplicado a modelos estatísticos
como, por exemplo, no tratamento de amostras aleatórias. O surgimento dos computadores ele-
trônicos tornou possível as aplicações do método de MC a uma grande variedade de problemas
e permitiu a sua ampla e rápida difusão [33].
Com base no livro do Sobol [33] apresentaremos nesta seção umabreve revisão dos con-
ceitos estatísticos básicos necessários para compreensãodo método de Monte Carlo. Para en-
tender como o MMC funciona, é primordial conhecer o significado de uma variável aleatória.
Matemáticamente, variável aleatória é uma variável que pode assumir diversos valores, com
uma determinada probabilidade. Se os possíveis valores queesta variável pode assumir com-
põem um conjunto inteiro de valores, então dizemos que esta éuma variável aleatória discreta.
Variáveis aleatórias que podem assumir qualquer valor num intervalo real, são ditas variáveis
aleatórias contínuas. Uma variável aleatória contínuaξ é caracterizada por um intervalo(a,b)
compreendendo o seu domínio e por uma funçãop(x) densidade de probabilidade. Neste caso,
a probabilidade deξ assumir um valor no intervalo(a′,b′), pertencente a(a,b) é igual a
8
2.1 O MÉTODO DE MONTE CARLO 9
P{a′ < ξ < b′} =∫ b′
a′p(x)dx.
Embora conhecida a probabilidade de ocorrência de cada valor, não é possível saber exa-
tamente qual deles a variável aleatória assumirá numa determinada amostragem. Contudo, o
caráter de seu comportamento pode ser previsto uma vez realizada uma série de experiências,
e esta informação é tão mais precisa quanto maior for a série de amostragens consideradas.
O valor esperado e a variância constituem as características numéricas básicas das variáveis
aleatórias. Considere uma variável aleatóriaξ que pode assumir com uma dada probabilidade,
uma sequência longa porém finita de valoresx1,x2, ...,xN. A média aritmética destes valores
será próxima do valor esperado,〈ξ 〉, isto é,(x1 +x2 + ...+xN)/N ≈ 〈ξ 〉. O grau de dispersão
dos valores assumidos porξ em torno de〈ξ 〉 é caracterizado pela variância,δ 2ξ , que se rela-
ciona com o desvio padrão porσ =√
(δ 2ξ ). A quantidadeσ caracteriza o erro do método de
medição adotado.
De acordo com a regra dos 3 sigmas, a probabilidade deξ assumir um valor entrem−3σ
em+3σ , é
P{m−3σ < ξ < m+3σ} =∫ m+3σ
m−3σp(x)dx= 0,997, (2.1)
ondem é o valor esperado. A equação (2.1) implica que, numa experiência, é muito provável
que o desvio entreξ e 〈ξ 〉 seja inferior a 3σ .
SejamN variáveis aleatórias independentes,ξ1,ξ2, ...,ξN, igualmente distribuídas. O fato
de tais variáveis obedecerem à mesma distribuição de probabilidades, implica em dizer que
estas possuem as mesmas esperanças e variâncias. Ou seja,〈ξ 〉i = meδ 2ξi = b2, para qualquer
que sejai.
No que diz respeito à independência de variáveis aleatórias, pode-se dizer que, duas variá-
veisξ eη são independentes se a distribuição deξ não depender do valor que terá assumidoη .
2.1 O MÉTODO DE MONTE CARLO 10
Logo, o valor esperado do produto entre elas será〈ξ η〉 = 〈ξ 〉〈η〉, e a variância da soma será
δ 2(ξ +η) = δ 2ξ +δ 2η .
Se considerarmos uma variável,ρN, tal queρN = ξ1+ξ2+ ...+ξN é a soma dasN variáveis,
então o valor esperado e a variância deρN serão respectivamente,〈ρ〉 = Nme δ 2ρ = Nb2.
Suponha finalmente uma variável aleatória,η , normalmente distribuída, e que seja carac-
terizada pelos parâmetrosa = Nme σ2 = Nb2. O teorema central do limite garante que, para
qualquer intervalo(a′,b′) eN suficientemente grande,
P{a′ < ρN < b′} ≈∫ b′
a′pηN(x)dx, (2.2)
o que implica em afirmar que a soma de um número elevado de variáveis aleatórias indepen-
dentes é aproximadamente normal, ou seja,(pρN(x) ≈ pηN(x)). Este teorema vale ainda para
variáveis aleatórias que não são independentes, nem normalmente distribuídas.
É de se esperar que, dado o grande número de fatores que influenciam o resultado de uma
amostra, o comportamento individual de cada fator não deve ser relevante. O elevado número
de termos da soma explica porque cada um deles não precisa, necessariamente, obedecer a uma
distribuição normal de probabilidades.
Quanto à validade do teorema para o caso de variáveis que não são independentes, suponha
que em uma determinada prova realizada paraN variáveis aleatórias, algumas delas apresen-
taram uma distribuição dependente do valor que alguma outraassumiu. Ou seja, existe uma
correlação entre os resultados de tais variáveis. Neste caso é possível selecionar apenas os da-
dos descorrelacionados. Denotemos porn, a extensão da correlação. Para cada conjunto de
dados correlacionados, apenas um elemento do grupo será considerado. Ou seja, ao invés deN
termos, a soma terá agoraN/n termos relevantes. Sen é finito eN grande o suficiente, então
a razãoN/n é ainda um número grande. Assim, a soma deN/n termos é aproximadamente
normal, o que valida o teorema central do limite.
No contexto de uma simulação Monte Carlo, a regra dos três sigmas e a aplicação do te-
2.1 O MÉTODO DE MONTE CARLO 11
orema central do limite são fundamentais. Suponha que, interessados em avaliar uma dada
grandezam, imaginemos uma variável aleatóriaξ , tal que〈ξ 〉 = me δ 2ξ = b2.
Para um grande número,N, de variáveis aleatórias independentes distribuídas igualmente a
ξ , temos que, de acordo com o teorema central do limite, a distribuição da soma,ρN, de tais
variáveis será aproximadamente normal, com os parâmetrosa = Nmeσ2 = Nb2. Da regra dos
três sigmas, equação (2.1), podemos escrever que
P
{
m− 3b√N
<ρNN
< m+3b√N
}
≈ 0,997 (2.3)
Lembrando queρN = ∑Ni ξi e usando as propriedades do módulo das funções de variáveis
reais, podemos reescrever (2.3) como
P
{∣
∣
∣
∣
∣
1N
N
∑i
ξi −m∣
∣
∣
∣
∣
<3b√N
}
≈ 0,997 (2.4)
A relação (2.4) mostra que, com alta probabilidade, o valor estimado param é muito pró-
ximo da média aritmética tomada sob asN variáveis aleatórias, e o erro cometido na estimativa
será no máximo igual a 3b/√
N. Exatamente emm= (∑Ni ξi)/N, temos o menor valor possível
para 3b/√
N. Quanto maior o valor deN, melhor será a estimativa dem e, consequentemente,
menor será o erro.
2.1.2 O Ensemble Canônico
Em Mecânica Estatística, a obtenção da função de partição deum sistema permite a cone-
xão com a termodinâmica e a obtenção dos valores médios dos observáveis de interesse. No
ensemble canônico o valor médio de um dado observável físicoé calculado para um determi-
nado valor fixo de temperatura, ou seja,
2.1 O MÉTODO DE MONTE CARLO 12
〈A〉β = ∑s
Pβ (s)A(s), (2.5)
onde
Pβ (s) =1
Zβe−βEs (2.6)
é a probabilidade de ocorrência de uma dada configuraçãos quando o sistema encontra-se
termalizado a uma temperaturaT, comβ = 1/kBT, e
Zβ = ∑s
e−βEs (2.7)
é a função de partição canônica.
Determinar a função de partição analiticamente nem sempre épossível. De fato poucos são
os modelos que têm soluções exatas [34]. Por outro lado, o cálculo numérico da função de par-
tição através do MMC seria proibitivo devido ao grande número de configurações necessárias
para determinarZβ com uma precisão aceitável.
Vimos que o método de Monte Carlo (MMC) nos fornece uma boa estimativa do valor
médio dos observáveis de interesse, sem necessariamente conhecer a função de partição do
sistema. Esta estimativa é uma medida simples da média aritmética dos valores obtidos da
grandeza em cada uma das configurações de equilíbrio geradasde acordo com a distribuição
de probabilidade canônica. O valor médio estimado será tão próximo do correto, quanto maior
for o número de configurações independentes empregadas no cálculo.
2.1.3 O Algoritmo de Metropolis
A primeira aplicação do método Monte Carlo no contexto da mecânica estatística ocorreu
em 1953 e deve-se à Nicolas Metropolis [16]. O método foi desenvolvido com o propósito
de realizar simulações de um fluido no ensemble canônico. O algoritmo de Metropolis é hoje
2.1 O MÉTODO DE MONTE CARLO 13
intensivamente utilizado no estudo de sistemas de naturezas bastante diversificadas [35].
Computacionalmente o algoritmo de Metropolis consiste na seguinte regra. Dada a to-
pologia da rede, partimos de uma configuração de spinss0 qualquer e a cada intervalo de
tempo geramos uma nova configuração. Apósτ varreduras de rede teremos um conjunto
s1,s2, ...,sk, ..,sτ de configurações, ou microestados do sistema.
Um novo microestado é gerado da seguinte forma. Suponha que no k-ésimo instante de
tempo, o spin doi-ésimo sítio visitado sejaσi . Calculamos o custo energético associado à
mudança de estado do spinσi paraσ ′i , δE = Es′k −Esk. No ensemble canônico a distribuição
de probabilidades (2.6) de uma dada configuração,s, é proporcional ao fator de Boltzmann
exp(−βEs). A probabilidade de transição de um microestadosk paras′k, dada pela razão entre
Pβ (s′k) ePβ (sk) é
Wβ (sk → s′k) = e−βδE. (2.8)
SeδE for menor ou igual a zero, então a mudança é aceita, eσi → σ ′i . Por outro lado, se
δE > 0, o novo estadoσ ′i será aceito de acordo com a probabilidade de transiçãoW(σi → σ ′i ),
(2.8). Neste caso geramos um número aleatório,ξ , distribuído uniformemente no intervalo
[0,1]. Seξ ≤ W entãoσi → σ ′i . Caso contrário, o estado do spin da partícula não se altera.
Uma nova configuração é obtida(k → k+ 1) após este procedimento ser repetidoN vezes,
ondeN é o número de spins da rede. Isto garante que cada spin seja atualizado uma vez, em
média.
Para garantir que estamos realizando as medidas sobre configurações de equilíbrio a um de-
terminado valor deT constante, a rede deve ser percorrida um número suficientemente grande
de vezes, nacessário para a equilibração do sistema. O tempoque o sistema leva para atingir
o equilíbrio térmico é chamado de tempo característico,τ∗, (ou tempo de relaxação). Após
passado esse tempo, cada varredura de rede corresponde à geração de uma nova configuração
de equilíbrio. Observe que, nas simulações, o tempo é medidoem número de passos Monte
2.1 O MÉTODO DE MONTE CARLO 14
Carlo (MCS). Uma unidade de tempo equivale a um MCS, o que corresponde a uma varredura
de rede.
A qualidade do resultado obtido para os valores médios, estimados na simulação Monte
Carlo, é tão boa quanto maior for o número de configurações de equilíbrio geradas. Para reduzir
o erro cometido na estimativa, devemos também realizar amostras diferentes para um mesmo
valor deT. Neste contexto, uma boa estatística consiste em gerar não somente um número
elevado de configurações de equilíbrio, mas realizar váriasamostras para cada conjunto de
configurações geradas sob um mesmo valor deT.
Para verificar se o sistema atingiu o estado de equilíbrio, podemos observar a evolução tem-
poral de grandezas macroscópicas, tais como a energia configuracional e a magnetização. Na
verdade, a relaxação da magnetização nos dá uma melhor estimativa do tempo característico,
pois as flutuações são mais pronunciadas para a magnetizaçãodo que para energia. Este fato
é ilustrado nas figuras (2.1) e (2.2) onde temos, respectivamente, a relaxação da energia e da
magnetização a partir de uma configuração inicial desordenada, para um filme com as dimen-
sões do planoL = 100 e espessuraLZ = 10 através de simulações Monte Carlo. Se olharmos
apenas para a relaxação da energia, na figura (2.1) é possívelque o tempo característico seja
subestimado. Apóst = τ∗ a energia oscila pouco em torno de um valor bem definido (valor
médio de E), no entanto acima da temperatura crítica,T > TC, vemos que o sistema atinge o
equilíbrio para um tempo característico ainda menor. Isso ocorre porque a energia do sistema é
uma grandeza mais bem comportada. Podemos ver na figura (2.2)que as flutuações são ainda
mais pronunciadas paraT = TC.
O tempo computacional necessário para realizar este procedimento é proporcional ao nú-
mero de MCS, ao tamanho do sistema, ao número de temperaturas sob as quais o sistema foi
simulado e ao número de amostras realizadas para cada valor de temperatura. Portanto, boas
simulações Monte Carlo são computacionalmente caras.
2.1 O MÉTODO DE MONTE CARLO 15
0 2000 4000 6000 8000 10000tempo (MCS)
-1.25
-1
-0.75
-0.5
Ene
rgia
(J)
T > Tc
T = Tc
T < Tc
τ*
Figura 2.1 Energia como uma função do tempo, para um filme com extensãoL = 100 e espessuraLZ = 10. As simulações foram realizadas nas temperaturasT = 3.025, 3.069 e 3.15, sendoTC = 3.069 atemperatura crítica deste sistema. O valor deτ∗ foi estimado com base na relaxação obtida emT = TC.
0 2000 4000 6000 8000 10000tempo (MCS)
0
0.1
0.2
0.3
0.4
0.5
Mag
netiz
acao
T < Tc
T = Tc
T > Tc
τ*
Figura 2.2 Magnetização como uma função do tempo, para o mesmo sistema apresentado nafigura(2.1) e simulações realizadas nas mesmas temperaturas. O tempo característico, τ∗, estimado neste casoé obtido com base na relaxação da magnetização emT < TC.
2.2 O MÉTODO DO HISTOGRAMA 16
2.2 O Método do Histograma
A partir de simulações Monte Carlo realizadas em alguma temperatura específica é possí-
vel extrair informações sobre o comportamento das grandezas de interesse do sistema físico
estudado em temperaturas vizinhas, através de extrapolações numéricas de acordo com o mé-
todo do histograma. O método do histograma simples, originalmente proposto por Ferrenberg
e Swendsen em 1988 [31], tem sido aplicado à Mecânica Estatística no estudo de sistemas
magnéticos. Os autores usaram o método do histograma no estudo do modelo de Ising em três
dimensões e obtiveram com alta precisão a temperatura na qual o calor específico é máximo e
uma boa estimativa para o expoente crítico associado a esta função resposta [31].
2.2.1 O método do histograma simples
Apresentamos nesta secção a essência do método do histograma na versão mais simples e
discutiremos a sua aplicação no estudo de sistemas magnéticos.
Da definição de valor médio expressa em (2.5) podemos agruparem{s} cada conjunto de
configurações que correspondem a um dado valor de energiaE e, dessa forma expressar〈A〉βem termos de somas sobre os possíveis valores de energia, ou seja,
〈A〉β = ∑E
ĀEPβ (E), (2.9)
onde
ĀE =1
g(E)
(
∑s
A(s)δE,Es
)
(2.10)
é a média microcanônica e
Pβ (E) =g(E)e−βE
Zβ(2.11)
2.2 O MÉTODO DO HISTOGRAMA 17
é a probabilidade de o sistema exibir uma dada energiaE quando se encontra termalizado à
temperaturaT, β = 1/kBT, g(E) é o número de microestados com energiaE. Analogamente
podemos expressar a função de partição canônica, equação (2.7), na forma
Zβ = ∑E
g(E)e−βE. (2.12)
Seja ainda
Pβ0(E) =g(E)e−β0E
Zβ0, (2.13)
a probabilidade do sistema exibir a mesma energiaE quando termalizado à temperaturaT0,
correspondente aβ0 = 1/kBT0.
Das equações (2.11) e (2.13) podemos extrair informações dePβ (E) a partir do conheci-
mento dePβ0(E), ou seja
Pβ (E) = e−(β−β0)E
(
Zβ0Zβ
)
Pβ0(E). (2.14)
Note que basta conhecer a razão(
Zβ0/Zβ)
para quePβ (E) seja obtida em termos dePβ0(E).
Somando a expressão (2.14) sobre todas as possíveis energias, temos
(
Zβ0Zβ
)
=1
∑E e−(β−β0)EPβ0(E). (2.15)
Substituindo (2.15) em (2.14) podemos finalmente escrever
Pβ (E) =e−(β−β0)EPβ0(E)
∑E e−(β−β0)EPβ0(E). (2.16)
A relação acima é uma expressão analítica da Mecânica Estatística, válida independente-
mente de como a função distribuição canônicaPβ0(E) é obtida. No entanto, fora do contexto
de uma simulação computacional a expressão (2.16) torna-seinútil, dado a impossibilidade de
cálculo analítico dePβ0(E).
2.2 O MÉTODO DO HISTOGRAMA 18
Uma estimativa paraPβ0(E) pode ser obtida da seguinte forma. Suponha que realizamos
uma amostra comn configurações descorrelacionadas, geradas de acordo com a probabilidade
Pβ0(E). SejaN(E), o histograma da energiaE, ou seja o número de configurações que resulta-
ram no valor de energiaE. A lei dos grandes números [36] garante que
limn→∞
N(E)n
= Pβ0(E). (2.17)
Para um grande número de configurações, isto é, paran suficientemente grande
Pβ0(E) ≈N(E)
n. (2.18)
A última igualdade em (2.18) nos permite relacionar o histograma medido emβ0, e a dis-
tribuição de probabilidadePβ (E), como
Pβ (E) =N(E)e−(β−β0)E
∑E N(E)e−(β−β0)E. (2.19)
Usando a expressão (2.9), podemos escrever o valor médio de um dado observávelA como uma
função contínua deβ
〈A〉β = ∑E
ĀEN(E)e−(β−β0)E
∑E N(E)e−(β−β0)E. (2.20)
A expressão (2.20) é a relação fundamental do método do histograma simples. Através
dela podemos extrapolar o valor médio de qualquer grandeza física para temperaturas vizinhas
à temperaturaT0, na qual se realizou as simulações Monte Carlo. Contudo, a região em que
será feita a extrapolação não deve afastar-se muito deT0. Na verdade, um dos desafios deste
método é a incerteza sobre a faixa de temperatura na qual podemos realizar a extrapolação,
pois a exatidão dos valores extrapolados decresce quando nos aproximamos dos extremos do
histograma [17]. Ou seja, temperaturas vizinhas que estão mais próximas do valor no qual foi
2.2 O MÉTODO DO HISTOGRAMA 19
2,8 2,85 2,9 2,95Temperatura (J/k
B)
-1,2
-1
-0,8
Ene
rgia
(J)
Figura 2.3 Energia em função da temperatura para um filme com espessuraLZ = 5 e tamanhoL =80. Os pontos são resultados de simulações Monte Carlo usuais. A linha contínua é a realização daextrapolação do resultado da simulação emT0 = 2.885, representado pelo símbolo cheio, via método dohistograma simples. A linha tracejada, que resulta da extrapolação dos resultados obtidos de todas assimulações, foi obtida através do método do histograma múltiplo.
realizada a simulação apresentam os valores iguais àquelesque seriam obtidos em simulações
usuais, enquanto que os valores médios das grandezas que foram extrapolados na região dos
extremos do histograma não coincidem. A figura (2.3), para a energia média por partícula,
ilustra como podemos estimar a largura da faixa de extrapolação através da comparação entre
o resultado de uma simulação usual com o resultado da extrapolação na vizinhança deT0.
Da figura (2.3) podemos ver que os pontos obtidos das simulações coincidem com a linha
contínua apenas numa região estreita em torno do valor deT0 (símbolo cheio). Os valores
extrapolados para a estimativa da energia média que estão fora desta faixa não são confiáveis,
o que nos permite concluir que não é possível estimar os observáveis de interesse numa faixa
ampla de temperatura, a partir da extrapolação do resultadode uma única simulação.
Todavia, ainda é possível obter uma estimativa para os observáveis numa faixa ampla, através
do método do histograma simples, pois uma vez conhecida a largura da faixa de temperatura
2.2 O MÉTODO DO HISTOGRAMA 20
para uma boa extrapolação, é possível realizar vários histogramas simples e superpor os resul-
tados das extrapolações. No entanto, os resultados obtidospara o valor médio das grandezas
observadas seria um tanto duvidoso, pois este procedimentonão garante o casamento perfeito
entre uma extrapolação e outra. Por outro lado, apresentamos ainda na figura (2.3) a estimativa
para a energia média numa faixa ampla deT (linha tracejada), obtida a partir da combinação
entre vários histogramas simples, realizados para cada valor deT representados pelos símbolos.
Note que todos os pontos coincidem com a linha tracejada. Este resultado foi obtido através do
método do histograma múltiplo, o qual será descrito de formadetalhada na seção a seguir.
2.2.2 O método do histograma múltiplo
O método de histograma múltiplo é uma generalização das idéias envolvidas no método
do histograma simples. Basicamente, ele combina eficientemente o resultado de simulações
realizadas em um dado conjunto de temperaturas em um único histograma, a partir do qual
as extrapolações são efetuadas. Através do método do histograma múltiplo é possível obter
a partir de extrapolações numéricas uma boa estimativa parao valor médio de um observável
numa faixa ampla de temperatura. A forma com a qual esta estimativa é feita é o que iremos
discutir nesta seção.
O procedimento realizado no método do histograma múltiplo não consiste simplesmente em
construir vários histogramas simples e depois superpor os resultados obtidos das extrapolações.
Este procedimento não é recomendável, uma vez que cada simulação fornece uma estimativa
confiável para o valor médio de um dado observável apenas em uma faixa de temperatura em
torno da temperatura na qual os dados foram obtidos. Logo, asestimativas para o valor médio
na vizinhança de uma temperaturaT, resultantes de extrapolações oriundas de simulações re-
alizadas em duas temperaturas distintas, podem diferir significativamente. De fato, considere
que produzimos dois histogramas a partir de simulações realizadas nas temperaturasT1 e T2,
respectivamente, e que as estimativas sejam realizadas numa temperatura na vizinhaça deT1.
2.2 O MÉTODO DO HISTOGRAMA 21
-16000 -14000 -12000 -10000Energia (J)
0
500
1000
H
isto
gram
a
E’
N2(E’)
N1(E’)
T = 2.80
T = 2.84
T = 2.88
Figura 2.4 Histogramas simples resultantes de simulações realizadas nas temperaturasT = 2.80, 2.84e 2.88 para um filme de tamanhoL = 50 com espessuraLZ = 5.
A estimativa obtida da extrapolação a partir do histograma gerado emT1 deve ser mais acurada
que aquela obtida a partir do outro histograma. Portanto, não podemos esperar o casamento
perfeito das curvas, que representam os valores médios dos observáveis, resultantes das duas
extrapolações.
Na figura (2.4) apresentamos os histogramas resultantes de simulações realizadas em três
temperaturas relativamente próximas uma das outras. Podemos observar faixas de energias nas
quais ocorre uma superposição significativa entre os histogramas e outras regiões onde a altura
de dois dos histogramas é praticamente nula.
No exemplo da figura (2.4) vemos que há intersecções entre os intervalos de energias ge-
radas durante as simulações. É intuitivo supor que a faixa detemperaturas em que é seguro
realizar a extrapolação a partir de um histograma possa se superpor com a correspondente faixa
de temperaturas de um outro histograma. Logo podemos obter várias estimativas para o valor
médio de um observável na região de intersecção. Uma combinação desses valores deve resul-
tar numa estimativa mais precisa para o valor médio de interesse. É esta a idéia subjacente ao
2.2 O MÉTODO DO HISTOGRAMA 22
método de histograma múltiplo.
Embora a discussão do parágrafo anterior sugira que pode-seobter boas estimativas para
o valor médio de um observável a partir de alguma combinação de valores extrapolados de
diversas simulações, esta não é a melhor forma de proceder. Como veremos mais abaixo, é
mais interessante combinar os dados das diversas simulações, tentar obter uma estimativa da
densidade de estados e realizar a extrapolação a partir desta.
Por exemplo, na figura (2.4) vemos que há duas alturas distintasN1(E′) e N2(E′) para o
mesmo valor de energiaE′, ou seja, o valor deE = E′ ocorreuN1 vezes no primeiro histograma
enquanto que no segundo foi registradoN2(< N1) vezes. A questão que se coloca é: Qual destes
dois histogramas deve ser considerado na hora de extrapolaros resultados? Melhor ainda,
como podemos aproveitar os dados oriundos das duas simulações? Devemos considerar os
dois resultados, ponderando-os adequadamente. Como fazê-lo constitui a tarefa da técnica de
histograma múltiplo e a figura (2.5) ilustra o resultado da combinação de várias simulações para
construir um único histograma a partir do qual valores médios podem ser obtidos. Observamos
no histograma múltiplo que todos os microestados são bem amostrados em uma ampla faixa de
energias.
Podemos ver da figura (2.5) que o histograma múltiplo apresenta uma altura mais elevada
que a dos histogramas simples, e uma largura que pode aumentar com o número de histogra-
mas simples. Se não há intersecção das faixas de energias cobertas nas simulações individuais,
ou em outras palavras, os valores de temperatura para os quais se realizou as simulações não
são próximos o suficiente, o histograma múltiplo apresentará regiões onde alguns microestados
são mais pobremente amostrados. Este último caso é ilustrado na figura (2.6) onde mostra-
mos histogramas gerados pela combinação de duas e três simulações. Podemos notar que, ao
acrescentarmos dados da simulação realizada na temperatura intermediária, o histograma forma
praticamente um plateau na faixa completa de energias visitadas.
O método do histograma múltiplo baseia-se numa observação relativamente simples e intui-
2.2 O MÉTODO DO HISTOGRAMA 23
-16000 -14000 -12000 -10000E (J)
0
500
1000
1500
2000
His
togr
ama
N(E
)
Histograma Multiplo
T =
2.80
T =
2.82
T =
2.84
T = 2.86
T =
2.88
T = 2.90
Figura 2.5 Histograma múltiplo resultante da combinação dos resultados de simulações paraum filmecomL = 50 eLZ = 5, em várias temperaturas correspondentes aos histogramas simples especificados nafigura.
tiva. Durante o curso de uma simulação Monte Carlo a probabilidade de se observar o sistema
em uma configuração com energiaE é dada pela expressão (2.11), a qual reproduzimos aqui
por conveniência,
Pβ (E) =g(E)e−βE
Zβ. (2.21)
Assim, uma estimativa para a densidade de estados pode ser obtida a partir de
g(E)Zβ
=N(E)
nωβ (E), (2.22)
com ω(E) = exp(−βE). Na última expressão tomamosPβ (E) ≈ N(E)/n, ondeN(E) é o
número de configurações com energiaE dentre asn configurações independentes geradas na
simulação. Observe que tantog(E) quantoZβ são funções desconhecidas. Em princípio, a
densidade de estados poderia ser obtida através de umaúnicasimulação, extremamente longa,
2.2 O MÉTODO DO HISTOGRAMA 24
-10000 -9000 -8000 -7000 -6000E
0
500
1000
1500
N(E
)
T = 2.75, 2.78 e 2.81
T = 2.75 e 2.81
Figura 2.6 Histogramas múltiplos gerados a partir da combinação dos resultados de simulações emduas e três temperaturas, para um filme com espessuraLZ = 4 e estensãoL = 50.
realizada em uma dada temperatura. Como vimos na discussão que se seguiu à apresentação
do método do histograma simples, este procedimento não funciona na prática. PoisN(E) decai
muito rapidamente para energias muito diferentes de〈E(β )〉.
Não obstante, podemos realizar um certo número de simulações Monte Carlo em um dado
conjunto de valores distintos de temperaturas. Neste caso,o mesmo valor de energia pode ser
observado um número arbitrário de vezes em qualquer uma das simulações, c.f. figura (2.5).
Suponha agora que, na i-ésima simulação de comprimentoni, realizada emβi, a energiaE
tenha sido observadaNi(E) vezes. Segue portanto que
gi(E)Zi
=Ni(E)
niωi(E), (2.23)
é uma estimativa para a densidade de estadosg(E). Observe que usamos a notaçãoZ(βi) = Zi
e ωβi = ωi.
É importante ressaltar que temos tantas estimativas para a densidade de estados quanto for o
número de temperaturas. Além disso, comog(E) não depende da temperatura e só depende do
2.2 O MÉTODO DO HISTOGRAMA 25
sistema em estudo, cada uma das estimativasgi(E) é de fato uma estimativa da mesma função.
O problema que se coloca agora é o seguinte. Dado que temos várias estimativas para a
mesma grandeza, como podemos construir a melhor estimativapara a densidade de estados
verdadeira? Retornemos por um instante ao exemplo da figura (2.5), onde mostramos os his-
togramasNi(E) medidos em seis simulações de um filme, realizadas em um conjunto de seis
temperaturas distintas. Vemos que cada histograma individual cobre uma certa faixa de ener-
gias, embora existam intersecções entre alguns intervalos. Esperamos, portanto, que a equação
(2.23) forneça uma boa estimativa para a densidade de estados na faixa de energias na qual o
correspondente histogramaNi(E) assuma um valor significativamente alto. E, ao contrário, a
estimativa (2.23) deve ser muito ruim ondeNi(E) assumir valores muito baixos. Gostaríamos,
portanto, de combinar as diversas estimativas (2.23) para obter a melhor estimativa parag(E)
em todo intervalo de energias coberto por nossas simulações, porém levando em considera-
ção a importância relativa de cada histograma. A maneira formal de proceder é realizar uma
média ponderada sobre as realizaçõesgi(E). Naturalmente, devemos atribuir maior peso ao
gi(E) cujo correspondente histogramaNi(E) estiver melhor amostrado neste particular valor
deE. De fato, a única fonte de erro estatístico associado a uma particular estimativagi(E) está
diretamente ligada ao erro cometido ao medirNi(E), uma vez queZi é apenas uma constante
de normalização quandoβi é fixo. Vamos considerar, agora,τ configurações independentes
geradas em uma simulação Monte Carlo comβ = βi . Sep é a probabilidade de se observar o
sistema simulado em uma configuração cuja energia sejaE, então a média da altura do histo-
grama, tomada sobre um número grande simulações, todas realizadas com o mesmoβi , é dada
por
Ni(E) = τ p (2.24)
e a variância é
δ 2Ni(E) = τ p(1− p). (2.25)
2.2 O MÉTODO DO HISTOGRAMA 26
Como o número de possíveis energias é muito grande, de fato cresce exponencialmente com o
número de partículas do sistema simulado,p≪ 1 e as equações (2.24) e (2.25) implicam que
δ 2Ni(E) = Ni(E). (2.26)
Ou seja, é razoável assumir que se considerarmos um número grande simulações, todas reali-
zadas com o mesmoβi , a distribuição deNi(E), paraE fixo, será uma distribuição de Poisson
com excelente aproximação. E segue da equação (2.23) que o histograma idealNi(E) está re-
lacionado com a densidade de estados exata. De fato, tomando-se a média da equação (2.23)
sobre o mesmo conjunto de simulações, temos
g(E) = gi(E) =Ni(E)Ziniωi(E)
, (2.27)
Para estimar a densidade de estadosg(E) vamos escrevê-la como uma combinação linear
degi , ou seja
g(E) = ∑i
αi(E)gi(E), (2.28)
com o seguinte vínculo para os pesos
∑i
αi = 1. (2.29)
Por exemplo, na figura (2.4) a densidade de estadosg estimada nas simulações emT1, T2 eT3 é
uma combinação linear deg1, g2 eg3 estimadas a partir deN1, N2 eN3, sendo considerados os
seus pesos estatísticos,α1, α2 e α3, respectivamente.
Estamos interessados na estimativa deg(E), cuja variância seja mínima com relação ao
conjunto de coeficientesαi(E). A variância deg(E) é dada por
2.2 O MÉTODO DO HISTOGRAMA 27
δ 2g = ∑i
α2i (g2i −gi2), (2.30)
queremos minimizar a grandeza (2.30) sujeito ao vínculo∑i αi = 1. Para tanto, vamos reescre-
ver (2.30) como
δ 2g = ∑i
αi2σi2, (2.31)
com σi2 = g2i − gi2. Estamos interessados nos valores deαi que minimizam o valor deδ 2g.
Para satisfazer o vínculo (2.29), vamos introduzir o multiplicador de Lagrangeλ , e minimizar
a funçãoG dada por
G(α1,α2, . . .) = ∑i
αi2σi2−λ (∑i
αi −1). (2.32)
A condição de mínima variância é satisfeita quando∂G/∂αi é nula. O que nos dá
αi =λ/2σi2
. (2.33)
Somando (2.33) emi e usando o vínculo (2.29) obtemos, portanto
αi =1/σi2
∑ j 1/σ j2. (2.34)
Substituindo (2.34) em (2.28), podemos finalmente escrevera estimativa para a densidade de
estados, como
g(E) =∑i gi(E)/σi2
∑i 1/σi2. (2.35)
Tomando a média do quadrado da expressão (2.23) e subtraindodo resultado o valor médio
da mesma equação, obtemos
σi2 =[
N2i (E)−Ni(E)2]
[
Ziniωi(E)
]2
=Ni(E)
n2i
[
Ziωi(E)
]2
. (2.36)
2.2 O MÉTODO DO HISTOGRAMA 28
Na última passagem usamos o fato estabelecido na equação (2.26). Observe que a última
igualdade em (2.36) envolve o histogramaideal, ou seja, o número médio de vezes que um dado
valor de energia é realizado em um conjunto enorme de simulações executadas nas mesmas
condições.
Neste caso, a equação (2.27) permite reescrever (2.36) como
σi2 =g2(E)
Ni(E), (2.37)
na qualg(E) é a densidade de estadosexata. A substituição deste último resultado na eq.
(2.35), fornece a nossa melhor estimativa para a densidade de estados como
g(E) =∑i gi(E)Ni(E)/g2(E)
∑i Ni(E)/g2(E)=
∑i gi(E)Ni(E)∑i Ni(E)
. (2.38)
Finalmente, aplicando outra vez a expressão (2.27) podemoseliminarNi(E), pois não conhe-
cemos seu valor, e escrever a densidade de estados como
g(E) =∑i Ni(E)
∑ j n je−β jE/Z(β j). (2.39)
A função de partição do sistema,Z(β ) = ∑E g(E)e−βE, é estimada através da expressão
Z(β ) = ∑E
∑i Ni(E)∑ j n je(β−β j )E/Z(β j)
. (2.40)
E o valor médio de qualquer observável físico é calculado através da seguinte equação
〈A〉β =1
Z(β ) ∑EĀEg(E)e
−βE
=1
Z(β ) ∑EĀE
∑i Ni(E)∑ j n je(β−β j)E/Z(β j)
(2.41)
2.2 O MÉTODO DO HISTOGRAMA 29
A equação (2.41), relação fundamental da técnica do histograma múltiplo, depende dos va-
lores da função de partição,Z(β j), nas temperaturas nas quais as simulações foram realizadas.
Para estimar tais valores, usamos a equação (2.40) paraβ = βk ondek = 1, 2, . . . rotula os
valores do parâmetroβ empregado nas simulações. Ou seja,
Z(βk) = ∑E
∑i Ni(E)∑ j n je(βk−β j )E/Z(β j)
. (2.42)
A implementação do método do histograma múltiplo envolve resolver iterativamente o sistema
de equações (2.42). O procedimento consiste em propor algumconjunto de valores iniciais
para osZk, inseri-los nas equações (2.42) e obter novos valores. Apósvárias repetições, o
novo conjunto deZk difere pouco do anterior. O conjunto de valores para o qual ocorre a
convergência constitui o conjunto solução do sistema.
CAPÍTULO 3
Resultados e Discussões
3.1 Procedimentos computacionais
Modelamos os filmes por estruturas cúbicas simples de dimensões linearesLX, LY e LZ,
comLX = LY = L, maior que a espessuraLZ. Associamos a cada vérticei da rede uma variável
de spin,σi, que pode assumir os valores−1,0 e +1. Os spins interagem de acordo com o
hamiltoniano
H = −J ∑〈i, j〉
σiσ j , (3.1)
onde a soma se estende sobre os pares de spins vizinhos mais próximos eJ > 0 é a interação
ferromagnética de troca. Quandoσi assume valor+1 ou−1 significa que o spin da partícula
aponta na direção do eixo de fácil magnetização do filme, com sinal positivo (negativo) indi-
cando spin para cima (baixo), e o valor deσi igual a zero indica que o spin aponta em qualquer
direção perpendicular ao eixo de fácil magnetização. A geometria dos filmes é ilustrada na
figura (3.1).
123 Ly
Lz
Lx Lx = Ly = L
Figura 3.1 Filmes magnéticos formados por planos atômicos, de ladoL e espessuraLZ. O valor deLZindica o número de planos atômicos da camada magnética.
30
3.1 PROCEDIMENTOS COMPUTACIONAIS 31
Adotamos condições de contorno periódicas nas direções longitudinais (para reduzir os efeitos
de borda nos planos dos filmes) e condições de contorno abertas na direção transversal ao
filme. Isto é, uma partícula que se encontra numa posição qualquer de um plano nas camadas
atômicas mais internas do filme interage com seis partículasvizinhas, enquanto que partículas
localizadas em pontos quaisquer das superfícies que delimitam o filme possuem cinco sítios
vizinhos mais próximos. Em particular, paraLZ = 1, caso em que temos apenas uma única
camada atômica, ou uma rede quadrada, cada uma das partículas possuem apenas quatro sítios
vizinhos. Quando o filme é composto por duas camadas atômicas, como no exemplo ilustrado
na figura (3.2), qualquer uma das partículas possue cinco sítios vizinhos mais próximos.
Figura 3.2 Estrutura cúbica simples em geometria de filme composta por duas camadas atômicas. Ospontos representam as partículas de spinS= 1 que residem nos sítios da rede.
Por conveniência, dividimos as simulações em três etapas. Na primeira etapa realizamos
o procedimento que gera configurações de equilíbrio do sistema de acordo com o algoritmo
de Metropolis e registramos em um arquivo para análise posterior a energia,E(s), e a magne-
tização,M(s), de cada configuração de equilíbrio. Para cada sistema de volume L× L× LZ,
realizamos simulações Monte Carlo em um conjunto de temperaturas em um certo intervalo
[Tmin : Tmax]. Inicialmente, a rede foi posta em uma configuração totalmente desordenada (pa-
râmetro de ordem nulo), na qual o spin de cada partícula assume um dos valores−1,0 e+1
com igual probabilidade. A rede foi, então, percorrida de forma aleatória, e os spins atualizados
3.1 PROCEDIMENTOS COMPUTACIONAIS 32
segundo a regra de Metropolis. As 104 primeiras configurações geradas foram descartadas para
permitir o sistema termalizar, sendo as medidas realizadasnos 106 MCS subsequentes. Nesta
fase, a cada varredura de rede calculamos a energia interna do sistema,E(s), de acordo com a
expressão (3.1), e a magnetização total,M(s), de acordo com
M(s) =N
∑i
σi, (3.2)
ondeN = L2LZ denota o número de spins do sistema. Com este procedimento temos o registro
dos valores deE(s) e M(s) para cada uma das configurações geradas em um dado valor de
temperatura. Assim, temos séries temporais compostas de resultados em 106 configurações, o
que nos dá a liberdade de escolher qualquer subconjunto delas para proceder nossa análise. Em
particular, podemos calcular as correlações temporais entre as configurações e estimar o tempo
de autocorrelação. O tempo de CPU empregado neste procedimento é proporcional ao número
de spins do sistema, ao número de temperaturas, ao número de MCS e depende fortemente do
tipo de processador utilizado. Algumas de nossas simulações duraram vários dias, enquanto
que outras não ultrapassaram uma hora.
Na segunda etapa, processamos os dados obtidos para realizar a análise de histograma múl-
tiplo. Nesta fase são gerados os histogramas e o cálculo das médias microcanônicas das gran-
dezas de interesse. Explicitamente, calculamos as médias microcanônicas, equação (2.10),
(MkEl )E com k, l = 0,1,2, ..., necessárias para o cálculo dos valores médios em função da
temperatura de acordo com a expressão fundamental do métodode histograma múltiplo (2.41).
Dessa forma, obtemos as correlações do tipo〈MkEl 〉 para qualquer temperatura no intervalo
[Tmin : Tmax]. Algumas destas correlações serão úteis para a análise do comportamento crítico
dos filmes.
Com os resultados obtidos na segunda etapa, calculamos finalmente os valores médios dos
3.1 PROCEDIMENTOS COMPUTACIONAIS 33
observáveis de interesse a partir das correlacões entre energia e magnetização, como
Ck,l = 〈MkEl 〉. (3.3)
3.1.1 As grandezas termodinâmicas
A energia média por partícula é obtida fazendok = 0 e l = 1 em (3.3), isto é,
〈E〉N
=C0,1N
=−JN
〈
∑i, j
σiσ j
〉
, (3.4)
onde a última igualdade foi obtida usando a expressão (3.1).De acordo com (3.2) e (3.3) para
k = 1 e l = 0, a magnetização médiam= 〈|M(s)|〉/N é dada por
m=C1,0N
=1N
〈∣
∣
∣
∣
∣
N
∑i
σi
∣
∣
∣
∣
∣
〉
, (3.5)
onde consideramos emC1,0, o módulo do somatório dos spins devido às configurações que
resultam em valores de magnetização com sinais contrários ocorrerem com igual probabilidade.
As médias em (3.4) e (3.5) são tomadas de acordo com o método dohistograma múltiplo através
da expressão (2.41).
A partir destas grandezas termodinâmicas obtivemos como funções contínuas deT, as fun-
ções respostas e algumas derivadas em relação à temperaturanecessárias para analisar os efei-
tos de tamanho finito e estimar os expoentes críticos. Por exemplo, estamos interessados na
derivada demem relação aT, dada por
∂m∂T
=1
NkBT2(C1,1−C1,0C0,1), (3.6)
uma vez que esta grandeza apresenta um máximo na região crítica, e o valor deT para o qual
este máximo ocorre está relacionado com a temperatura pseudo-crítica associada à magnetiza-
3.1 PROCEDIMENTOS COMPUTACIONAIS 34
ção.
Outras grandezas de interesse, que envolvem momentos da magnetização, são o cumulante
de Binder [37] e a sua derivada em relação à temperatura, dadaspelas expressões
u = 1− 〈M4〉
3〈M2〉2 = 1−C4,03C22,0
(3.7)
e∂u∂T
=
(
1NkBT2
)
(1−u){
1C4,0
(C4,1−C4,0C0,1)−2
C2,0(C2,1−C2,0C0,1)
}
. (3.8)
Estas grandezas, como veremos nas secções seguintes, serãoúteis para estudar o comporta-
mento crítico dos filmes.
3.1.2 As funções respostas
Através das funções respostas medimos as flutuações das grandezas macroscópicas. Por
exemplo, o calor específicoc = ∂ 〈E/N〉/∂T está relacionado com as flutuações em energia e
pode ser obtido a partir da relação
c =1
NkBT2(C0,2−C20,1). (3.9)
Já a susceptibilidade magnéticaχ, que fornece as flutuações associadas à magnetização é cal-
culada a partir da relação
χ =1N
(〈M2〉−〈M〉2) = 1N
(C2,0−C21,0). (3.10)
Calculamos também a derivada da susceptibilidade em função da temperatura, conforme a
expressão
3.2 COMPORTAMENTO TERMODINÂMICO 35
∂ χ∂T
=1
N2kBT2[C2,1−C2,0C0,1−2(C0,1C1,1−C20,1C1,0)]. (3.11)
A partir da análise destas grandezas, obtivemos o comportamanto crítico e termodinâmico
dos filmes em função da espessura. Em particular, estudamos como o valor da temperatura
crítica para um filme se aproxima do valor da temperatura crítica do sistema volumétrico.
3.2 Comportamento termodinâmico
As curvas apresentadas nesta subsecção para o comportamento termodinâmico das gran-
dezas físicas, como funções contínuas da temperatura, foram obtidas através do método do
histograma múltiplo, de acordo com a expressão (2.41), conforme mencionamos na secção
anterior.
Como vimos, o cálculo da função de partição presente na equação (2.41) consiste na re-
solução do sistema de equações em (2.42). Adotamos uma formaiterativa para resolver este
sistema de equações. Ou seja, propomosZ(0)i = 1 parai = 1,2, . . ., como solução. Uma solução
melhorada,Z(1)i é obtida inserindo-se a solução corrente do lado direito do sistema de equa-
ções (2.42). O processo é iterado até que a convergêcia seja atingida ou o número máximo de
iterações seja ultrapassado. Especificamente, fixamos o número máximo de iterações em 200 e
uma tolerância|Z(n+1)i /Z(n)i −1| < ε = 10−7.
No geral a convergência é rápida (ocorre bem antes das 200 iterações), e o tempo de pro-
cessamento empregado nesta etapa é, no pior dos casos, da ordem de segundos. No caso de não
haver a convergência em duzentas iterações, significa que astemperaturas nas quais realizamos
as simulações não são próximas o suficiente, e a própria análise do histograma indica este fato.
Nestes casos, realizamos simulações em temperaturas intermediárias e procedemos uma nova
análise dos dados.
3.2 COMPORTAMENTO TERMODINÂMICO 36
2,8 2,85 2,9 2,95Temperatura (J/k
B)
-1,2
-1
-0,8
Ene
rgia
med
ia (
J)
L = 50L = 60L = 70L = 80L = 90L = 100
LZ = 5
(a)
3,14 3,16 3,18 3,2Temperatura (J/K
B)
-0,85
-0,8
-0,75
-0,7
Ene
rgia
med
ia (
J)
L = 50L = 60L = 70L = 80L = 90L = 100
LZ = 18
(b)
Figura 3.3 Energia média em função da temperatura para vários tamanhos de filmes com espessurasLZ = 5 eLZ = 18.
Apresentamos apenas os resultados obtidos para filmes com espessurasLZ = 5 eLZ = 18.
Os resultados das simulações para outros valores deLZ são qualitativamente muito semelhantes
e deixamos para as tabelas e diagrama de fases a inclusão destes.
Para os filmes comLZ = 5, realizamos simulações MC para valores deT entre 2.8 e 2.96,
com intervalo de variação∆T = 0.01. ParaL = 70,80,90,100 simulamos o filme também
em temperaturas intermediáriasT = 2.875,2.885,2.891,2.892 e 2.895, com os dois últimos
valores de temperatura simulados apenas para os tamanhos defilmes comL = 90 e 100. Os
valores intermediários deT em que foram realizadas as simulações para os tamanhos maiores
são pontos pertencentes à região crítica. De fato, boas simulações MC exigem a realização de
simulações em um número maior de pontos na região crítica, devido às grandes flutuações nesta
região. Para os filmes comLZ = 18, realizamos simulações no intervalo deT = 3.13 aT = 3.20,
com intervalo de variação∆T = 0.01, e em temperaturas intermediáriasT = 3.135,3.145 e
3.155 para todos os tamanhos de filmeL.
Na figura (3.3) apresentamos a energia média como uma função da temperatura, paraLZ = 5
e 18, com valores deL variando entre 50 e 100. Observamos uma dependência muito pequena
com a extensãoL dos filmes, o que está de acordo com o que é observado para redesregulares
quadradas e cúbicas. Esta dependência com o tamanhoL, ocorre apenas em uma faixa estreita
3.2 COMPORTAMENTO TERMODINÂMICO 37
de temperatura e, se olharmos numa escala mais apropriada para esta região, podemos perceber
que ela se torna mais acentuada à medida que a extensão do filmeaumenta. Este fato é melhor
observado para o casoLZ = 18. A forma com a qual ocorre a dependência comL significa que
as flutuações em energia aumentam à medida queL cresce, e isto tem repercussão direta no
comportamento do calor específico apresentado a seguir. Obtivemos, como esperado, resulta-
dos semelhantes do comportamento da energia média em funçãodeT para filmes com outras
espessuras, sendo a única diferença entre os resultados a localização da faixa de temperaturas
na qual o sistema apresenta dependência com o tamanho.
2,8 2,85 2,9 2,95Temperatura (J/k
B)
10
15
20
25
30
35
40
45
Cal
or e
spec
ifico
L = 50L = 60L = 70L = 80L = 90L = 100
LZ = 5
(a)
3,14 3,16 3,18Temperatura (J/K
B)
10
20
30
40
Cal
or e
spec
ifico
L = 60L = 70L = 80L = 90L = 100
LZ = 18
(b)
Figura 3.4 Calor específico em função da temperatura para vários tamanhosL de filmes com espessurasLZ = 5 eLZ = 18.
Na figura (3.4) apresentamos o comportamento do calor específico como uma função deT,
com vários valores deL, para os mesmos filmes. As faixas de temperaturas nas quais a energia
apresenta uma dependência com o tamanho dos filmes, figuras (3.3-a) e (3.3-b), correspondem
às regiões de máximo do calor específico nos dois casos. O valor de T para o qual o calor
específico apresenta o seu valor máximo é o que chamamos de temperatura pseudo-crítica
associada ao calor específico,T∗c , a qual se aproxima deTC no limite termodinâmico. Para
um sistema magnético Ising por exemplo, é observado um deslocamento do máximo do calor
específico em função do tamanho do sistema [38]. Isto é, conforme aumentamos o tamanho do
3.2 COMPORTAMENTO TERMODINÂMICO 38
sistema, o valor máximo do calor específico desloca-se para aesquerda.
Note da figura (3.4) que o calor específico apresenta para cadavalor deL um máximo,
cuja amplitude aumenta com o tamanho dos filmes. Todavia, esta dependência comL é muito
pequena em ambos os casos, de acordo com a dependência observada para a energia média.
Este comportamento é típico de um sistema bidimensional. É bem conhecido para um sistema
Ising-2D [39] que o expoente crítico associado ao calor específico é nulo, ou seja, a divergência
em TC é logarítmica,c ∼ ln(L). Temos então uma forte indicação da natureza bidimensional
dos filmes.
Na figura (3.5) mostramos o comportamento da magnetização emfunção da temperatura
para filmes com espessurasLZ = 5 eLZ = 18, para alguns valores deL. Vemos que a magne-
tização do sistema apresenta uma dependência com o tamanho do filme, apenas acima de uma
dada temperatura. A região de temperatura na qual a magnetização apresenta uma dependência
considerável com a temperatura corresponde à região crítica. Efeitos de tamanho finito são ob-
servados, de forma que o decaimento das curvas é bem mais lento para os tamanhos menores.
Estes comportamentos são mais facilmente observados nas curvas da susceptibilidade. O valor
deT correspondente ao ponto de inflexão nas curvas das figuras (3.5-a) e (3.5-b) depende do
tamanho dos filmes, de forma que quanto maior é o valor deL, menor é o valor deT no qual
ocorre a inflexão. Logo, quanto maior é o tamanho do filme, paraambos os valores deLZ,
mais próximo estamos do valor da temperatura crítica no limite termodinâmico. O ponto de
inflexão nas curvas da magnetização, correspondente a cada valor deL, para ambos os filmes
está associado à temperatura pseudo-crítica,T∗m.
Apresentamos na figura (3.6) a susceptibilidade magnética em função da temperatura para
vários tamanhos de um filme com espessuraLZ = 5 e 18. A susceptibilidade magnética, que
como vimos fornece as flutuações da magnetização, apresentauma grande dependência comT
na região crítica, e quase nenhuma dependência fora desta região. Ainda nesta região os filmes
apresentam uma forte dependência com o tamanhoL, como observado para a magnetização.
3.2 COMPORTAMENTO TERMODINÂMICO 39
2,8 2,85 2,9 2,95Temperatura (J/k
B)
0
0,1
0,2
0,3
0,4
0,5
Mag
netiz
acao
L = 50L = 60L = 70L = 80L = 90L = 100L
Z = 5
(a)
3,14 3,16 3,18 3,2Temperatura (J/k
B)
0
0,05
0,1
0,15
0,2
0,25
Mag
netiz
acao
L = 50L = 60L = 70L = 80L = 90L = 100
LZ = 18
(b)
Figura 3.5 Magnetização como uma função contínua deT para vários tamanhos de filmes com espes-surasLZ = 5 eLZ = 18.
2,8 2,85 2,9 2,95Temperatura (J/k
B)
0
100
200
300
400
500
600
700
800
Sus
cept
ibili
dade
L = 50L = 60L = 70L = 80L = 90L = 100
LZ = 5
(a)
3,14 3,16 3,18 3,2Temperatura (J/k
B)
0
200
400
600
800
Sus
cept
ibili
dade
L = 50L = 60L = 70L = 80L = 90L = 100
LZ = 18
(b)
Figura 3.6 Susceptibilidade magnéticaχ em função da temperatura para vários valores deL para filmescom espessurasLZ = 5 eLZ = 18.
Se compararmos as figuras (3.6-a) e (3.6-b), podemos ver que,para um dadoL, o filme com
espessura maior apresenta um valor máximo da susceptibilidade maior. A temperatura pseudo-
crítica,T∗χ , é o valor deT no qualχ é máximo para cadaL.
3.3 COMPORTAMENTO CRÍTICO 40
3.3 Comportamento crítico
Na subsecção anterior discutimos o comportamento termodinâmico das grandezas macros-
cópicas e das funções respostas para vários tamanhosL de filmes com espessurasLZ = 5 e 18.
Observamos os efeitos de tamanho finito para estas grandezase definimos de forma qualitativa
a temperatura pseudo-crítica associada