ads:
ESTUDO TEÓRICO DA DEPOSIÇÃO DE UMA CAMADA DE
NITROGÊNIO SOBRE UMA SUPERFÍCIE (100) DE ARSENETO
DE GÁLIO
ALDA DE PAIVA CASTRO
ads:
Livros Grátis
http://www.livrosgratis.com.br
Milhares de livros grátis para download.
i
Estudo Teórico da Deposição de uma Camada de Nitrogênio
sobre uma Superfície (100) de Arseneto de Gálio
Alda de Paiva Castro
DISSERTAÇÃO SUBMETIDA AO CORPO DOCENTE DA COORDENAÇÃO DO
PROGRAMA DE PÓS-GRADUAÇÃO MULTIDISCIPLINAR EM FÍSICA,
QUÍMICA E NEUROCIÊNCIAS DA UNIVERSIDADE FEDERAL DE SÃO JOÃO
DEL REI COMO PARTE DOS REQUISITOS NECESSÁRIOS PARA A OBTENÇÃO
DO GRAU DE MESTRE EM CIÊNCIAS, ÁREA DE FÍSICA, EM TEORIA
QUÂNTICA DE MATERIAIS.
Aprovada por:
SÃO JOÃO DEL REI, MG – BRASIL
NOVEMBRO DE 2005
ads:
ii
CASTRO, ALDA DE PAIVA
Estudo Teórico da Deposição de uma
Camada de Nitrogênio sobre uma superfície
(100) de Arseneto de Gálio [Minas Gerais] 2005
VII, 84 p. 29,7 cm (UFSJ, M.Sc., Física,
2005)
Dissertação – Universidade Federal de São
João del Rei, Multidisciplinar em Física,
Química e Neurociências
1. Nitridação de uma Superfície (100) de GaAs
I. UFSJ II. Título (série)
iii
Agradecimentos
Agradeço a todas as pessoas que, direta ou indiretamente, contribuíram para a
conclusão desta dissertação de mestrado. De modo especial, gostaria de agradecer:
• A Deus, que esteve sempre presente em minha vida, pela luz, pela força e
perseverança ao longo deste período;
• Ao meu orientador, pela orientação e oportunidade de realizar este trabalho;
• A minha família, pelo apoio, aconchego e incentivo;
• Ao Ronaldo, meu grande companheiro, por todo apoio, paciência e amor;
• A todos os professores do DCNAT, pelo conhecimento que me forneceram
desde o período da graduação;
• Aos meus colegas de mestrado, pela amizade, conversas e trocas de
conhecimentos;
• À CAPES, pelo apoio financeiro.
iv
Resumo da Dissertação apresentada à UFSJ como parte dos requisitos necessários para
a obtenção do grau de Mestre em Ciências (M.Sc.).
ESTUDO TEÓRICO DA DEPOSIÇÃO DE UMA CAMADA DE
NITROGÊNIO SOBRE UMA SUPERFÍCIE (100) DE ARSENETO
DE GÁLIO
Alda de Paiva Castro
Novembro/2005
Orientador: Horácio Wagner Leite Alves
Programa: Multidisciplinar em Física, Química e Neurociências
Neste trabalho usamos a Teoria do Funcional da Densidade na aproximação da
densidade local, expansão da função de onda em ondas planas e o método do
pseudopotencial para determinar os aspectos estruturais e energéticos da superfície
(100) de GaAs antes e depois da adsorção de N sobre esta superfície. As energias totais
obtidas no processo de adsorção e no de difusão do N sobre esta superfície foram usadas
como parâmetros de uma simulação Monte Carlo Cinético deste processo no MBE,
baseado no modelo de gás de rede.
Os nossos resultados mostram que a deposição do N sobre a superfície (100)
(2x1) de GaAs ocorre de duas maneiras. Na primeira, o N se liga ao dímero de As e,
quando chega um segundo átomo de N, realiza uma troca iônica com os átomos de As
na superfície que são dessorvidos. Na outra maneira, o N é adsorvido ao redor do
dímero de As. Neste caso, o N se difunde muito lentamente na superfície, mesmo
quando a temperatura de crescimento é alta. Há também uma troca iônica entre dois
átomos de N com o dímero de As, facilitando sua dessorção. Finalmente, a superfície
nitridada resultante é bem rugosa, devido à formação de várias ilhas antes da ocorrência
da nucleação e formação das camadas depositadas de N, seguindo o modelo de Volmer-
Weber de crescimento.
v
Abstract of Dissertation presented to UFSJ as a partial fulfillment of the requirements
for the degree of Master of Science (M.Sc.)
THEORICAL STUDY OF THE DEPOSITION A MONO LAYER OF
NITROGEN OVER A GALLIUM ARSENIDE (100) SURFACE
Alda de Paiva Castro
November/2005
Advisor: Horácio Wagner Leite Alves
Program: Multidisciplinary in Physics, Chemistry and Neuroscience
We have used in this work the Density Functional Theory in the Local Density
Approximation, plane-wave expansion of the wavefunction and the pseudopotential
method, in order to obtain the structural and energetic aspects of the GaAs (100) surface
before and after the adsorption of N over this surface. The obtained total energies for the
N adsorption and diffusion over the surface were used as parameters of a kinetic Monte
Carlo the simulation this process, in the MBE chamber based in the lattice gas model.
Our results, show that the deposition of N over GaAs (100) (2x1) surface occur
in the two following ways: the N binds to the As dimers and, when the second N
arrives, it occurs the anion exchange between the N and As atoms, culminating with the
As
2
dessorption. In the other way, the N is adsorbed around As dimer. In this case the N
diffuses very slow over the surface, even high growth temperatures. Also, there is an
anion exchange between the two N atoms and the As dimer, favouring their dessorption.
Finally, the resulting surface is very rough, due to the formation of the several islands
before nucleation and formation of the deposited layer of the N, as described by the
Volmer-Weber model of crystal growth.
vi
Índice
1 Introdução.............................................................................................................. 01
2 Métodos Teóricos e Metodologia
2.1 Introdução....................................................................................................... 09
2.2 Teoria do Funcional da Densidade.................................................................. 10
2.2.1 Aproximação da Densidade Local...................................................... 12
2.2.2 Método das Ondas Planas................................................................... 13
2.2.3 Método do Pseudopotencial................................................................ 15
2.2.3.1 Pseudopotencial de Troullier-Martins.................................... 17
2.2.4 Método do Gradiente Conjugado........................................................ 18
2.2.5 Amostragem da Zona de Brillouin..................................................... 19
2.2.6 Forças de Helmann-Feynman............................................................. 21
2.3 Método de Monte Carlo................................................................................ 22
2.3.1 O Método de Monte Carlo e o Algoritmo de Metropolis................... 23
2.3.2 O Método de Monte Carlo Cinético................................................... 26
2.3.3 Modelos de Crescimento Cristalino................................................... 28
2.3.4 Simulações Monte Carlo Cinético para o Crescimento por MBE...... 30
2.4 – Metodologia................................................................................................. 35
3 Reconstrução e Relaxação da Superfície (100) de GaAs
3.1 Superfície de GaAs (100)................................................................................ 42
3.1.1 Propriedades de Equilíbrio Cristalino do GaAs.................................. 43
3.1.2 Aspectos Estruturais da Superfície (100) de GaAs, na reconstrução
(2x1)................................................................................................. 44
3.2 Adsorção de um átomo de N sobre a superfície (100) de GaAs, terminada em
As, na reconstrução (2x1)............................................................................... 46
vii
3.3 Adsorção do segundo N sobre a superfície (100) de GaAs, terminada em As,
na reconstrução (2x1).................................................................................. 50
3.4 Energia de Formação da Adsorção de N na superfície (100) de GaAs,
terminada em As............................................................................................ 53
3.5 Difusão do N................................................................................................... 55
3.6 Simulações Monte Carlo Cinético da Nitridação............................................ 57
4 – Conclusões .......................................................................................................... 64
Apêndices
A Pseudopotenciais para o H, As, Ga e N...................................................................... 66
B Propriedades Eletrônicas da Superfície (100) de GaAs, terminada em As, na
reconstrução (2x1).......................................................................................................... 73
C Propriedades Eletrônicas da Superfície (100) de GaAs, terminada em As, na
reconstrução (2x1) com um Átomo de N Adsorvido..................................................... 75
Referências Bibliográficas.................................................................................... 77
1
Capítulo 1
Introdução
Os nitretos do grupo III (BN, AlN, GaN e InN) e suas ligas têm sido muito
usados na construção de dispositivos eletrônicos que operam sob altas temperaturas
(transistores e diodos de alta potência), e em dispositivos optoeletrônicos que trabalham
na região do visível do espectro eletromagnético, tais como lasers (LDs), diodos
emissores de luz (LEDs), além de detetores de raios X e de alta voltagem. Toda esta
gama de aplicações é possível pois o “gap” de energia (diferença de energia entre o
estado mais alto ocupado e o primeiro estado desocupado) destes materiais é grande e a
transição entre estas bandas é direta. Os valores para estes “gaps” são 0.75 eV
(wurtzita) ou 1.3 eV (cúbica) [1] para o Nitreto de Índio (InN), 3.4 eV (wurtzita) ou 3.2
eV (cúbica) para o Nitreto de Gálio (GaN) e 6.2 eV (wurtzita) ou 5.3 eV (cúbica) para o
Nitreto de Alumínio (AlN), o que inclui toda a região do visível estendendo-se até a
região do ultravioleta (UV). Desta forma, podemos combinar estes materiais para
formar ligas ternárias e/ou quaternárias que, em princípio, torna possível a existência de
qualquer valor para o “gap” de energia dentro deste intervalo, conforme está descrito na
figura 1.1 [2]. Além disto, suas ligações químicas são bem fortes, apresentam valores
para constante dielétrica baixos, condutividade térmica elevada, alto ponto de fusão,
elevada dureza e módulo de compressibilidade elevado [3].
Figura 1.1: Variação do “gap” de energia de GaN, InN e AlN hexagonais e cúbicos e suas ligas em
função do parâmetro de rede.
2
Os primeiros lasers semicondutores foram fabricados no início da década de 60 e
eram constituídos de uma junção do tipo p-n de Arseneto de Gálio (GaAs). As faces
perpendiculares a esta junção eram polidas, formando a cavidade do laser, que foi
denominada de laser de emissão lateral emitindo luz em um comprimento de onda de
0.88 µm. Contudo, o seu uso prático foi inviável, uma vez que sua operação necessitava
de uma corrente elétrica muito alta. Com o aperfeiçoamento da técnica de epitaxia por
fase líquida (LPE, liquid-phase epitaxy) tornou-se possível tanto o crescimento de
heteroestruturas, favorecendo desta forma, a melhoria das características dos materiais
para estes lasers, como também o seu uso, pela primeira vez, em sistemas de
telecomunicações.
Devido ao desenvolvimento dos métodos de crescimento de sólidos, como a
Epitaxia por Feixe Molecular (Molecular Beam Epitaxy - MBE), Deposição por Vapor
Químico de Organo-Metálicos (Metal-Organic Chemical Vapor Deposition - MOCVD),
dentre outros [4], bem como dos métodos de caracterização para controle do
crescimento de materiais, como a Microscopia de Varredura por Tunelamento
(Scanning Tunneling Microscopy – STM) [5], houve um grande progresso na
construção de lasers nos últimos 20 anos.
A técnica MBE possibilita o crescimento ordenado de camadas cristalinas
através do fluxo unidirecional de moléculas ou átomos dos elementos constituintes de
um filme a ser crescido. Possui um controle preciso (em escala atômica) do fluxo dos
feixes, devido às baixas taxas de deposição e das condições de deposição quando a
amostra é colocada numa câmara mantida em ultra-alto-vácuo (UHV) (sob pressões da
ordem de 10
-11
torr). Logo, esta tem sido a mais sofisticada técnica de crescimento.
O princípio de funcionamento do MBE consiste, então, na evaporação de fontes
sólidas altamente purificadas em alto-vácuo, produzindo feixes moleculares
direcionados sobre a superfície aquecida do substrato. O processo básico da evaporação
pela técnica de epitaxia por feixe molecular está esquematicamente ilustrado na figura
1.2. Fornalhas especializadas, denominadas células de Knudsen (ou células-K, ou
células de efusão), induzem um fluxo de moléculas. Este fluxo de moléculas entra na
câmara através de um orifício do obturador apropriado. A temperatura de cada célula-K
é que controla o fluxo de seus constituintes. Um feixe molecular é formado e, pelo
controle deste feixe na câmara e da temperatura do substrato, os materiais podem ser
crescidos com composição e espessura específicas. A técnica MBE permite o controle
adequado da espessura da superfície, possibilitando variações abruptas na rede, uma
3
exigência para o crescimento de heteroestruturas [6,7]. Nesta técnica, os componentes
do material a ser depositado são colocados em células de efusão, e evaporados (ou
sublimados, dependendo do elemento) através de um sistema de filamentos resistivos. A
fim de manter a temperatura constante das células, gerando um fluxo estável de
material, são utilizados controladores eletrônicos de potência e sensores (termopares).
Os obturadores mecânicos, por sua vez, estão localizados na frente das células
permitindo que os feixes sejam interrompidos quando necessário e possibilitando o
crescimento de camadas alternadas de diferentes materiais com interfaces abruptas. As
células devem estar dispostas na câmara de ultra-alto vácuo de modo que todos os feixes
se cruzem na posição em que o substrato é colocado.
(a) (b)
Figura 1.2: Representação esquemática do processo MBE: a) obturadores fechados e b) um obturador é
aberto.
Na superfície do substrato ocorrem vários processos físico-químicos
relacionados ao crescimento, tais como a adsorção de átomos ou moléculas, a difusão
superficial (migração), dissociação de moléculas, incorporação de átomos na rede
cristalina, e a dessorção térmica das espécies não incorporadas. A superfície do
substrato (ou última camada epitaxial crescida) possui várias regiões com altas
probabilidades de interação, as quais são denominadas sítios de superfície, criados por
ligações químicas pendentes ou vacâncias. É nestes sítios que ocorre o processo de
incorporação de átomos, ou moléculas, no substrato conforme é mostrado na figura 1.3.
O substrato tem uma forte influência sobre o filme a ser crescido. Geralmente, o filme
depositado segue, quando possível, a orientação e as características cristalográficas do
substrato. A dessorção térmica das espécies não incorporadas também é de grande
importância no processo de crescimento, uma vez que garante a estequiometria do
material que está sendo crescido.
4
Figura 1.3: Representação esquemática dos processos envolvidos no MBE [8].
A técnica MOCVD está ilustrada na figura 1.4. Como mostrado nesta figura, os
diferentes gases que irão constituir o material, passam através de uma câmara que
contém o substrato aquecido. A composição do material depositado é controlada pela
reação química entre os gases e pela temperatura do substrato. Os gases passam através
da câmara por convecção forçada, favorecendo a precipitação do material a ser crescido
sobre o substrato. Desta forma, a MOCVD consiste em uma técnica de síntese de
materiais semicondutores baseada na reação química de espécies químicas denominadas
precursores. Estes precursores são vapores organo-metálicos, tais como tri-metil-gálio
[Ga(CH
3
)
3
] ou tri-metil-alumínio [Al(CH
3
)
3
], que irão reagir com a amônia (NH
3
) ou
tri-hidreto de Arsênio (AsH
3
), os quais são transportados para dentro de um reator
usando gás portador de hidrogênio, que é mantido a uma alta temperatura, onde os
precursores são decompostos.
Figura 1.4: Representação esquemática da técnica MOCVD.
Após a reação química, os átomos liberados recombinam-se formando o
composto semicondutor. Todo este processo ocorre em substratos colocados sobre um
susceptor radiativamente aquecido dando origem ao crescimento do filme. A equação
5
1.1 mostra uma reação usando o Ga(CH
3
)
3
e o AsH
3
, que é utilizado para o crescimento
do GaAs [6].
(CH
3
)Ga + AsH
3
→ GaAs ↓ + 3CH
3
(1.1)
Durante o crescimento por MOVCD, uma série de reações superficiais podem
ocorrer, como está descrito na figura 1.5.
Figura 1.5: Representação esquemática das reações superficiais que podem ocorrer durante o
crescimento por MOCVD.
Enquanto na MOCVD os precursores são vapores de compostos organo-
metálicos, no MBE, estes são de origem atômica.
A STM é uma técnica que permite a observação, no espaço real, de uma
estrutura superficial com resolução a nível atômico, devido à dependência exponencial
entre a corrente de tunelamento e a distância da ponta para a amostra. Utilizando essa
dependência, a varredura é usualmente executada, mantendo a corrente constante. Além
disto, uma vez que a corrente depende da voltagem da amostra, pode-se também obter
informações sobre os estados eletrônicos da superfície em função da voltagem aplicada.
Na figura 1.6, mostramos uma representação esquemática do aparato STM.
Figura 1.6: Representação esquemática do aparato STM.
6
Usando estas técnicas de crescimento, Dupuis e colaboradores [9] fizeram, em
1978, o primeiro laser à base de poços quânticos de GaAs (como meio ativo), que
necessitavam de menor corrente de operação e possuíam uma eficiência superior aos
lasers de heteroestrutura dupla da época.
Há cerca de dez anos, a fabricação de diodos emissores de luz azul altamente
eficientes [10] e de um laser que emitia luz azul foi conseguida por Nakamura e
colaboradores. Estes usavam nitretos do grupo III (GaN, InN, AlN) como meio ativo
[11]. Neste caso, GaN cúbico foi crescido ou por MOCVD sobre substratos de safira.
Mais tarde, o mesmo tipo de dispositivo foi obtido sobre substratos de Carbeto de
Silício (SiC), Silício (Si) ou GaAs por MBE [12-14]. Além disto, LEDs de
heteroestruturas de c-InGaN/GaN, crescidas sobre GaAs (001), foram obtidos através de
MOCVD [15]. Além do crescimento destes dispositivos, houve um grande progresso na
obtenção das ligas cúbicas de InGaN, InAlN e AlGaN [16-19].
Entretanto, apesar do avanço tecnológico na concepção e fabricação destes
dispositivos optoeletrônicos, pouco se sabe sobre os mecanismos fundamentais do
crescimento epitaxial e como os átomos se distribuem sobre a superfície do substrato.
Apenas em 1997 micrografias de varredura por tunelamento atomicamente resolvidas
(“atomically resolved scanning tunneling micrographs”) foram obtidas para superfícies
do GaN hexagonal [20]. Torna-se necessário, portanto, um conhecimento detalhado das
propriedades e estruturas destas superfícies para a compreensão das características dos
dispositivos e melhoria do seu desempenho.
A fase cristalina mais estável dos nitretos do grupo III é a estrutura wurtzita [2].
Porém, estes materiais podem se cristalizar tanto na estrutura wurtzita quanto na
zincblende. A obtenção destes materiais na estrutura cúbica (zincblende), que levariam
grande vantagem sobre a wurtzita por serem isotrópicos, só é possível através de
fenômenos fora do equilíbrio termodinâmico, que só devem ocorrer no crescimento
epitaxial via MBE [21]. Ainda não há modelos teóricos que consigam descrever com
precisão, o crescimento cristalino por epitaxia, que poderiam auxiliar de sobremaneira,
o melhor uso destas técnicas para o crescimento destes materiais.
GaN cúbico (zincblende) tem sido crescido sobre substratos de GaAs (001)
cúbico [22-24]. FEUILLET e colaboradores [25], em um trabalho sobre o crescimento
via MBE da fase cúbica do GaN sobre a superfície de SiC (001) e de GaAs (001),
observaram um padrão de reconstrução (4x1) para a superfície rica em Ga, e também a
estrutura de reconstrução ideal (1x1) para a superfície rica em N. Padrões de
7
reconstruções (2x2) e c(2x2) foram observadas por outros grupos [24]. FEUILLET e
colaboradores [25] sugerem que estas reconstruções são obtidas devido ao resíduo de As
na câmara do MBE. Para o crescimento de GaN cúbico sobre substratos de GaAs (001)
eles obtiveram superfícies (2x2) ou c(2x2), embora a orientação do crescimento
utilizada fosse a mesma. Estas observações são resumidas na figura 1.7.
Figura 1.7: Resumo das observações experimentais verificadas para o crescimento de GaN cúbico em
dois substratos diferentes com ou sem a presença de As.
Com este cenário, o problema científico e tecnológico de se desenvolver e
aprimorar técnicas para o crescimento dos nitretos na fase cúbica tem atraído a atenção
da comunidade acadêmica nestes últimos anos [25-28]. Em particular, o crescimento do
GaN sobre substratos de GaAs (001) ou de SiC (001), tornou-se um problema
interessante uma vez que as reconstruções na superfície de GaN se alteram
drasticamente dependendo da atmosfera usada em sua deposição, utilizando tanto o
MBE como a MOCVD.
Mesmo com o atual desenvolvimento da Teoria do Funcional da Densidade
(DFT), os resultados teóricos só prevêem a reconstrução (4x1) do GaN (001) como a
mais estável, apesar dos cálculos para as outras reconstruções vistas experimentalmente
mostrarem que suas energias de formação estão bem próximas [29].
Do ponto de vista teórico, o crescimento de materiais por técnicas de epitaxia em
escala atômica, pode ser estudado com alta precisão empregando a DFT por duas vias:
i) do equilíbrio termodinâmico do sistema, em que as reconstruções da superfície do
material crescido são estudadas e comparadas, ii) da cinética de crescimento do sistema,
onde procura-se entender os mecanismos de adsorção, dessorção e difusão de processos
simples, através de cálculos das barreiras de energia de um átomo sobre a superfície em
questão [30, 31, 32, 33]. Embora tenha havido grande progresso no estudo do
crescimento de materiais usando a Teoria do Funcional da Densidade (DFT), estes
cálculos usam até algumas centenas de átomos e as escalas de tempo usadas nas
simulações usando DFT com dinâmica molecular são várias ordens de magnitude
8
menores que as vistas nos processos reais [30]. Devido em grande parte a estas
limitações, têm sido desenvolvidas técnicas de simulação de Monte Carlo Cinético para
descrever o crescimento epitaxial de materiais [30, 31].
Sob este prisma, apresentaremos neste trabalho um estudo da deposição de
Nitrogênio (N) sobre substratos de GaAs (100), através de simulações de Monte Carlo
Cinético [30-35] conjugado com métodos de primeiros princípios [36, 37] para
entendermos como seria possível o crescimento de GaN cúbico. Neste trabalho não
iremos usar a reconstrução β2(2x4) como ponto de partida para a simulação e sim a
(2x1) pois, conforme a figura 1.8, na temperatura de crescimento do GaN sobre GaAs
(~700 °C), a superfície do substrato sofre transição de fase para o padrão de
reconstrução (2x1) [38].
Figura 1.8: Diagrama de fases experimental da superfície (001) de GaAs em câmaras MBE, após
LABELLA e colaboradores [38] em função da presença de As.
Iremos, então, verificar o que acontece quando somente a temperatura é variada
na nitridação da superfície (100) de GaAs, que é considerada o primeiro passo para o
crescimento de Nitretos sobre este substrato.
A apresentação de nosso trabalho está organizada da seguinte maneira. O
capítulo 2 apresenta uma descrição da Teoria do Funcional da Densidade, do método do
pseudopotencial e do método de Monte Carlo. Fazemos também neste capítulo uma
exposição da metodologia por nós usada. No capítulo 3, apresentamos nossos
resultados, sua análise e sua discussão. Finalmente, o capítulo 4 é destinado às
conclusões e considerações finais do trabalho.
9
Capítulo 2
Métodos Teóricos e Metodologia
2.1 Introdução
A determinação das propriedades de um sistema quântico qualquer,
caracterizado por um hamiltoniano H, é baseada na solução da equação de Schrödinger,
que é dada por
Ψ
=
Ψ
E
H
)
, (2.1)
onde
Ψ
é a função de onda do sistema e E é a energia do sistema.
Para átomos, moléculas e sólidos, o hamiltoniano H é escrito da seguinte forma:
∑∑∑∑∑∑∑ ∑
= == >= >= = −
−
−
+
−
+∇−−∇−=
N
A
M
A
Ai
A
M
A
M
AB
ji
BA
N
i
N
j
ji
A
N
i
N
A
A
i
Rr
Z
RR
ZZ
rr
M
H
1 111 1
2
1 1
2
1
2
1
2
1
r
r
rr
rr
, (2.2)
em que M
A
é a razão entre a massa do núcleo A e a massa do elétron, Z
A
é o número
atômico de A. {R
i
}denotam as posições dos núcleos e {r
i
}, as posições dos elétrons. Os
dois primeiros termos representam as energias cinéticas dos elétrons e dos núcleos,
respectivamente. O terceiro representa as repulsões entre os elétrons e o quarto termo
representa as repulsões entre os núcleos. O último termo corresponde à atração entre
elétrons e núcleos. Nesta equação foram usadas unidades atômicas, onde a carga e a
massa do elétron e a constante de Planck foram consideradas unitárias (energia dada em
hartree, e a distância em raios de bohr).
Para resolver a equação 2.1, a primeira aproximação que fazemos é a
aproximação de Born-Oppenheimer, que é baseada na hipótese de que os elétrons se
movem muito mais rápido do que os núcleos, de maneira que os núcleos podem ser
10
considerados aproximadamente fixos. Essa suposição se justifica pela grande diferença
entre as massas dos elétrons e dos núcleos. Na prática, isso significa que a estrutura
eletrônica é calculada para um arranjo fixo e os núcleos, então, se movem de maneira
quase clássica.
Para a obtenção dos autovalores do hamiltoniano (2.2) utilizaremos a Teoria do
Funcional da Densidade (DFT), na Aproximação da Densidade Local (LDA).
2.2 Teoria do Funcional da Densidade
A DFT é utilizada para calcular a estrutura eletrônica de átomos, moléculas e
sólidos. Ela é baseada nos teoremas de Hohenberg e Kohn [39] e no de Kohn e Sham
[40].
O teorema de Hohenberg e Kohn mostra que a energia de um sistema de N
elétrons interagentes sob a ação de um potencial externo
(
)
rv
v
é um funcional único da
densidade eletrônica,
(
)
rn
r
. A energia do estado fundamental E, é em termos
matemáticos, um funcional
(
)
[
]
rnE
r
de
(
)
rn
r
ou seja,
(
)
[
]
(
)
[
]
(
)
(
)
rdrvrnrnFrnE
r
r
r
r
r
∫
+= ,
(2.3)
onde
(
)
[
]
rnF
r
é um funcional universal da densidade de carga,
(
)
rn
r
. De acordo com o
princípio variacional, a energia do estado fundamental é o valor mínimo para a
densidade. Desta forma, reduzimos o problema de N-corpos exatamente à determinação
de uma função tridimensional
(
)
rn
r
que minimiza um funcional
(
)
[
]
rnE
r
. Apesar do
grande avanço obtido pelo teorema de Hohenberg e Kohn, ele não oferece uma maneira
prática para o cálculo da densidade de carga do estado fundamental de um sistema.
Um ano mais tarde, Kohn e Sham reformularam o problema e desenvolveram
um método simples para a execução de cálculos através da DFT. Primeiramente, o
sistema de elétrons interagentes é mapeado sobre um sistema fictício de elétrons não-
11
interagentes tendo a mesma densidade de carga do estado fundamental
(
)
rn
r
. Isto é feito
pela introdução dos orbitais de Kohn-Sham
(
)
r
i
r
Ψ para N elétrons de modo que
( ) ( )
∑
=
Ψ=
N
i
i
rrn
1
2rr
. (2.4)
De acordo com a conservação da carga, os orbitais de Kohn-Sham devem obedecer a
condição de ortonormalização, ou seja:
(
)
(
)
(
)
ijji
rdrr δ=ΨΨ
∫
r
r
r
*
. (2.5)
Assim, o funcional energia E pode então ser reescrito como:
( )
[ ]
(
)
(
)
( )
[ ]
( ) ( )
rdrvrnrnErdrd
rr
rnrn
rnTE
XC
rrrrrr
rr
r
r
r
∫∫
++
−
+= '
'
'
2
1
0
, (2.6)
onde o primeiro termo é a energia cinética de um gás de elétrons não-interagentes e é
dada por
( )
[ ]
( ) ( )
( )
'2
1
*
0
2
1
rdrrrnT
i
N
i
i
rr
r
rr
Ψ∇Ψ−=
∑
∫
=
, (2.7)
O segundo termo é a energia eletrostática dos elétrons, ou energia de Hartree. O terceiro
termo é a energia de troca e correlação que contém, principalmente, a contribuição não-
clássica para as interações elétron-elétron.
Minimizando a equação (2.6) com vínculo de manter constante o número de
partículas ou seja,
(
)
Nrdrn =
∫
r
r
, obtemos:
( ) ( ) ( )
rrrv
iiiKS
rrr
r
Ψ=Ψ
+∇− ε
2
2
1
, (2.8)
12
onde
i
ε são os multiplicadores de Lagrange e
(
)
rv
KS
r
é o potencial efetivo de Kohn-
Sham, dado por
( ) ( )
(
)
( )
rvrd
rr
rn
rvrv
XCKS
rr
rr
r
rr
∫
+
−
+= '
'
'
, (2.9)
com
( )
(
)
[
]
( )
r
rnE
rv
XC
XC
r
r
r
δ
δ
= . (2.10)
A equação (2.8) é denominada equação de Kohn-Sham e somente pode ser
resolvida através de um procedimento auto-consistente.
2.2.1 Aproximação da Densidade Local
A LDA [41] é usada para a determinação da energia de troca e correlação. Nesta
aproximação, supõe-se que a energia de troca e correlação por elétron no ponto
r
r
,
(
)
r
XC
v
ε seja igual à energia de troca e correlação por elétron em um gás homogêneo de
elétrons que tenha a mesma densidade
(
)
rn
r
em
r
r
. Dessa forma, a energia
[
]
)(rnE
XC
r
é
definida por
(
)
[
]
(
)
(
)
(
)
rdrnrnrnE
XC
r
r
r
r
∫
= ε ,
( )
( )( ) ( )
(
)
( )
rnn
XCXC
XC
dn
nd
nnrn
rn
E
r
r
r
=
−=≅
ε
εµ
δ
δ
, (2.11)
com
(
)
(
)
[
]
rnr
XCXC
r
r
hom
εε = .
Existem várias parametrizações para
(
)
[
]
rn
XC
r
hom
ε . A que é mais usada, é aquela
dada por Perdew e Zunger [42]. Esta parametrização é baseada em cálculos Monte
13
Carlo Quântico de Ceperley e Alder [43] de um gás de elétrons homogêneo em
diferentes densidades.
A LDA tem se mostrado bastante útil em cálculos de sólidos, apresentando bons
resultados para as suas propriedades estruturais e eletrônicas. Em geral, os cálculos que
utilizam a aproximação LDA superestimam a energia de ligação entre átomos e com
isto, subestimam os valores para as distâncias de ligação. Como conseqüência, nos
cálculos de estrutura de bandas, o valor do “gap” de energia normalmente apresenta um
valor menor que o experimental. Apesar disso, a forma e a dispersão dos níveis
eletrônicos para a banda de condução e a de valência se mostram muito parecidos com
os valores experimentais [44, 45].
2.2.2 Método das Ondas Planas
O método mais utilizado para o cálculo da energia total de um sólido consiste na
expansão das funções de onda de Kohn-Sham em uma base de ondas planas. Desta
forma, de acordo com o teorema de Bloch [46], a função de onda de um elétron em um
potencial periódico é escrita como
(
)
rki
j
kj
eru
r
r
r
r
.
,
=Ψ , (2.12)
onde j refere-se à banda e k
r
é um vetor de onda da primeira Zona de Brillouin. Uma vez
que
(
)
ru
j
r
é uma função periódica, podemos expandi-la em termos de um conjunto
discreto de bases de ondas planas ou seja,
∑
=
G
rGi
Gj
j
ecu
v
r
v
v
.
,
, (2.13)
na qual G
v
é o vetor da rede recíproca definido por mRG π2. =
r
v
, sendo m um inteiro,
R
r
é
um vetor de translação da rede cristalina e
Gj
c
r
,
é o coeficiente da expansão em ondas
14
planas. A função de onda eletrônica pode consequentemente, ser descrita como uma
combinação linear de ondas planas:
∑
+
+
=Ψ
G
rGki
Gkj
kj
ec
r
r
r
r
rrv
).(
,
,
. (2.14)
O uso de um conjunto de bases de ondas planas na expansão das funções de
ondas eletrônicas, conduz a uma formulação particularmente simples para as equações
de Kohn-Sham.
Em princípio, é necessário um conjunto infinito de ondas planas para a expansão
da função de onda eletrônica. Entretanto, os coeficientes relacionados a ondas planas
com energia cinética mais baixa são tipicamente mais importantes que os coeficientes
associados a ondas planas de alta energia cinética. Podemos então introduzir uma
energia cinética de corte, E
cut
, a fim de obtermos um conjunto de base finito. A energia
cinética de corte é definida como
2
2
2
Gk
m
E
cut
rr
h
+= . (2.15)
Logo, na expansão por ondas planas, as equações de Kohn-Sham podem ser
escritas como
(
)
(
)
(
)
Gkj
i
Gkj
G
HXCion
GG
ccGGVGGVGGVGk
rrrr
r
rr
rrrrrrrr
++
=⋅
−+−+−++
∑
,,
'
'
2
'''
2
1
εδ , (2.16)
que foram obtidas pela substituição da equação (2.14) na equação (2.8). De acordo com
essa equação, a representação do espaço recíproco da energia cinética é diagonal,
enquanto os potenciais são descritos em termos de componentes de Fourier. Em
princípio, esta equação secular pode ser resolvida pela simples diagonalização da matriz
hamiltoniana e o tamanho da matriz é determinado pela escolha da energia de corte,
E
cut
.
15
2.2.3 Método do Pseudopotencial
O método do pseudopotencial [41.] tornou-se bastante difundido nos cálculos das
propriedades estruturais e eletrônicas dos semicondutores. Este método atua na região
do caroço iônico (constituído pelos elétrons fortemente ligados), eliminando a
necessidade de incluir os estados do caroço. Desta forma, a função de onda nessa região
é substituída por uma função suave e sem nós, denominada pseudofunção de onda, a
qual pode ser descrita por um conjunto de ondas planas muito menor do que o conjunto
de ondas planas usado para descrever as funções de onda do sistema, devido à
existência de nós na região próxima ao núcleo. Embora os potenciais na região do
caroço possam não ser descritos corretamente, as funções de onda são corretamente
descritas na região de valência, como podemos observar na figura 2.1, que mostra o
conceito do pseudopotencial. Nela, as linhas sólidas mostram a função de onda de todos
os elétrons,
(
)
r
AE
r
Ψ , e o potencial iônico,
(
)
rv
AE
ion
r
, enquanto as linhas pontilhadas
mostram a pseudofunção de onda correspondente
(
)
r
PS
r
Ψ , dada pelo pseudopotencial,
(
)
rv
PS
ion
r
. Todas as quantidades são mostradas como uma função da distância, r, a partir
do núcleo atômico. O raio de corte r
c
marca o ponto além do qual as funções de onda de
todos os elétrons e as pseudofunções de onda tornam-se idênticas.
Figura 2.1: Ilustração esquemática do conceito de pseudopotencial.
Neste método, um potencial iônico forte
(
)
rv
ion
r
na região do caroço é
substituído por um pseudopotencial mais fraco
(
)
rv
PS
ion
r
. O conjunto correspondente de
pseudofunções de onda
(
)
r
PS
r
Ψ e as funções de onda de todos os elétrons
(
)
r
AE
r
Ψ são
idênticas para distâncias maiores que um raio de corte r
c
, escolhido de tal maneira que
16
represente o caroço. Com esta substituição,
(
)
r
PS
r
Ψ não possui as estruturas nodais que
causam as oscilações para distâncias menores que r
c
. Assim, a expansão em ondas
planas passa a ter um número razoável de termos.
A maioria dos pseudopotenciais usados em cálculos DFT são gerados a partir de
cálculos atômicos de todos os elétrons através da solução auto-consistente da equação
radial de Schrödinger,
( ) ( ) ( ) ( )
rrrvrv
r
Z
r
ll
dr
d
AE
ll
AE
lXCH
Ψ=Ψ
++−
+
+− ε
22
2
2
)1(
2
1
, (2.17)
onde
(
)
rv
H
é o potencial de Hartree,
(
)
rv
XC
é o potencial de troca e correlação e
AE
l
Ψ é
a função de onda de todos os elétrons com componente do momento angular l.
Convencionalmente, os pseudopotenciais são, então, construídos obedecendo a vários
critérios. São eles: i) os autovalores reais e os pseudo-autovalores devem concordar
entre si para uma determinada configuração escolhida, ii) as funções de onda reais e as
pseudofunções de onda devem concordar entre si dentro de uma distância escolhida
(raio de corte), iii) as integrais de 0 a r das densidades de carga real e das pseudofunções
de carga devem concordar entre si, para valores de r maiores que o raio de corte, iv) as
derivadas logarítmicas das funções de onda reais e das pseudofunções de onda,
juntamente com as suas primeiras derivadas com relação à energia devem concordar
entre si para valores de r maiores que o raio de corte.
Os pseudopotenciais que atendem a essas quatro propriedades são ditos
pseudopotenciais de primeiros princípios e de norma conservada [42].
Neste trabalho adotaremos a forma para o pseudopotencial elaborado por
Troullier e Martins [47], que será descrita no próximo item.
17
2.2.3.1 Pseudopotencial de Troullier-Martins
Troullier e Martins [47] em 1991, apresentaram um método diferente para a
construção das pseudofunções de onda. De acordo com este método, as pseudofunções
de onda podem ser definidas como
( )
(
)
( )
=
rpl
AE
nlPP
l
er
rR
rR
cl
cl
rr
rr
<
>
se,
se,
, (2.18)
onde )(rR
l
é a parte radial da função de onda com momento angular l. Os termos PP e
AE representam as pseudofunções de onda e as funções de onda de todos os elétrons,
respectivamente. O índice n nas funções de onda com todos os elétrons indicam o nível
de valência. A distância, além da qual as funções de onda com todos os elétrons e as
pseudofunções de onda são iguais,
cl
r , é também dependente de l. O polinômio, p(r), é
dado por
12
12
10
10
8
8
6
6
4
4
2
20
)( rcrcrcrcrcrccrp ++++++= (2.19)
Os coeficientes de p(r) são ajustados impondo a conservação da norma, a continuidade
das pseudofunções de onda e de suas quatro primeiras derivadas em
cl
rr = .
A grande vantagem desse método é a diminuição do tamanho da base de ondas
planas necessária para a convergência da energia total.
18
2.2.4 Método do Gradiente Conjugado
O gradiente conjugado [41] é um método usado na obtenção do ponto de mínimo
de uma função geral, quando conhecemos também o seu gradiente. Consiste, no entanto,
em uma técnica eficiente para a determinação do mínimo do funcional energia de Kohn-
Sham em função das posições dos átomos, uma vez que podemos calcular as forças
entre os átomos envolvidos.
Nesse método, o processo de busca do ponto de mínimo se inicia em um ponto
arbitrário, ou seja, o ponto 1 da figura 2.2, e depois calcula-se o gradiente da função.
Figura 2.2: Ilustração esquemática da busca do mínimo usando a aproximação do gradiente conjugado.
A direção inicial é dada pelo negativo do gradiente no ponto inicial. A nova
direção de busca contém informação, não somente do gradiente naquele ponto, mas
também do gradiente anterior que é da forma
1−
+=
mmmm
dgd
r
r
r
γ
, (2.20)
onde d
m
é a direção de busca no passo m e
m
g
r
é a direção do declive máximo (steepest
descent) no ponto m
(
)
m
xx
m
x
xF
g
=
∂
∂
−=
(2.21)
19
O coeficiente γ
m
é definido por
(
)
( )
11 −−
⋅
⋅
=
mm
mm
m
gg
gg
rr
r
r
γ (2.22)
com a condição inicial 0
1
=γ .
Stick e colaboradores [48] usaram um método do gradiente conjugado com uma
aproximação indireta para a determinação dos auto-estados do hamiltoniano de Kohn-
Sham. De acordo com esse método, na solução das equações de Kohn-Sham, a função F
a ser minimizada é o funcional energia de Kohn-Sham. As funções de onda {Ψ
i
}
desempenham o papel da variável x. O operador hamiltoniano H define o operador
gradiente e consequentemente a direção do declive máximo g. Inicia-se o procedimento
então, usando uma função de onda tentativa, calculando o seu gradiente conjugado. As
posições dos átomos são alteradas e as novas funções de onda são obtidas com o seu
gradiente conjugado. O procedimento é repetido até que se obtenha a função de onda
que minimiza o funcional Kohn-Sham. A ortonormalidade das funções de onda, que
deve ser mantida, é alcançada modificando adequadamente a direção do declive
máximo.
É interessante ressaltar que a implementação do método do gradiente conjugado
deve além de aumentar a velocidade computacional, reduzir a memória exigida nos
cálculos a fim de que os cálculos não sejam limitados pela memória computacional
disponível.
2.2.5 Amostragem da Zona de Brillouin
Baseado no teorema de Bloch, qualquer integral do espaço real sobre um sistema
periódico com uma extensão infinita, pode ser substituída por uma integral no espaço
recíproco sobre a primeira zona de Brillouin (finita). Contudo, isto ainda requer cálculos
de funções periódicas em um número infinito de pontos no espaço recíproco,
denominados pontos k
r
. Este problema pode ser superado em sistemas de alta simetria,
20
através de somatórios em apenas uma parte da zona de Brillouin como descrito por
Baldereschi [49]. Posteriormente, Chadi e Cohen [50], Monkhorst e Pack [51], baseados
nestes argumentos de simetria, descreveram métodos sistemáticos para gerar conjuntos
de pontos especiais que representam adequadamente a zona de Brillouin. Dessa forma,
qualquer função integrada
(
)
rf
r
, tais como a densidade ou a energia total, pode ser
calculada como uma soma discreta, e ponderada
(
)
( )
j
j
j
BZ
kFwkdkF
∑
∫
Ω
=
1
r
r
, (2.23)
em que
(
)
kF
r
é a transformada de Fourier de
(
)
rf
r
, Ω é o volume da célula e w
j
são os
fatores pesos.
As posições dos pontos k
r
dentro da zona de Brillouin devem ser
cuidadosamente selecionadas uma vez que uma escolha prudente resultará em uma
descrição eficiente de um sistema particular, conduzindo a uma redução significativa do
tempo computacional.
Neste trabalho, os nossos cálculos foram baseados no método de Monkhorst-
Pack [51]. De acordo com este método, os pontos k
r
são distribuídos homogeneamente
através do espaço em linhas e colunas que seguem a forma da zona de Brillouin ou seja,
332211
bxbxbxk
jjjj
r
r
r
r
++= , (2.24)
onde
1
b
r
,
2
b
r
,
3
b
r
são vetores da rede recíproca, e
j
i
ij
n
l
x = , j = 1, ..., n
j
(2.25)
onde l
i
são os comprimentos dos vetores da rede recíproca e n
j
caracteriza o número de
pontos especiais no conjunto.
O número de pontos k
r
pode variar dependendo do tamanho da supercélula no
espaço real. Para uma célula grande, um único ponto k
r
pode ser suficiente para
21
descrever o sistema. Para uma célula pequena, um número maior de pontos k
r
deve ser
usado.
2.2.6 Forças de Hellmann-Feynman
Uma importante conseqüência do caráter variacional da DFT é a possibilidade
de calcular as forças de Hellmann-Feynman que atuam sobre os átomos [52, 53]. Estas
forças são definidas como a derivada da energia total em relação às posições atômicas
i
R
r
, ou seja
( ) ( )
( )
( )
( ) ( )
rdrnrV
rn
E
RErdrVrnEF
iiii
R
II
RRR
i
rr
r
r
r
rr
rr
r
r
rr
rrrr
∇
+−∇−∇−=∇−=
∫∫
δ
δ
, (2.26)
onde E
II
é a energia de interação íon-íon.
Existem dois termos que contribuem para o gradiente da energia total, um
proveniente da dependência explícita sobre as posições atômicas através do potencial
(
)
{ }
(
)
rVrV
R
r
r
r
≡ , e o outro proveniente de uma dependência implícita através da
densidade de carga do estado fundamental
(
)
rn
r
. A derivada explícita e os termos de
interação íon-íon não são difíceis de serem calculadas. O mesmo não pode ser dito para
a derivada implícita. Na verdade, precisamos conhecer como a densidade de carga varia
em função da variação da posição dos átomos. O último termo na equação (2.26)
contém a variação de primeira ordem do funcional energia em torno da energia do
estado fundamental, que se anula quando o estado fundamental é um mínimo. Em
termos matemáticos, a minimização da energia com a restrição de manter a carga total
constante implica que a variação da função
(
)
[
]
(
)
(
)
(
)
(
)
{
}
∫∫
−−+ NrdrnrdrVrnrnf
r
r
r
r
r
r
λmin (2.27)
22
em relação a uma variação arbitrária
(
)
rn
r
δ , deva anular. Nesta função, λ é um
multiplicador de Lagrange e N é o número total de elétrons.
Usando cálculos variacionais elementares, temos que:
( )
( ) ( )
0=
−+
∫
rdrnrV
rn
F rrr
r
δλ
δ
δ
, (2.28)
ou ainda,
( )
( )
λ
δ
=
+ rV
rdn
F r
r
(2.29)
Logo, forças são simplesmente um elemento de matriz do estado fundamental do
gradiente do potencial externo adicionado a um termo de interação íon-íon
(
)
(
)
(
)
RErdrVrnF
II
RR
i
ii
r
r
r
r
r
r
r
rr
∇−∇−=
∫
. (2.30)
2.3 Método de Monte Carlo
A expressão Monte Carlo é usada para métodos de simulação computacional
baseados na utilização de números aleatórios [54]. O nome Monte Carlo surge do
caráter aleatório ou probabilístico do método, e é considerado uma alusão à cidade de
Monte Carlo no principado de Mônaco, famosa pelos seus cassinos e jogos de azar.
Este método é utilizado para simular flutuações estatísticas a fim de gerar
numericamente as distribuições de probabilidades, e foi primeiro usado por Metropolis e
colaboradores [55] para estudar um sistema de esferas duras. Esta técnica é geralmente
usada no estudo das propriedades de equilíbrio de um sistema de átomos. Para o estudo
de processos de não-equilíbrio ou de relaxação, é usado o método de Monte Carlo
Cinético.
Nesta seção, apresentamos uma breve descrição do método de Monte Carlo e o
algoritmo de Metropolis antes de descrevermos o método de Monte Carlo Cinético.
23
Posteriormente, são apresentados os modelos de crescimento cristalino utilizados neste
trabalho.
2.3.1 O Método de Monte Carlo e o algoritmo de
Metropolis
A essência da simulação Monte Carlo consiste na obtenção das médias de
qualquer propriedade através de técnicas de amostragem, podendo ser usado na solução
de problemas determinísticos. Neste método, o problema é tratado a partir de um
conjunto finito de pontos ou eventos, que seja estatisticamente representativo (uma
amostragem). Como exemplo, podemos citar a determinação do número π, em que um
problema determinístico é substituído por um problema estocástico envolvendo
probabilidades da razão entre a área de um quadrado e o seu maior círculo inscrito.
A geração de números aleatórios é, então, uma característica importante desse
método. Contudo, em termos computacionais, é impossível obter um gerador de
números verdadeiramente aleatório. Torna-se necessário, o uso de números pseudo-
aleatórios. O termo pseudo refere-se ao fato de que um determinado valor é usado para
inicializar o processo de geração.
Métodos de Monte Carlo são utilizados como ferramentas computacionais em
muitas áreas da física. Suas aplicações vão desde a mecânica estatística à física de
partículas elementares. O marco histórico no desenvolvimento deste método é o
trabalho de Metropolis e colaboradores, em 1953 [55]. Para descrevermos a idéia
principal do algoritmo, faremos uma breve revisão da fundamentação teórica do método
de Monte Carlo [34].
Consideremos que o estado de um determinado sistema S = {C} consista de
todas as possíveis configurações, C, e que a distribuição de probabilidades das
configurações no tempo t seja P (C). A evolução temporal de P(C,t) é feita pela
Equação Mestra, que é dada por
(
)
( ) ( ) ( ) ( )
∑∑
→+→=
∂
∂
''
,'','
,
CC
TCPCCWTCPCCW
t
TCP
, (2.31)
24
onde
(
)
'CCW → é a taxa de transição do estado C para o estado C’. No equilíbrio,
obtemos uma descrição única para P(C) que satisfaça a seguinte condição estacionária:
(
)
0, =∂∂ ttCP . Na prática, isto pode ser substituído pela condição denominada
“balanço detalhado”, que é
(
)
(
)
(
)
(
)
TCPCCWTCPCCW ,'',' →=→ . (2.32)
Num ensemble canônico, a distribuição de probabilidades é dada por:
( )
T
k
E
B
C
e
Z
tCP
−
=
1
, , (2.33)
em que E
C
é a energia da configuração C, T é a temperatura, k
B
é a constante de
Boltzmann e Z é a função partição. A probabilidade P(C,t) é geralmente desconhecida
devido à função partição Z, que aparece no denominador. No entanto, podemos
estabelecer uma cadeia de Markov pela geração de uma seqüência de configurações
convergindo para a distribuição de equilíbrio P(C) desejada. Calculando a razão
(
)
(
)
tCPtCP ,,' , notamos que as funções partição se cancelam, dando origem à seguinte
forma para a condição do “balanço detalhado”:
(
)
( )
∆−
=
→
→
T
k
E
B
e
CCW
CCW
'
'
, (2.34)
onde ∆E = E
C’
– E
C
. De acordo com essa equação, a razão das probabilidades de
transição depende somente da diferença de energia entre os estados C e C’. Qualquer
taxa de transição que satisfaça à condição do “balanço detalhado” é aceitável. A escolha
mais comum foi introduzida por Metropolis e colaboradores [55], que é dada por:
( )
=→
−
∆−
−
1
0
1
0
'
τ
τ
T
k
E
B
e
CCW
casosoutrospara
Epara
,
0,
>
∆
(2.35)
25
Nesta equação τ
0
refere-se ao tempo entre movimentos de tentativas (freqüentemente é
colocado igual a um). Em alguns sistemas, o parâmetro τ
0
pode ser colocado igual a
alguma escala de tempo característica do sistema e, consequentemente, o tempo de
Monte Carlo pode ser interpretado como uma medida do “tempo físico real”. Esta
interpretação torna-se bem justificada através de exemplos de simulações de processos
que são aleatórios por natureza, como a difusão.
Considerando um sistema de N partículas interagentes, as interações entre as
partículas são descritas por um potencial )(
N
rU , onde r
N
denotam as coordenadas de
todas as partículas. Então, o algoritmo de Metropolis pode ser realizado de acordo com
o diagrama descrito na figura 2.3.
Escolha um estado inicial
Escolha aleatoriamente uma partícula, i
Calcule um deslocamento aleatório para
a
partícula,
δ+=
ii
rr
v
v
'
Calcule a variação da energia
correspondente,
( ) ( )
NN
rUrUU
rr
−=∆
'
∆
U < 0 ?
Sim
Gere um número
aleatório r tal que
0 < r <1
Não
∆−
<
T
k
U
B
er
Aceite o
movimento
Mantenha a
configuração
antiga
Sim
Não
Figura 2.3: Fluxograma para o algoritmo de Metropolis.
Torna-se necessário ressaltar que neste algoritmo, o tempo é medido através de
passos de Monte Carlo, ou seja, através da contagem do número de configurações
tentativas geradas no processo.
?
26
2.3.2 O Método de Monte Carlo Cinético
O método Monte Carlo Cinético é adequado para o estudo de fenômenos
dependentes do tempo, tais como a difusão e os processos de crescimento de cristais, ou
seja, fenômenos cinéticos [56].
Neste método, as mudanças das configurações correspondem a eventos reais.
Dessa forma, é possível definir um tempo real, que se baseia em tempos característicos
do sistema. Vários algoritmos foram feitos para este método [56], possibilitando a
otimização computacional das simulações. Para descrevermos como um destes
algoritmos para esse método é construído, utilizaremos como exemplo, o crescimento
cristalino. Supondo que N, seja o número de eventos possíveis ocorrendo com uma certa
probabilidade por unidade de tempo, R
α
(α = 1, ... , N), numa dada configuração C, em
que N e o conjunto {R
α
} são dependentes da configuração C, podemos, então, definir a
taxa total das transições como:
( )
∑
=
==
N
RCQQ
1α
α
(2.36)
As probabilidades de transição entre as configurações C e C’ são formalmente
escritas como:
( ) ( )
∑
=
→=→
N
CCRCCW
1
'V'
α
α
α
(2.37)
onde
(
)
'V CC →
α
é uma matriz estocástica que especifica se a transição C → C’ pode
ser realizada pelo evento α, ou não. A probabilidade de ocorrência do evento α é
R
α
/Q(C), ou seja, os eventos são selecionados com probabilidades proporcionais às suas
taxas de ocorrência.
Um dos algoritmos para o Método de Monte Carlo Cinético, freqüentemente
usado nas simulações do crescimento cristalino, é o algoritmo “N-fold way” proposto
por Bortz, Kalos e Lebowitz, também denominado algoritmo BKL [57], que se baseia
no Modelo de Ising. Neste tipo de algoritmo, o processo de proposição e teste de
27
configurações é substituído por uma lista contendo todos os eventos (ou transições)
possíveis, ponderados pelas suas respectivas probabilidades. Depois de cada evento, é
feita uma atualização da lista de eventos, pois novas transições são possíveis, assim
como outras transições são impossibilitadas de ocorrerem. Este procedimento resulta no
aumento do desempenho computacional, principalmente nos estágios próximos ao
equilíbrio, na qual a “taxa de rejeição” de novos eventos passa a ser muito alta.
De uma maneira simplificada, a seqüência dos passos executados pelo algoritmo
BKL está descrita na figura 2.4.
Inicialize osistema no tempo t = 0
“
loop
” dependente do tempo ou do número de
passos Monte Carlo Cinético da simulação
Forme uma lista de todas as taxas ri de todas as
transições possíveis w
i
no sistema
Calcule a função cumulativa
∑
=
=
i
j
ji
rR
1
, para i = 1, ..., N
Denote R = R
N
e obtenha umnúmero aleatório
u ª [0,1]
Execute o evento i, obtenha todas w
i
e recalcule todos
r
i
que podemter sofrido variações devido à transição
Obtenha um novo número aleatório
u’ ª [0,1]
Atualize o tempo com t = t +∆t,onde
∆
t = - log u’/R
e o cálculo do passo Monte Carlo Cinético
Final do “loop”? Finalize
SimNão
Figura 2.4: Fluxograma da seqüência de passos executados pelo algoritmo BKL.
28
2.3.3 Modelos de Crescimento Cristalino
Simulações de Monte Carlo Cinético para o crescimento de cristais são baseadas
em modelos simplificados deste fenômeno [56]. Para a construção destes modelos, a
estratégia usada consiste em considerar apenas os processos provavelmente mais
importantes, desconsiderando os detalhes menos significativos. Dessa forma, algumas
aproximações devem ser utilizadas durante as simulações. Neste trabalho, as nossas
simulações foram baseadas no modelo “Solid-on-Solid” (SOS), que fornece uma
descrição qualitativa do crescimento por MBE, que se baseia no modelo do gás de rede.
Neste modelo, um número grande de partículas (átomos ou moléculas) se move ao redor
dos vértices de uma rede que representa o substrato [34]. A superfície, então, é
representada por uma rede bidimensional de N x N sítios, cuja geometria corresponde à
sua estrutura de relaxação ou reconstrução (quadrada, retangular, triangular, hexagonal,
etc.). As partículas, por sua vez, podem ocupar somente posições discretas sobre esta
rede. O sistema é, deste modo, descrito por uma rede com alguns sítios ocupados por
partículas e alguns sítios vazios. As regiões de alta concentração de átomos
correspondem ao sólido e as regiões de baixa concentração de átomos, ao vapor,
conforme descrito na figura 2.5.
Figura 2.5: Configuração do gás de rede com a região inferior correspondendo ao sólido e a região
superior ao gás. Átomos são representados por quadrados hachurados.
Na figura 2.5, somente algumas propriedades estruturais foram levadas em
consideração e muitos dos aspectos dos cristais reais foram desprezados. Os átomos são
tratados como objetos clássicos estáticos e, portanto, as propriedades mecânicas,
químicas, e outras foram desconsideradas.
Outras aproximações, além das utilizadas em um modelo de gás de rede, podem
também ser usadas. Podemos, por exemplo, supor que os processos dentro de um fluido
29
podem ser desprezados e que o crescimento possa ser bem descrito apenas por
processos na interface entre o fluido e o sólido (Figura 2.6).
Figura 2.6: Posições dos átomos do sólido são considerados somente, novos átomos são adicionados
através de regras de crescimento definidas no modelo.
No mesmo nível de aproximações, podemos também desprezar a difusão dos
átomos dentro do cristal. Na figura 2.6, por exemplo, estamos interessados apenas nos
movimentos dos átomos representando o sólido e o crescimento ocorre pela adição de
átomos nas posições determinadas pelas regras de crescimento. Esta aproximação é bem
justificada pelo crescimento balístico, quando as interações mútuas de partículas num
fluido são desconsideradas. Um importante exemplo é o crescimento por MBE, na qual
a superfície é diretamente bombardeada por átomos. A omissão dos processos dentro
das fases pode ser problemática para o crescimento por vapor próximo ao equilíbrio,
uma vez que as proximidades da interface podem sofrer variações devido, por exemplo,
à influência da temperatura e da pressão [58]. Átomos dentro do fluido estão se
difundindo através da interface onde o crescimento ocorre. Novamente, podemos
esperar que, usando regras de crescimento apropriadas, um bom modelo qualitativo
pode ser formulado, mas a descrição correta dos processos dentro das fases deve ser
importante somente em alguns casos.
No caso do crescimento com interface planar, e sob condições onde há muito
poucas vacâncias e/ou saliências (Figura 2.6) no material em crescimento, podemos usar
a aproximação SOS (Figura 2.7).
Figura 2.7: Configurações superficiais na aproximação SOS.
30
De acordo com essa aproximação cada átomo é colocado sobre um outro átomo.
A superfície é, então, descrita por uma única função, h(x), da coordenada x do substrato.
Em princípio, todas as estruturas de rede relevantes podem ser implementadas em um
modelo SOS. O nível adequado, ou necessário, de sofisticação dependerá do problema a
ser investigado.
Além dos aspectos geométricos dos modelos de crescimento, devemos também
considerar a definição das regras de crescimento dinâmico, que também exige algumas
aproximações. É geralmente suposto que o movimento de partículas individuais (átomos
ou moléculas) ocorre instantaneamente e que os movimentos são independentes e
Markovianos. As partículas podem se mover de um sítio para outro com uma freqüência
proporcional às taxas dos processos. As taxas de processos elementares não são muito
bem conhecidas, e são deduzidas com o uso de todas as informações teóricas e
experimentais disponíveis.
2.3.4 Simulações Monte Carlo Cinético para o
crescimento por MBE.
As simulações de crescimento de materiais sobre superfícies semicondutoras
ainda são um grande desafio para a comunidade científica, uma vez que a maioria destas
superfícies sofre reconstruções e, consequentemente, podem apresentar vários sítios
para a adsorção de átomos [59]. Argumentos simples baseados na contagem do número
de ligações de um átomo adsorvido e seus vizinhos, apesar de terem alguma justificativa
física para a adsorção de um átomo sobre uma superfície metálica, não são aplicadas a
superfícies semicondutoras. Nas superfícies reconstruídas dos semicondutores do grupo
IV e III-V, com a estrutura cristalina do diamante ou zinc-blende, os átomos da
superfície, na maioria das vezes, têm três ligações, e os átomos adsorvidos podem ter
duas ou três ligações. Por esta razão, a contagem de ligações não é suficiente para
classificar os diferentes sítios de ligação e as barreiras de difusão de um átomo
31
adsorvido sobre a superfície. Então, a aproximação utilizada neste trabalho explora o
fato de que podemos representar as configurações atômicas que ocorrem durante o
crescimento sobre uma rede cristalina. Substituímos, então, a dinâmica newtoniana dos
átomos da superfície pela dinâmica discreta. Os átomos podem se mover de um sítio da
rede para outro, ou aparecer e desaparecer sobre a superfície (adsorção e dessorção)
através de saltos discretos. Após cada salto, o sistema entra em equilíbrio na nova
geometria antes que uma transição para outro estado ocorra devido a uma flutuação
térmica. Deste modo, nas temperaturas usadas durante a deposição do filme, estes saltos
correspondem a eventos raros sobre as escalas de tempo das vibrações atômicas.
Concentrando-nos sobre os saltos, eliminamos a informação superficial sobre a
dinâmica dos átomos entre estes eventos raros, ou seja, as vibrações dos estados
discretos e as suas escalas de tempo associadas. Não desprezamos esta informação, mas
a descrevemos em termos de um prefator probabilidade. Podemos demonstrar isto
através de um exemplo simples: o passeio de uma única partícula sobre uma superfície
[60]. A nível atômico, o movimento das partículas é governado por uma superfície de
energia potencial (PES), que é a energia potencial experimentada pela difusão do átomo
adsorvido na superfície, ou seja,
( )
{ }
=
→
I
adadad
tot
RZ
adad
PES
RZYXEYXE
Iad
,,,min,
,
, (2.38)
onde
→
I
adadad
tot
RZYXE ,,, é a energia do estado fundamental do sistema
“superfície mais átomo adsorvido” (também denominada energia total) na configuração
atômica
→
I
adadad
RZYX ,,, . A PES representa a energia total mínima em relação à
coordenada z do átomo adsorvido, Z
ad
, e todas as coordenadas dos átomos do substrato,
→
I
R . Se desconsiderarmos, por um momento, as contribuições vibracionais para a
energia livre, o mínimo da energia potencial superficial representa os sítios estáveis e
meta-estáveis dos átomos adsorvidos.
32
A descrição cinética do crescimento define os processos atomísticos básicos, tais
como difusão ou dessorção, em termos das probabilidades com que eles ocorrem. Para a
taxa,
J
Γ , de um processo molecular, J, (difusão, por exemplo), podemos escrever,
( )
TkF
B
J
BJ
e
h
Tk
∆−
=Γ , (2.39)
onde
J
F∆ é a diferença na energia livre de Helmholtz entre o máximo (ponto de sela) e
o mínimo (geometria inicial) da PES ao longo da trajetória de reação do processo J, T é
a temperatura, k
B
é a constante de Boltzmann e h é a constante de Planck. A equação
seguinte fornece a energia livre de ativação,
J
F∆ , que o sistema precisa para mover de
uma posição inicial para o ponto de sela:
vib
JJJ
STUF ∆−∆=∆ . (2.40)
A energia interna, U, nesta expressão, consiste da energia total (estática), E
J
, que
é obtida a partir de cálculos DFT, e a energia vibracional, U
vib
.
J
U∆ é a diferença da
energia interna do sistema com a partícula no ponto sela e no ponto de mínimo.
vib
J
S∆ é
a diferença análoga na entropia vibracional. Freqüentemente, a taxa de um processo, J,
pode ser descrita pela lei de Arrhenius,
(
)
TkE
JJ
BJ
e
∆−
Γ=Γ
0
, (2.41)
onde
( )
TkUkS
B
J
B
vib
JB
vib
J
e
h
Tk
∆−∆
=Γ
0
(2.42)
é a freqüência tentativa. Para modos vibracionais harmônicos clássicos, podemos
mostrar que a freqüência tentativa torna-se independente da temperatura. Na equação
2.41, tanto
0
J
Γ quanto
J
E∆ podem ser obtidos de cálculos DFT. Para cálculos de difusão
de uma única partícula, podemos ler
J
E∆ diretamente da PES. Neste caso, encontramos
que
0
J
Γ é da mesma ordem, aproximadamente 10
13
s
-1
, para vários processos. Como
uma primeira aproximação, podemos então evitar recalcular
0
J
Γ para cada processo,
33
pois variações em
0
J
Γ por um fator 5 são tipicamente irrelevantes quando comparadas às
variações de
J
E∆ no expoente. Entretanto, quando muitos átomos se movem
coletivamente, o prefator pode ser significantemente diferente. Para processos
envolvendo várias partículas (nucleação, formação de ilhas, etc), devemos levar em
consideração as interações entre as partículas. Estas quantidades são acessíveis em
cálculos DFT, bem como pelo cálculo da energia potencial superficial ou parte dela,
para todas as situações que ocorrem durante o crescimento. Contudo, isto envolve a
realização de cálculos para muitos arranjos atômicos diferentes. A fim de descrevermos
a ligação química local de um átomo da superfície, é geralmente suficiente considerar o
ambiente local de seus vizinhos mais próximos, contanto que as condições de contorno
sejam propriamente consideradas. Para obtermos estas configurações das energias de
ligação e das barreiras, realizamos cálculos DFT, e os resultados são coletados. Depois,
um programa computacional, que realiza uma simulação do crescimento Monte Carlo
Cinético, lê as energias destas barreiras a partir de um arquivo, e calcula uma lista de
taxas de transição de acordo com a equação 2.41. Cada passo da simulação seleciona
aleatoriamente um processo particular com uma probabilidade proporcional à sua
contribuição para a taxa total. O tempo da simulação é ajustado por um incremento dado
pelo inverso da taxa total multiplicado por um número aleatório retirado a partir de um
decaimento exponencial da distribuição de probabilidade.
Deste modo, um mapeamento fisicamente significativo entre o tempo de
simulação e o tempo físico é estabelecido, uma vez que os processos microscópicos
podem ser tratados, durante um intervalo de tempo pequeno, como processos
independentes de Poisson [61]. No total, uma simulação típica de crescimento pode
compreender bilhões de tais passos de simulação.
Como uma maneira conveniente de representar os dados sobre as interações
atômicas obtidos a partir de cálculos DFT, é algumas vezes útil introduzir um
Hamiltoniano de gás de rede. Este conceito já foi utilizado para estudar as interações de
adsorvatos e seus efeitos sobre sua dessorção em superfícies [62]. Recentemente, foi
aplicado também no contexto da nucleação de ilhas e do crescimento cristalino [63].
Introduzindo um Hamiltoniano de gás de rede [63; 64] vamos convenientemente
representar os dados de interações atômicas obtidos de cálculos DFT. Neste trabalho
usaremos esta aproximação para o crescimento, por MBE, do GaN sobre a superfície
(100) de GaAs.
34
A fim de descrevermos as energias de ligação dos átomos de N nos vários sítios
da superfície reconstruída e a ligação dos átomos de Ga adsorvidos nos vários sítios,
propomos um Hamiltoniano de gás de rede que consiste de uma contribuição de um
único sítio, H
0
, adicionado a um termo de interação:
int0
HHH += . (2.43)
Seja k um sítio da rede e n
k
uma variável discreta que tem o valor 1 se um átomo
de N ocupa o sítio, ou 0 se o sítio está vazio. Analogamente, m
k
é uma variável de
ocupação que descreve os átomos de Ga adsorvidos no sítio. Então, H
0
tem a forma
∑
∑
∑
−
++=
k
kk
NGa
k
k
Ga
b
k
k
N
b
mnVmEnEH
0
. (2.44)
A superfície de GaAs (100) reconstruída apresenta algumas complexidades por
causa da necessidade de distinguir os vários sítios da rede não equivalentes devido às
reconstruções. Para esclarecer a idéia básica, deixaremos estas complicações de lado por
um momento. Variações das energias totais, devido às diferentes reconstruções e
diferentes arranjos possíveis de átomos de Ga, estão descritas na parte de interação do
Hamiltoniano de gás de rede,
...
3
1
2
1
2
1
2
1
2
1
'
,,
,
2
,
1
,
2
,
1int
'
++
+++
++=
+
++
+
+
+
+
+
+
+
+
+
∑
∑∑
∑∑
ak
akaak
akk
Ga
trio
bkk
bkk
Ga
n
akk
akk
Ga
n
bkk
bkk
As
n
akk
akk
As
n
mmmV
mmVmmV
nnVnnVH
, (2.45)
onde k+a e k+b indicam os sítios vizinhos e de próximos relativos a k. Os parâmetros
de interação
As
n
V
1
,
As
n
V
2
,
Ga
n
V
1
,
Ga
n
V
2
e
Ga
trio
V são determinados através de cálculos DFT
realizados para um conjunto de configurações atômicas. A adição de mais tipos de
interação refina a descrição do sistema através do Hamiltoniano de gás de rede.
Para iniciar uma simulação Monte Carlo Cinético, precisamos das energias de
ligação bem como das barreiras de energia para o cálculo das taxas. Estas barreiras
mostram diversidade considerável, pois em cada processo específico, a partir de um
35
estado inicial com um dado ambiente local através de uma transição para o estado final
com um ambiente local diferente, exige uma barreira específica. Devemos ter o cuidado
para não violarmos o princípio do “balanço detalhado”. Este princípio fundamental da
mecânica estatística exige que o quociente das taxas de quaisquer processos, obtidos
pela alternação dos estados iniciais e finais, seja igual a,
(
)
( )
( )
TkG
if
fi
b
e
∆−
=
Γ
Γ
, (2.46)
onde G
∆
é a diferença na energia livre de Gibbs no estado final e inicial. Se as
contribuições vibracionais para a entropia bem como os efeitos de pressão podem ser
desprezados, a diferença G
∆
reduz-se à diferença de energia de ligação,
(
)
(
)
1,00,1 ==+===∆≈∆
fifib
nnHnnHEG . A aproximação usando um
Hamiltoniano de gás de rede, além de ser considerada uma maneira conveniente de
armazenar e recuperar informação sobre as energias de ligação, oferece a vantagem de
que a relação do balanço detalhado (2.46) seja sempre satisfeita, se as energias de
ligação dos estados final e inicial são tomadas a partir do Hamiltoniano de gás de rede.
2.4 Metodologia
Neste trabalho, usamos a Teoria do Funcional da Densidade dentro da
aproximação local para a densidade, expansão da função de onda em ondas planas e o
método pseudopotencial, através do programa ABINIT [37]. Para a energia de troca e
correlação usamos a aproximação proposta por Ceperley e Alder [43] como
parametrizada por Perdew e Zunger [42]. O programa foi utilizado para a obtenção tanto
dos aspectos estruturais e eletrônicos das superfícies (100) de GaAs, na reconstrução
(2x1), terminadas em Ga e em As, como os aspectos energéticos da adsorção de
Nitrogênio sobre estas superfícies, que foram utilizados no Método de Monte Carlo
usado para simular o processo de Nitridação desta superfície através do MBE.
Utilizamos os pseudopotenciais de Troullier-Martins, que foram construídos com o
36
código fhi98PP [36]. Uma descrição detalhada do processo de preparação dos
pseudopotenciais, suas propriedades e convergência está descrita no Apêndice A. Antes
de realizar os cálculos para a “slab supercell”, calculamos a energia total para vários
valores do parâmetro de rede a fim de obtermos uma boa descrição das propriedades de
equilíbrio cristalino do GaAs, checando a eficiência do método e metodologia
escolhidos. Usamos, para este fim, a técnica do gradiente conjugado, enquanto que a
autoconsistência foi realizada através da técnica de Broyden. A partir do ponto de
mínimo obtido, calculamos a variação da energia total do bulk de GaAs em função da
variação do volume da célula unitária. Os resultados obtidos foram ajustados pela
equação de Murnaghan [65] e foram convergidos para uma energia máxima da
expansão de ondas planas de 140 ou 80 Ry, dependendo da inclusão ou não dos elétrons
3d do Ga, respectivamente. Foi usado um conjunto de 10 pontos k
r
especiais através de
uma rede de Monkhorst-Pack (4 4 2) para a amostragem da zona de Brillouin.
Em nossos cálculos, para a simulação da superfície (100) de GaAs na
reconstrução (2x1), usamos o modelo de lâminas [66], que consiste em montar uma
supercélula contendo átomos e vácuo que se repetem periodicamente no espaço. A
camada de vácuo deve ser colocada em uma das direções, tornando o cristal finito na
direção de crescimento. Assim, teremos nas outras direções, um cristal infinito e, desta
forma, a periodicidade na direção de crescimento é retirada, que é o nosso objetivo. O
número de átomos e sua distribuição na célula é que indicará o tipo de superfície a ser
estudada, se é (100), (110) ou (111) e, o tamanho do vácuo é tal que dificulte a interação
entre as duas superfícies da lâmina. A superfície (100) de GaAs foi, então, modelada
utilizando-se uma supercélula com cinco camadas atômicas e uma região de vácuo
equivalente a cinco camadas atômicas. Na figura 2.4.1, mostramos a lâmina (slab) usada
para estudar os aspectos estruturais da superfície (100) de GaAs, na reconstrução (2x1),
utilizando cálculos ab initio.
37
Figura 2.6: Modelo de lâmina (slab) que descreve a superfície (001) de GaAs.
A fim de determinarmos as estruturas atômica e eletrônica destas superfícies, a
superfície inferior foi passivada com átomos de Hidrogênio para evitar tanto as
transferências de cargas, como a indução de dipolos entre as superfícies. O passo inicial
de nossos cálculos consistiu na otimização da posição dos átomos de Hidrogênio que
saturam a superfície, sendo que, durante esse processo, todos os outros átomos da
superfície de GaAs foram mantidos fixos nas posições cristalinas. Posteriormente, as
três camadas superiores foram relaxadas, enquanto a outra camada de Ga e a camada de
Hidrogênio foram mantidas congeladas nas posições de equilíbrio, a fim de simular as
reconstruções observadas. Os resultados convergem para uma E
cut
na expansão de ondas
planas de 60 Hartrees (120 Ry). A energia total da célula foi minimizada e a posição dos
átomos foi relaxada até que a força sobre eles em todas as direções fosse menor que 1
mHartree/Bohr. Um conjunto de 8 pontos k
r
foi usado através de uma rede de
Monkhorst-Pack (4 4 2) para a amostragem da zona de Brillouin. Os efeitos dos orbitais
3d do Ga não são levados em consideração.
Para a simulação da adsorção do átomo de Nitrogênio sobre a superfície (100) de
GaAs, na reconstrução (2x1), fizemos o mapeamento da energia total de um átomo de N
adsorvido na superfície a fim de verificarmos qual será o sítio mais favorável para a
deposição desse átomo sobre o substrato de GaAs. Fixamos, então, tanto a posição (x,y)
do átomo adsorvido como a de todos os átomos da superfície, deixando apenas que o
átomo adsorvido relaxasse na direção z. Após encontrado o valor da altura ∆z do átomo
adsorvido, que corresponde a um valor de energia total mínima, as coordenadas (x,y) do
átomo adsorvido e as posições dos átomos das duas primeiras camadas da superfície
foram otimizadas. Obtivemos também o valor da energia de adsorção, E
ads
, através da
38
relação: E
ads
= E(N/GaAs) – E(GaAs) – E(N), onde E(N/GaAs) representa a energia total
da superfície (100) de GaAs, com um átomo de N adsorvido, E(GaAs) é a energia total
da superfície (100) de GaAs limpa e E(N ) é a energia total de um átomo de N.
Analisamos as variações nos comprimentos e nos ângulos de ligação, bem como
as energias totais envolvidas antes e após a deposição.
A fim de simularmos a difusão de um átomo de N sobre a superfície (100) de
GaAs, terminada em As, na reconstrução (2x1), calculamos a energia total da superfície
para um átomo de N adsorvido nesta superfície. Plotamos, então, a variação da energia
total como uma função do deslocamento de um átomo de N, iniciando a partir do sítio
de adsorção E, ao longo da direção [-1 1 0]. O mesmo procedimento foi utilizado ao
longo da direção [1 1 0] e sobre o dímero de As. Obtivemos, assim, os respectivos
valores para as barreiras de difusão.
Simulamos também a adsorção de um segundo átomo de N sobre a superfície
(100) de GaAs, na reconstrução (2x1). Nesta simulação utilizamos a mesma
metodologia que foi aplicada à adsorção de um átomo de N. Contudo, consideramos o
primeiro átomo de N adsorvido no sítio de adsorção de menor energia, ou seja, ao sítio
de adsorção mais estável.
Como resultado das análises dos cálculos obtidos, encontramos tanto os sítios de
adsorção metaestáveis para o átomo de N adsorvido, como as barreiras de energia de
difusão desse átomo sobre a superfície (100) de GaAs. Essas barreiras de energias
determinam a probabilidade por unidade de tempo para que um certo evento de difusão
possa ser termicamente ativado e foram usadas como parâmetros de entrada nas
simulações Monte Carlo Cinético. Nestas simulações, a geometria do sistema envolvido
é simplificada usando a aproximação Solid-on-Solid [31]. As partículas podem mover-
se de um sítio da rede para outro através de saltos discretos, que ocorrem com uma taxa
que, de acordo com a teoria do estado de transição [67] na aproximação harmônica [68],
tem a forma de Arrhenius, ou seja,
(
)
TkE
if
Bd
e
−
Γ=Γ
0
, (2.47)
onde
0
Γ é o prefator usualmente interpretado como uma freqüência tentativa de escape
de um átomo adsorvido no sítio meta-estável, T é a temperatura do sistema e E
d
é a
energia de ativação para a difusão do sítio i para o sítio f. As barreiras de energia E
d
são
39
obtidas a partir de cálculos DFT. Usamos, neste trabalho, um prefator comum igual a
10
13
s
-1
para todos os processos [31].
40
Capítulo 3
Reconstrução e relaxação da superfície (100)
de GaAs
O cristal de GaAs é caracterizado por ligações covalentes de caráter
parcialmente iônico entre os seus átomos. Logo, cada átomo tem seus orbitais
contraídos na forma de hibridizações sp
3
. Quando um cristal infinito é clivado de forma
a expor a sua superfície, os átomos da superfície passam a apresentar ligações pendentes
(“dangling bonds”) sp
3
parcialmente ocupadas. Este tipo de configuração eletrônica é
energeticamente desfavorável e a superfície sofre uma transição, através da mudança de
posições dos seus átomos, para um estado de energia mais baixa. Quando isto ocorre, ou
os átomos da superfície se relaxam, ou podem se reconstruir, dependendo da alteração
das suas novas coordenadas, que passam a minimizar a energia total do sistema. Na
relaxação, os átomos da superfície deslocam-se de suas posições originais porém,
conservam a simetria da superfície. Na reconstrução, por sua vez, a mudança nas
posições atômicas provoca não apenas uma relaxação, mas também uma mudança na
simetria da superfície [69].
Consideremos, então, uma superfície produzida ou por clivagem, ou por um
processo de crescimento. A disposição de seus átomos deverá corresponder a um
mínimo da energia livre da superfície. Para sistemas no estado de carga neutro, a
relaxação e a reconstrução podem ser explicadas por três princípios básicos: i) a energia
da superfície pode ser reduzida através da relaxação atômica, conduzindo o sistema à
sua natureza semicondutora. Neste caso, há uma transferência de carga entre os orbitais
pendentes do cátion e do ânion, aumentando a diferença de energia entre o orbital mais
alto ocupado (HOMO) e o mais baixo desocupado (LUMO). Esta troca conduz a uma
re-hibridização dos átomos da superfície, eliminando o possível caráter metálico desta
[70]; ii) a reconstrução da superfície pode reduzir a energia pela formação de novas
ligações, mudando o seu caráter para semicondutor, sendo a formação de dímeros uma
de suas características. Estes dímeros, por sua vez, tendem a se rearranjarem na
41
superfície. Para reduzir ainda mais a energia do sistema ocorre uma transferência de
carga de um dos átomos do dímero para o outro, fazendo com que o dímero sofra uma
inclinação com relação à superfície (“buckling”) , sendo que o átomo que recebe carga
ficará numa posição superior ao outro; iii) tanto a relaxação, como a reconstrução da
superfície, obedece à regra da contagem de elétrons. De acordo com esta regra, os níveis
de energia das ligações pendentes dos elementos eletropositivos (Ga, por exemplo)
estão na banda de condução e devem, portanto, estar vazios. Por outro lado, os níveis de
energia das ligações pendentes do elemento eletronegativo (As, por exemplo), estão na
banda de valência e devem estar preenchidos. Assim, a estrutura de energia mínima é
obtida através da transferência de elétrons das ligações pendentes dos elementos
eletropositivos para as ligações pendentes dos elementos eletronegativos.
Em particular, quando a regra de contagem de elétrons é aplicada a uma
superfície, como a (100) de GaAs, devemos levar em consideração algumas suposições
referentes à natureza de suas reconstruções [71]: a superfície se reconstrói na forma
(2xN), ii) na qual a periodicidade 2 surge da formação de dímeros na superfície e a
periodicidade N origina-se da ausência de dímeros na superfície, onde há apenas D
dímeros por célula unitária, onde D ≤ N. Para entendermos melhor como esta regra
funciona, consideremos um exemplo em que uma célula unitária da superfície apresenta
a reconstrução 2xN [71]. Supondo que a superfície seja terminada pelo elemento
eletronegativo (ânion), contaremos primeiro o número de elétrons exigidos para esta
reconstrução. Cada dímero de ânion da superfície precisa de 6 elétrons perfazendo um
total de 6D elétrons por célula unitária. Além disso, 8D elétrons são necessários para
ligar estes D dímeros aos cátions da segunda camada; dois elétrons para cada uma das
quatro ligações entre um dímero na camada superior e os cátions na segunda camada.
Logo, o total de elétrons necessários será 6D + 8D. Agora, contando o número de
elétrons disponíveis nas camadas reconstruídas, temos que o número de elétrons
disponíveis na camada superficial é 2V
n
D, onde V
n
é o número de elétrons de valência
dos elementos eletronegativos. O número de elétrons disponíveis na segunda camada é
2V
p
x N/2 = V
p
/N (uma vez que a metade dos elétrons totais de 2V
p
N da segunda camada
estão envolvidos com os átomos da terceira camada), na qual V
p
é o número de elétrons
de valência dos elementos eletropositivos. Para uma estrutura estável, o número de
elétrons necessários deve ser igual ao número de elétrons disponíveis ou seja:
6D + 8D = 2V
n
D + V
p
N.
42
Para GaAs, V
n
= 5, V
p
= 3, então, 4D = 3N. A menor célula unitária que satisfaz
a esta condição é a célula unitária 2x4, que tem três dímeros (D = 3), e o quarto está
ausente (N = 4). Contudo, cálculos realizados por Northrup e Froyen [72] mostram que
a reconstrução mais estável contém dois dímeros de As na camada superior, sendo
confirmados por observações realizadas através de STM de alta resolução [73]. Desta
forma, na aplicação da regra de contagem de elétrons, devemos supor a validade da
hibridização sp
3
para os estados de ligações pendentes. Caso contrário, esta regra não
funciona. Como exemplo, podemos citar a superfície (111) B de GaAs induzida por Sb,
pois as ligações pendentes re-hibridizadas não são ligações tipo sp
3
e a regra é violada
[74]. Um outro exemplo é a superfície (001) de GaSb, que é fracamente metálica e,
assim, a suposição da regra de contagem de elétrons é também violada [75].
3.1 Superfície de GaAs (100)
Superfícies semicondutoras, especialmente a superfície (100) de GaAs, são
tecnologicamente importantes devido às suas aplicações em dispositivos eletrônicos e
optoeletrônicos, além de serem convenientemente preparadas por MBE. A fim de
compreendermos as propriedades dessa superfície torna-se necessário o estudo de sua
estrutura.
A superfície (100) de GaAs é polar, ou seja, os planos paralelos à superfície
consistem de somente átomos de As, ou de somente átomos de Ga. Como conseqüência,
esta superfície exibe uma variedade de reconstruções, as quais dependem da sua
composição [76]. Experimentalmente, esta superfície é geralmente terminada com a
camada de As, devido ao uso excessivo deste material na câmara de MBE. Para
superfícies de GaAs (100) terminadas em As, a formação de dímeros de As reduz o
número de ligações pendentes por um fator igual a dois, tornando possível a
interpretação de imagens STM, através de padrões (2x4), conforme descrito no item
anterior, utilizando a regra de contagem de elétrons. Estudos teóricos baseados em
cálculos de energia total indicam que o padrão de reconstrução mais estável é a β2(2x4).
Embora o modelo mais estável para as superfícies (100) de GaAs, terminadas em As,
43
seja a reconstrução β2(2x4), o padrão (2x1) foi escolhido em nossas simulações pois,
durante as condições de crescimento do GaN sobre substratos de GaAs, que é a 700 °C,
este sistema sofre uma transição de fase estrutural para este modelo de reconstrução
[37]. Neste capítulo, mostraremos em detalhe, os resultados obtidos para as
propriedades de equilíbrio do GaAs e de sua superfície (100) na reconstrução (2x1).
Investigamos as conseqüências previstas nos resultados obtidos para a sua energia total
a partir de cálculos DFT, e discutimos a cinética da formação de ilhas durante a
deposição de uma monocamada de N, através da nitridação desta superfície de GaAs,
baseadas em simulações Monte Carlo Cinético.
3.1.1 Propriedades de Equilíbrio Cristalino do GaAs
Antes de proceder os cálculos para a superfície, calculamos a variação da
energia total do cristal de GaAs com relação a vários valores do parâmetro de rede para
obtermos uma boa descrição de suas propriedades de equilíbrio, além de checar a
precisão e a acurácia do método escolhido. Como discutido no capítulo 2, calculamos a
variação da energia total do bulk de GaAs em função da variação do volume da célula
unitária. Os resultados obtidos foram ajustados através da equação de Murnaghan e
obtivemos a
0
= 5.63 Å e B
0
= 0.75 Mbar como valores de equilíbrio para o parâmetro de
rede (a
0
) e para o Bulk Modulus (B
0
), respectivamente, quando os elétrons 3d do Ga
foram incluídos como estado de valência. Quando os elétrons 3d do Ga não foram
levados em consideração, os valores obtidos para o parâmetro de rede e para o Bulk
Modulus foram a
0
= 5.33 Å e B
0
= 0.75 Mbar, respectivamente. Estes resultados estão
em boa concordância com os dados experimentais (a
0
= 5.65 Å e B
0
= 0.77 Mbar) [77].
Os resultados foram convergidos para uma energia máxima da expansão de ondas
planas de 140 ou 80 Ry, dependendo da inclusão ou não dos elétrons 3d do Ga,
respectivamente, e foi usado um conjunto de 10 pontos k
r
especiais através da rede de
Monkhorst-Pack (4 4 2) usada para a amostragem da zona de Brillouin.
44
3.1.2 Aspectos Estruturais da superfície (100) de GaAs,
na reconstrução (2x1)
Antes de iniciarmos o nosso estudo da adsorção, analisamos a reconstrução
(2x1) da superfície (100) de GaAs, terminada em Ga e, em As, usando o modelo de
lâminas. As forças foram convergidas dentro do critério que elas sejam menores que 1
mHartree/Bohr. Na figura 3.1, mostramos, usando o modelo de lâminas, a geometria
atômica final obtida, após a reconstrução dos átomos da superfície (100) de GaAs e a
relaxação dos átomos da primeira e segunda camada abaixo desta. Como podemos
observar nesta figura, as superfícies se reconstroem criando dímeros. Na superfície
terminada em As, como descrito na figura 3.1 a, os dímeros de As não apresentam
inclinação com relação à superfície, ou seja, são simétricos. Na superfície terminada em
Ga, os dímeros apresentam uma certa inclinação, logo, estes são assimétricos, como
descrito na figura 3.1 b.
(a) (b)
Figura 3.1: Geometria atômica da superfície (100) de GaAs na reconstrução (2x1), vistas de lado, para
as superfícies: (a) terminada em As e (b) terminada em Ga.
Na figura 3.2, mostramos o resultado do teste de convergência da expansão em
ondas planas para a energia total destas superfícies, usando o método do
pseudopotencial. Notemos que a diferença na energia total, quando a energia de corte
varia de 55 para 60 Hartrees, é da ordem de 10
-5
Hartrees. Esta mesma diferença é
obtida quando a energia de corte varia de 60 para 65 Hartrees. Desta forma,
consideramos que, com a energia de corte na expansão da função de onda da ordem de
45
55 Hartrees a energia total que descreve o sistema converge com uma precisão de 0.027
meV/átomo para a superfície terminada em As. Para a superfície terminada em Ga,
analogamente, ela é de 0.023 meV/átomo.
20 30 40 50 60 70
1E-6
1E-5
1E-4
1E-3
0.01
0.1
1
10
Log (∆E) (Hartree)
Energia de Corte (Hartree)
20 30 40 50 60 70
1E-6
1E-5
1E-4
1E-3
0.01
0.1
1
10
Log (∆E) (Hartree)
Energia de Corte (Hartree)
(a) (b)
Figura 3.2: Testes de convergência para a superfície (100) de GaAs na reconstrução (2x1) para as
superfícies: (a) terminada em As e (b) terminada em Ga, usando o pseudopotencial de Troullier-Martins.
A energia total, em Hartree, é convergida com relação à energia cinética de corte.
Na tabela 3.1, mostramos os valores obtidos para os parâmetros estruturais
calculados para as superfícies reconstruídas que estão definidos na figura 3.3.
Analisando os dados da tabela 3.1, para a superfície terminada em Ga, os átomos de Ga
tentam formar orbitais híbridos sp
2
, pois o ângulo α tende a 120°. Por outro lado, para a
superfície terminada em As, os átomos de As procuram manter o caráter de orbital p
puro, visto que α não é muito diferente de 90
º
. Além disto, estes dímeros são simétricos,
pois ω ≈ 0°. O comprimento de ligação do dímero calculado é d
1
= 2.457 Å que está em
concordância com os valores experimentais, que se situa na faixa 2.2 < d
1
< 2.73 Å. É
interessante mencionar que, no modelo de reconstrução β2(2x4), os dímeros de As são
ligeiramente inclinados (ω = 0.5°) [78].
Para a superfície terminada em Ga, os dímeros são assimétricos, com um ângulo
de inclinação ω = 22.9° e o seu comprimento de ligação vale d
1
= 2.377 Å. Notamos
também, que há uma diferença entre os comprimentos de ligação Ga-As, indicados na
figura 3.3 por d
2
e d
3
. Para a superfície terminada em As, d
2
e d
3
assumem valores
iguais.
46
As propriedades eletrônicas da superfície (100) de GaAs, terminada em As, na
reconstrução (2x1) estão descritas no apêndice B, uma vez que estamos preocupados
neste trabalho, com a energia total e com os aspectos estruturais destas superfícies.
Figura 3.3: Representação esquemática dos parâmetros que descrevem a geometria atômica da
superfície (100) de GaAs, na reconstrução (2x1), obtidos através do modelo de lâminas.
Tabela 3.1: Parâmetros estruturais calculados, segundo a definição dada na figura 3.3, para as
superfícies (100) de GaAs comparados com resultados experimentais.
Superfície
α (°) ω(°)
d
1
(Å) d
2
(Å) d
3
(Å) d
1
(Exp)(Å) [79]
Terminada em Ga
127.7 22.9 2.38 2.39 2.47 -
Terminada em As
101.6 0.0 2.46 2.33 2.33 2.2 – 2.73
3.2 Adsorção de um átomo de N sobre a superfície
(100) de GaAs, terminada em As, na reconstrução
(2x1)
Estudamos a adsorção do átomo de N nos possíveis sítios da superfície
de GaAs (100) na reconstrução (2x1) que estão descritos na figura 3.4. Na figura 3.5,
mostramos os nossos resultados para a variação da energia de adsorção do N nestes
sítios, em função da altura ∆z deste sobre a superfície. Como referência marcamos a
posição do dímero de As na superfície como ∆z = 0, uma vez que eles são simétricos, e
fizemos ∆E
ads
= 0 para a energia de adsorção do N quando este está a
47
aproximadamente 4 Å do dímero, para o qual não há variação sensível da E
Total
do
sistema em função do aumento da altura ∆z do N sobre esta superfície.
Figura 3.4: Vista do topo dos possíveis sítios de adsorção de N sobre superfícies (100) de GaAs,
terminadas em As, na reconstrução (2x1).
-5 -4 -3 -2 -1 0 1 2 3 4
-3
-2
-1
0
1
2
Sítios
A
B
C
D
E
F
Energia de Adsorção (eV)
∆Z (Å)
Figura 3.5: Variação da energia de adsorção em função da altura ∆z (em Å) do átomo de N adsorvido
sobre a superfície (100) de GaAs, terminada em As, na reconstrução (2x1).
Na figura 3.5, podemos verificar que o átomo de N prefere ser adsorvido
nos sítios C, D e E, respectivamente, que são os sítios desocupados da superfície. É
interessante notar que a adsorção sobre os dímeros de As (no sítio A) é 0.75 eV menos
estável que sobre o sítio C, que tem uma energia de adsorção de 3.35 eV. Também, a
adsorção nos sítios D e E são 0.54 e 0.58 eV menos estáveis que no sítio C,
respectivamente. Notamos que, diferente do que foi observado para a adsorção de N
sobre o sítio de menor energia da superfície (110) de GaAs [80], não existem barreiras
48
de ativação antes do início do processo de adsorção do N. Também verificamos que as
energias envolvidas no processo de nitridação são relativamente maiores quando
comparadas com as energias da homoepitaxia de GaAs, que tem aproximadamente 2 eV
em módulo [31].
Na tabela 3.2, mostramos os valores calculados para os parâmetros estruturais da
lâmina (como definidos na figura 3.3) para a adsorção de N nos sítios A, B, C, D, E e F,
após considerarmos a relaxação do N e das duas camadas de vizinhos mais próximos,
abaixo da superfície.
Tabela 3.2: Parâmetros estruturais calculados, como definidos na figura 3.3, para a adsorção de N
sobre superfícies (100) de GaAs, na reconstrução (2x1), nos sítios A, B, C, D, E e F.
Sítios
α (°) ω(°)
d
1
(Å) d
2
(Å) d
3
(Å) d
N-Ga
(Å) d
N-As
(Å)
A
93.6 0.0 3.48 2.34 2.34
- 1.74
B
91.5 0.0 2.52 2.33 2.33
1.66 1.73
C
102.9 11.8 2.44 2.88 2.31
1.64 2.03
D
103.6 0.0 2.39 2.34 2.34
2.98 1.90
E
112.6 0.0 2.38 2.30 2.30
1.70 1.83
F
105.9 0.0 2.63 2.37 2.33
- 1.70
A partir dos resultados obtidos, notamos que, somente para a adsorção no sítio
E, ocorreu uma ligeira distorção no ângulo α em direção a uma configuração planar das
ligações de As. Nos sítios A e B, ao contrário, houve uma variação do ângulo α na
direção de uma configuração p pura das ligações do As, isto é, 90
º
. Somente no sítio C,
detectamos uma inclinação do dímero de As, causada pela presença do N, que afasta um
dos As de seu átomo de Ga primeiro vizinho, devido ao aumento do valor do parâmetro
d
2
. Exceto para a adsorção nos sítios A, D e F, o N prefere permanecer próximo aos
átomos de Ga. Porém, a distância de ligação entre o Ga e o N na adsorção nos sítios B,
C e E é, aproximadamente, 13,6% menor que a mesma distância quando consideramos o
GaN bulk (a
0
~ 4.52 Å). Finalmente notamos que, exceto para o sítio C, há pouca
variação na distância de ligação entre o As e o Ga no substrato (parâmetros d
2
e d
3
). A
fim de compreendermos as variações dos comprimentos de ligação e ângulos de ligação,
mostramos na figura 3.6, as posições atômicas relaxadas do N e de seus vizinhos após a
sua adsorção nos sítios A, C, D e E.
49
(a) (b)
(c) (d)
Figura 3.6: Posições atômicas relaxadas após a adsorção de N nos sítios: (a) A, (b) C, (c) D e (d) E da
superfície (100) de GaAs, na reconstrução (2x1).
É interessante observar que o N, na figura 3.6 b, a fim de permanecer próximo
aos átomos de Ga, provoca uma torção nos dímeros de As, empurrando-os para fora da
superfície.
Logo, baseado em nossos resultados, podemos concluir que os átomos de N
preferem os sítios desocupados da superfície, próximos ao dímero de As, para provocar
uma possível troca iônica entre os átomos de N e As.
As propriedades eletrônicas da superfície (100) de GaAs, terminada em As, na
reconstrução (2x1) com um átomo de N adsorvido estão descritas no apêndice C, uma
vez que estamos preocupados neste trabalho, com a energia total e com os aspectos
estruturais destas superfícies.
50
3.3 Adsorção do segundo N sobre a superfície (100) de
GaAs, terminada em As, na reconstrução (2x1)
Antes de iniciarmos as simulações Monte Carlo da nitridação das superfícies
(100) (2x1) de GaAs terminadas em As, simulamos a adsorção de um segundo átomo de
N, mantendo o primeiro N no sítio C, a fim de verificarmos como se alteram as energias
de adsorção do sistema, agora mais rico em N.
Mostramos, na figura 3.7, os nossos resultados para a variação da energia de
adsorção do segundo N nos sítios de adsorção definidos na figura 3.4, em função da
altura ∆z deste sobre a superfície (100) (2x1) de GaAs. Para este caso, também
investigamos a energia de adsorção no sítio C’, na posição oposta ao sítio C, ou seja, no
outro lado do dímero a fim de verificarmos se a troca iônica entre os átomos de N e As é
possível.
Para a adsorção do segundo átomo de N sobre a superfície (100) de GaAs,
verificamos que o sítio mais estável para a adsorção deste átomo é, neste caso, tanto o
sítio D, quanto o sítio C’, cujas energias de adsorção são 3.12 e 3.11 eV,
respectivamente. Comparando com o caso anterior, em que, para a adsorção de um
átomo de N, o sistema iria precisar de uma energia de 3.35 eV para executar este
processo, as energias necessárias para a adsorção do segundo N são, em média, 0.17 eV
menores, conforme descrito na tabela 3.3, em que mostramos um resumo dos valores
obtidos para a energia de adsorção, E
ads
, de um átomo de N e de dois átomos de N
adsorvidos respectivamente, bem como a altura ∆z destes átomos sobre esta superfície.
Uma vez que este valor está próximo de 3.12 e 3.11 eV, escolhemos 3.35 eV
como uma boa aproximação para a energia de adsorção dos diversos átomos de N que
utilizaremos em nossas simulações Monte Carlo Cinético para a nitridação.
51
Figura 3.7: Variação da energia de adsorção em função da altura ∆z (em Å) do segundo átomo de N
adsorvido sobre a superfície (100) de GaAs, terminada em As, na reconstrução (2x1), nos sítios definidos
na figura 3.3.
Na figura 3.7, também notamos que a adsorção do segundo N no sítio A (sobre
os dímeros de As) é 0.65 eV menos estável que o sítio D. Os sítios de adsorção E e C
são 0.32 e 0.53 eV, menos estáveis que o sítio D, respectivamente.
Na tabela 3.3, podemos observar que para os sítios D, E e F, a altura ∆z dos
átomos de N assumem valores aproximadamente iguais. Quando a posição do átomo
adsorvido e as posições dos átomos das duas primeiras camadas da superfície são
otimizadas, a seqüência de estabilidade dos sítios sofre uma certa alteração, mas o sítio
C ainda continua sendo o sítio mais estável. Desta forma, sistema gasta mais energia
para relaxar nos sítios A, D e F. Nos sítios C, E e B esta energia é menor, logo, a
hipótese de que a deposição dos átomos deva começar em cima do sítio C torna-se
plausível. Na figura 3.8, mostramos as posições relaxadas do segundo átomo de N e de
seus vizinhos após a sua adsorção nos sítios mais estáveis, ou seja, sítio D e sítio C’. Na
figura 3.8 b, podemos verificar que quando um segundo átomo de N é adsorvido no sítio
C’, este empurra os átomos de As para fora da superfície, podendo ocorrer, desta forma,
a troca iônica entre os átomos de N e As, com uma conseqüente dessorção de As
2
.
52
(a) (b)
Figura 3.8: Posições atômicas relaxadas após a adsorção do segundo átomo de N: (a) no sítio D e (b) no
sítio C’, da superfície (100) de GaAs, na reconstrução (2x1).
Tabela 3.3: Valores obtidos para a energia de adsorção, E
ads
, para um átomo de N e para dois átomos de
N adsorvidos, como definidos na figura 3.3, sobre superfícies (100) de GaAs, na reconstrução (2x1), nos
sítios A, B, C, D, E e F, bem como a energia de adsorção do sistema relaxado, E
rel
para a adsorção de
um átomo de N.
GaAs(100) : N GaAs(100) : N
2
Sítios
E
ads
(eV) E
Rel
(eV)
∆z (Å)
E
ads
(eV)
∆z (Å)
A
-2.60 -3.17 1.33 -2.47 -1.87
B
-2.47 -3.36 1.07 -2.30 0.27
C
-3.35 -3.72 0.36 -2.59 1.33
C’
- - - -3.11 0.03
D
-2.81 -2.83 -2.40 -3.12 -2.4
E
-2.77 -3.71 -0.80 -2.80 -0.80
F
-2.22 -2.32 1.60 -2.23 1.73
53
3.4 Energia de Formação da adsorção de N na
superfície (100) de GaAs, terminada em As
A fim de verificarmos se a formação de dímeros da N na adsorção é possível e,
além disto, estável, calculamos a energia de formação em duas situações: i) superfície
de GaAs terminada com um dímero N-As, ii) terminada com um dímero N-N. Para tal,
calculamos as energias de formação, γ, destes sistemas, que foram definidas como[81]:
µµγ ∆−−−= )(
'
NGaeet
nnnEA ,
onde
'
t
E é a energia de coesão da supercélula contendo n
Ga
partículas de Ga, n
N
partículas de N e n
e
é o número de elétrons de um reservatório com um potencial
químico µ
e
. As restrições nos potenciais químicos são:
0 < µ
e
≤ E
gap
e -∆H
f
≤ ∆µ ≤ ∆H
f,
sendo que E
gap
é a energia do “gap” e ∆H
f
é o calor de formação do GaN e vale 1.14
eV. Uma vez que estamos lidando com estados neutros de carga, podemos considerar
que n
e
= 0. A área da superfície (100) de GaAs, terminada em As, na reconstrução (2x1)
é dada por 28.41 Å
2
.
Na figura 3.9 mostramos a variação da energia de formação, γ, em função do
potencial químico para as superfícies terminadas com o dímero N-As e com o dímero
de N-N, comparados com o dímero de As
2
, todos na reconstrução (2x1). Podemos
verificar nesta figura que, para ambientes ricos em N, o dímero N-N é mais estável que
os dímeros N-As e As
2
, respectivamente. Para ambientes pobres em N, o dímero As
2
é
mais estável que os dímeros N-N e N-As, respectivamente. Além disto, a situação em
que temos o dímero N-As não é possível de ser obtida (o sistema passa de um dímero de
As
2
para N
2
diretamente).
54
-1.5 -1.0 -0.5 0.0 0.5 1.0 1.5
-6
-5
-4
-3
-2
-1
dímero N-As
dímero N-N
dímero As-As
γ
(eV.átomos
-1
.Å
-2
)
∆µ(eV)
Figura 3.9: Variação da energia de formação γ em função do potencial químico ∆µ para a superfície
(100) de GaAs, na reconstrução (2x1), terminada em As, para os dímeros N-As e N-N.
Na tabela 3.3, mostramos os valores calculados para os parâmetros estruturais,
como definidos na figura 3.3, para os dímeros N-N e N-As. Da análise desta tabela
podemos observar que para o dímero N-N, os átomos de N tentam formar orbitais
híbridos sp
2
, pois o ângulo α tende a 120°. Por outro lado, para o dímero N-As, os
átomos de As e N procuram manter o caráter de orbital p puro, pois o ângulo α tende a
90°. O comprimento de ligação d
1
do dímero é 1.94 Å. Para o dímero N-As, os dímeros
são assimétricos, com um ângulo de inclinação ω = 23.1° e comprimento de ligação
d
1
= 1.25 Å. Para ilustrarmos as variações dos comprimentos de ligação e dos ângulos
de ligação, mostramos na figura 3.10, posições atômicas relaxadas dos dímeros N-As e
N-N.
Tabela 3.3: Valores obtidos para os parâmetros estruturais para os dímeros N-As e N-N.
Dímeros
α (°) ω(°)
d
1
(Å) d
2
(Å) d
3
(Å)
N-As
106.3 23.1 1.94 1.98 2.44
N-N
118.4 0.0 1.25 2.25 2.25
55
(a) (b)
Figura 3.10: Posições atômicas relaxadas dos dímeros N-As e N-N.
3.5 Difusão do N
Durante a nitridação, outros processos, além da adsorção, podem ocorrer. Dentre
eles podemos citar a difusão do átomo de N sobre a superfície de GaAs cujo estudo será
descrito neste item.
Na figura 3.11, mostramos a variação da energia total como uma função do
deslocamento do N, iniciando a partir do sítio de adsorção E, ao longo da direção [-1 1
0], simulando a difusão de N ao longo desta direção.
Figura 3.11: Variação da energia total como uma função do deslocamento do N, iniciando a partir do
sítio de adsorção E, ao longo da direção [-1 1 0], simulando a difusão de N sobre a superfície (100) de
GaAs, na reconstrução (2x1). As letras denotam os sítios de adsorção na superfície que o N irá passar.
56
Na figura 3.11, notamos que, para uma difusão completa ao longo da direção [-1
1 0], o N deve sobrepor uma barreira de energia de 7.34 eV. Se a difusão ocorre ao
longo da direção [ 1 1 0], ou sobre o dímero de As, esta barreira de energia é 5.50 eV,
ou aproximadamente 60.0 eV (a barreira é muito alta), respectivamente, como podemos
observar na figura 3.12 (a) e (b), que descreve a variação da energia total como uma
função do deslocamento de um átomo de N, iniciando a partir do sítio de adsorção E, ao
longo da direção [1 1 0] e sobre o dímero de As, respectivamente. Assim, podemos
inferir que os átomos de N têm baixa mobilidade quando comparados com a difusão de
Ga ou As
2
, durante o processo de nitridação, mesmo na temperatura deste processo, T ~
526°C (800 K). Devido a esta baixa mobilidade, esperamos que, após o processo da
nitridação, a superfície resultante seja rugosa.
Desta forma, dois cenários de crescimento podem, então, ocorrer: i) o N pode ser
adsorvido sobre o dímero da superfície e, neste caso, nenhuma difusão ocorre; e ii) o N
pode ser adsorvido sobre os sítios vazios da superfície. Neste caso, o N pode difundir
lentamente sobre a superfície. Em ambos os casos há a possibilidade da obtenção de
superfícies rugosas.
Figura 3.12: Variação da energia total como uma função do deslocamento do N, iniciando a partir do
sítio de adsorção E: a) ao longo da direção [1 1 0] e b) sobre o dímero de As, simulando a difusão de N
sobre a superfície (100) de GaAs, na reconstrução (2x1). As letras denotam os sítios de adsorção na
superfície que o N irá passar
57
3.6 Simulações Monte Carlo Cinético da Nitridação
Tendo como objetivo verificar como a baixa mobilidade do N e sua alta energia
de adsorção afetam o processo de nitridação, realizamos simulações Monte Carlo
Cinético da deposição do N sobre a superfície de GaAs (100), terminada em As, na
reconstrução (2x1), baseadas na aproximação “Solid-on-solid” [31], conforme descrito
no capítulo 2. Para a simulação do processo de nitridação, levamos em consideração
somente os processos de adsorção de N e de difusão de N, com um prefator comum de
10
13
s
-1
[31]. Os nossos resultados foram analisados quantitativamente pelo estudo da
variação da densidade de ilhas tomadas em função da variação do fluxo de N, nas
temperaturas de crescimento iguais a 600, 650, 700, 750 e 800 K, mantendo fixa a
cobertura em 0.5 ML (monocamada), cujo resultado está descrito na figura 3.13. Estas
temperaturas de crescimento equivalem à 327, 377, 427, 477 e 527 °C, sendo que a
partir da última temperatura, 527 °C o procedimento experimental é realizado.
0.05 0.10 0.15 0.20 0.25 0.30 0.35
5.800
5.805
5.810
5.815
5.820
800 K
750 K
700 K
650 K
600 K
log (densidade de ilhas)
Fluxo
N
(ML/s)
Figura 3.13: Variação do logaritmo da densidade de ilhas em função da variação do fluxo de N, obtidos
da simulação Monte Carlo Cinético do processo de nitridação obtidos.
Analisando a figura 3.13, notamos que, para valores de fluxo do N mais baixos,
a densidade de ilhas é muito alta. Observamos também que independente da
temperatura usada durante a simulação da deposição do N, a densidade de ilhas diminui
58
quando o fluxo de N é alto, ou seja, a partir de 0.30 ML/s, favorecendo o aumento da
rugosidade da superfície nitridada.
Nas figuras 3.14 e 3.15 mostramos a variação da densidade de ilhas formadas em
função da cobertura de N, nas temperaturas de crescimento iguais a 300 e 700 K,
respectivamente, usando taxas de deposição 0.15, 0.25 e 0.35 ML/s.
0.1 0.2 0.3 0.4 0.5
4.6
4.8
5.0
5.2
5.4
5.6
5.8
6.0
T = 300 K
Fluxo = 0.15 ML/s
Fluxo = 0.25 ML/s
Fluxo = 0.35 ML/s
log (densidade de ilhas)
Θ
N
(ML)
Figura 3.14: Variação calculada do logaritmo da densidade de ilhas como uma função da cobertura de
N (a temperatura de crescimento é 300 K).
Tanto para uma temperatura de crescimento de 300 K quanto para uma
temperatura de 700 K, como mostrado nas figuras 3.14 e 3.15, um importante aspecto a
ser observado nestas figuras é a saturação da densidade de ilhas após a deposição de
uma certa quantidade de material. Como podemos observar nestas figuras, a densidade
de ilhas se satura para uma cobertura de N da ordem de 0.3 ML.
59
0.1 0.2 0.3 0.4 0.5
4.6
4.8
5.0
5.2
5.4
5.6
5.8
6.0
T = 700 K
Fluxo = 0.15 ML/s
Fluxo = 0.25 ML/s
Fluxo = 0.35 ML/s
log (densidade de ilhas)
Θ
N
(ML)
Figura 3.15: Variação calculada do logaritmo da densidade de ilhas como uma função da cobertura de
N (a temperatura de crescimento é 700 K).
Na figura 3.16, mostramos os resultados obtidos para a variação da densidade de
ilhas em função da variação da temperatura de crescimento. Nestas simulações a
quantidade de material depositado foi igual a 0.5 ML e as taxas de deposição variaram
de 0.05, 0.15, 0.25 e 0.35 ML/s
600 650 700 750 800
5.802
5.804
5.806
5.808
5.810
5.812
5.814
5.816
5.818
5.820
θ
N
= 0.5 ML
Fluxo = 0.05 ML/s
Fluxo = 0.15 ML/s
Fluxo = 0.25 ML/s
Fluxo = 0.35 ML/s
Log(densidade de ilhas)
Temperatura (K)
Figura 3.16: Variação calculada do logaritmo da densidade de ilhas como uma função da temperatura
de crescimento. Todas as simulações foram realizadas com uma cobertura igual a 0.5 ML.
60
Na figura 3.16, notamos que para valores de fluxos de N altos, menor é a
densidade de ilhas formadas. Observamos também um decréscimo na densidade de ilhas
quando aumentamos a temperatura de 750 para 800 K. Para uma temperatura de
aproximadamente 650 K e um fluxo de N entre 0.05 e 0.15 ML/s, as condições são
favoráveis para uma distribuição uniforme de ilhas sobre o substrato
Nas figuras 3.17 (a) e (b), mostramos a variação do logaritmo da densidade de
ilhas em função do fluxo de N, nas temperaturas de crescimento 600 e 800 K,
respectivamente. Neste caso consideramos, em nossas simulações o “annealing” da
amostra após a nitridação à temperatura de 500 K, ou 227 °C. Nestas simulações a
quantidade de material depositado foi igual a 0.5 ML.
(a)
(b)
Figura 3.17: Variação calculada do logaritmo da densidade de ilhas como uma função do fluxo de N. Em
(a) a temperatura de crescimento usada foi de 600 K e em (b) a temperatura de crescimento usada foi de
800 K.
61
Na figura 3.17, podemos verificar que, para um fluxo de N baixo (da ordem de
0.05 ML/s), a densidade de ilhas é mais alta nas amostras que sofreram “annealing”
independente da sua temperatura de crescimento.
Nas figuras 3.18 e 3.19, mostramos os resultados obtidos de uma simulação
Monte Carlo Cinético da nitridação da superfície (100) (2x1) de GaAs, escolhida ao
acaso. A cobertura de N utilizada foi 0.5 ML, o fluxo de N foi igual a 0.05 ML/s e
temperatura de crescimento igual a 800 K, não levando em consideração o “annealing”
e considerando o “annealing” da amostra após a nitridação à temperatura de 500 K, ou
227°C, respectivamente. O tamanho da rede é (100 x 100) Å.
(a)
(b)
Figura 3.18: Resultado, escolhido ao acaso, de uma simulação Monte Carlo Cinético da nitridação da
superfície (100) (2x1) de GaAs. A cobertura usada foi igual a 0.5 ML e a temperatura de crescimento
igual a 800 K. Neste caso, o “annealing” não foi levado em consideração. O tamanho da rede é (100 x
100) Å. (a) Vista do topo e; (b) Vista em perspectiva.
62
(a)
(b)
Figura 3.19: Resultado, escolhido ao acaso, de uma simulação Monte Carlo Cinético da nitridação da
superfície (100) (2x1) de GaAs. A cobertura usada foi igual a 0.5 ML, o fluxo de N igual a 0.05 ML/s e a
temperatura de crescimento igual a 800 K. Neste caso, consideramos o “annealing” da amostra após a
nitridação à temperatura de 500K, ou 227°C. O tamanho da rede é (100 x 100) Å. (a) Vista do topo e; (b)
Vista em perspectiva.
Analisando as figuras 3.18 e 3.19, podemos notar o efeito da rugosidade da
superfície (100) (2x1) de GaAs nitridada. Os átomos de N caem sobre outros átomos de
N, podendo, assim, formar, primeiramente, várias ilhas, para depois ocorrer a nucleação
e a formação das camadas depositadas de N.
Com base nos nossos resultados, mostramos que a densidade de ilhas é muito
alta para valores de fluxo de N mais baixos, favorecendo, desta forma, uma distribuição
de N mais uniforme sobre o substrato. Para valores de fluxo de N mais altos, ocorre um
decréscimo na densidade de ilhas, independente da temperatura usada durante a
simulação da deposição do N, o que favorece o aumento da rugosidade da superfície
nitridada. Verificamos, também, que a densidade de ilhas é mais alta quando o
63
“annealing” é usado no crescimento, para temperaturas de crescimento mais altas,
principalmente, quando o fluxo de N é baixo.
Finalmente, os nossos resultados mostram que, na nitridação da superfície (100)
(2x1) de GaAs terminada em As, não há indícios da formação de “wetting layers” que
caracterizaria o modo de crescimento tipo camada-a-camada ideal (modo Frank-van-
der-Merwe) e nem do modo de crescimento de Stranski-Krastanov, em que após se
terem formado algumas camadas completas, se inicia o crescimento de aglomerados
tridimensionais de partículas. Ao contrário, eles apontam que o modo Volmer-Weber
para o crescimento deve ocorrer, na qual primeiro deve haver a formação de várias
ilhas, para depois haver a nucleação e a formação das camadas depositadas de N [4]. Na
figura 3.20, mostramos estes três modos de crescimento.
Figura 3.20: Modos de crescimento de superfícies.
64
Capítulo 4
Conclusões
Apresentamos neste trabalho um estudo sistemático da deposição de N sobre
substratos de GaAs (100) usando a Teoria do Funcional da Densidade na aproximação
LDA, expansão da função de onda em ondas planas e o método do Pseudopotencial, em
conjunto com simulações de Monte Carlo Cinético em um modelo de gás de rede.
Iniciamos o nosso estudo calculando os parâmetros estruturais para o GaAs bulk
e paras a superfícies (100) de GaAs na reconstrução (2x1). Os valores calculados para o
seu parâmetro de rede e o seu “Bulk Modulus” estão em boa concordância com os
valores experimentais. Nossos resultados para as propriedades estruturais das
superfícies (100) de GaAs, na reconstrução (2x1), também estão em concordância com
os dados experimentais disponíveis.
Para a adsorção de um átomo de N na superfície (100) de GaAs, nossos
resultados mostram que este átomo prefere ser adsorvido sobre os sítios desocupados da
superfície. Na adsorção do segundo átomo de N sobre a superfície (100) de GaAs,
verificamos que estas tendências são mantidas, porém as energias de adsorção são em
média 0.23 eV menores. Além disto, a situação de menor energia é aquela em que os
dois átomos de N podem fazer uma troca iônica com os átomos de As da superfície.
Quanto à energia de formação, para ambientes ricos em N, o dímero N-N é mais
estável que os dímeros N-As e As
2
. Para ambientes pobres em N, o dímero As
2
é mais
estável que os dímeros N-N e N-As. Os resultados corroboram a observação de que há
uma troca iônica por pares entre o N
2
e o As
2
verificada quando estudamos os aspectos
estruturais deste sistema.
Os nossos resultados mostram que dois cenários de crescimento podem, então,
ocorrer para a deposição do átomo de N: i) o N pode ser adsorvido sobre o dímero da
superfície, neste caso, nenhuma difusão ocorre; e ii) o N pode ser adsorvido sobre os
65
sítios vazios da superfície. Neste caso, o N pode difundir lentamente sobre a superfície,
ocorrendo assim a possibilidade da obtenção de superfícies rugosas.
Para a simulação do processo de nitridação, os nossos resultados mostram que
independente da temperatura usada durante a simulação da deposição do N, a densidade
de ilhas diminui quando o fluxo de N é alto. Verificamos, também, um decréscimo na
densidade de ilhas quando a temperatura é aumentada de 750 para 800 K e para um
fluxo de N baixo, da ordem de 0.05 ML/s de acordo com nossos cálculos, a densidade
de ilhas é mais alta quando o “annealing” é levado em consideração, tanto para a
temperatura de crescimento igual a 600 K quanto para 800 K.
Finalmente, os nossos resultados mostram que, na nitridação da superfície (100)
(2x1) de GaAs terminada em As, não há indícios da formação de “wetting layers” que
caracterizaria o modo de crescimento tipo camada-a-camada ideal (modo Frank-van-
der-Merwe). Ao contrário, eles apontam que o modo Volmer-Weber para o crescimento
deve ocorrer, na qual primeiro deve haver a formação de várias ilhas, para depois haver
a nucleação e a formação das camadas depositadas de N [4].
Uma continuidade do trabalho aqui apresentado é o estudo da rugosidade, ou
seja, da função correlação das diferenças de alturas da deposição de átomos sobre a
superfície, bem como verificar a ocorrência da troca iônica entre os átomos de N e As,
quando um segundo átomo de N chega à superfície.
66
Apêndice A
Pseudopotenciais para o H, As, N e Ga
O primeiro passo para o cálculo das propriedades estruturais e eletrônicas de um
determinado sistema é a escolha do pseudopotencial que descreve de forma eficaz as
propriedades dos átomos que compõe o sistema a ser estudado.
Neste apêndice, descreveremos os procedimentos usados na determinação dos
pseudopotenciais para os átomos de H, As, N e Ga, os quais foram utilizados em nossos
cálculos. Para tal, construímos pseudopotenciais, segundo o método proposto por
Troullier-Martins [47], através do programa computacional fhi98PP [36].
Segundo o método de Troullier-Martins, devemos ajustar o valor para o raio de
corte (r
c
) de cada orbital de valência do elemento em questão. Quanto maior r
c
,
menor
será a expansão em ondas planas para o cálculo de energia no método de ondas planas.
De modo geral, r
c
é escolhido de forma a se situar entre o último nó e o último pico da
função de onda “all electron” do orbital em questão, obtido de um cálculo atômico
relativístico. Para o átomo de H, que possui uma configuração atômica do tipo 1s
1
,
utilizamos raio de corte igual a 1.276 raios de bohr para todas as componentes do
momento angular (l = 0, 1 e 2). As funções de onda “all electron” e as pseudofunções
de onda para as componentes do momento angular l = 0, 1 e 2 para o átomo de H estão
mostradas na figura A.1.
Para o átomo de As, que possui uma configuração atômica do tipo [Ar] 3d
10
4s
2
4p
3
, usamos raios de corte iguais a 2.064 raios de bohr para l = 0, 2.167 raios de bohr
para l = 1 e 2.448 raios de bohr para l = 2. A comparação entre as funções de onda “all
electron” e as pseudofunções de onda são mostradas na figura A.2.
Analogamente, para o átomo de N, cuja configuração atômica é do tipo [He] 2s
2
2p
3
, o raio de corte usado foi 1.416 raios de bohr para todas as componentes do
momento angular (l = 0, 1 e 2). Na figura A.3, mostramos a comparação entre as
67
funções de onda “all-electron” e as pseudofuncões de onda para as componentes do
momento angular l = 0, 1 e 2 para o átomo de N.
Figura A.1: Função de onda radial “all electron”(linhas pontilhadas) e pseudofunções de onda obtidas
para as componentes l = 0, ou orbital 1s (linha preta), l = 1, ou orbital 2p (linha azul) e l = 2, ou
orbital 3d (linha verde), para o átomo de H.
Figura A.2: Função de onda radial “all electron”(linhas pontilhadas) e pseudofunções de onda obtidas
para as componentes l = 0, ou orbital 1s (linha preta), l= 1, ou orbital 2p (linha azul) e l = 2, ou orbital
3d (linha verde), para o átomo de As.
68
Figura A.3: Função de onda radial “all electron”(linhas pontilhadas) e pseudofunções de onda obtidas
para as componentes l = 0, ou orbital 1s (linha preta), l = 1, ou orbital 2p (linha azul) e l = 2, ou
orbital 3d (linha verde), para o átomo de N.
Finalmente, para o átomo de Ga, com as configurações atômicas do tipo [Ar
3d
10
] 4s
2
4p
1
e [Ar] 3d
10
4s
2
4p
1
, os raios de corte utilizados foram 2.042 raios de bohr
para l = 0 e 2.251 raios de bohr para l = 1 . Neste cálculo não levamos em
consideração o orbital 3d do Ga. A comparação entre as funções de onda “all-electron”
e as pseudofunções de onda podem ser vistas na figura A.4. Quando o orbital 3d do Ga
é levado em consideração, os raios de corte utilizados foram 2.042 raios de bohr para
l = 0, 2.251 raios de bohr para l = 1 e 2.042 raios de bohr para l = 2. Na figura A.5,
mostramos a comparação entre as funções de onda “all-electron” e as pseudofunções de
onda para as componentes do momento angular l = 0, 1 e 2.
A fim de verificarmos a eficiência e a transferabilidade do pseudopotencial
escolhido torna-se necessário a realização de alguns testes. O primeiro teste a ser feito
em um pseudopotencial é a comparação entre a derivada logarítmica da função de onda
“all electron” com a derivada logarítmica da pseudofunção de onda, as quais devem ser
semelhantes na região dos estados de referência, ou seja, para valores de energia em
torno dos autovalores de cada orbital. Mostramos a comparação das derivadas
logarítmicas nas figuras A.6, A.7, A.8, A.9 e A.10, as derivadas logarítmicas para os
átomos de H, As, N e Ga , respectivamente.
69
Figura A.4: Função de onda radial “all electron” (linhas pontilhadas) e pseudofunções de onda obtidas
para as componentes l = 0, ou orbital 1s (linha preta), l = 1, ou orbital 2p (linha azul) para o átomo de
Ga.
Figura A.5: Função de onda radial “all electron”(linhas pontilhadas) e pseudofunções de onda obtidas
para as componentes l = 0, ou orbital 1s (linha preta), l = 1, ou orbital 2p (linha azul) e l = 2, ou orbital
3d (linha verde), para o átomo de Ga.
70
Figura A.6: Derivada logarítmica para orbitais s, p e d do átomo de H.
Figura A.7: Derivada logarítmica para orbitais s, p e d do átomo de As.
71
Figura A.8: Derivada logarítmica para orbitais s, p e d do átomo de N.
Figura A.9: Derivada logarítmica para orbitais s e p do átomo de Ga.
72
Figura A.10: Derivada logarítmica para orbitais s, p e d do átomo de Ga.
Uma vez que o pseudopotencial foi gerado e testado é preciso que testes com
sistemas reais (como moléculas ou sólidos) sejam realizados e que os valores obtidos
para as propriedades estruturais e eletrônicas sejam comparados com os valores obtidos
experimentalmente. Normalmente, aceitam-se erros da ordem de 5 % para o parâmetro
de rede e de 10 % para o “Bulk Modulus”. Os nossos resultados para estes testes foram
apresentados na seção 4.1.1.
73
Apêndice B
Propriedades eletrônicas da superfície (100)
de GaAs, terminada em As, na reconstrução
(2x1)
Mostramos na figura B.1, a estrutura eletrônica da superfície (100) de GaAs,
terminada em As, na reconstrução (2x1). A área hachurada representa a projeção da
estrutura eletrônica do cristal perfeito nesta superfície e a zona de Brillouin para esta
superfície, está descrita na figura B.2.
Analisando a figura B.1, observamos que a superfície (100) de GaAs, terminada
em As, na reconstrução (2x1), de acordo com os nossos resultados, apresenta onze
estados de superfície. O primeiro e o segundo estados desocupados, C
1
e C
2
, estão a
0.70 e 1.25 eV, respectivamente, abaixo da banda de condução, ao redor do ponto
K
. O
estado C
3
está cerca de 1.31 eV abaixo da banda de condução ao longo da direção
'JK − . O estado C
4
, por sua vez, situa-se a 2.57 eV abaixo da banda de condução ao
longo da direção 'JK − , cruzando nessa direção com o estado V
1
, que está a 1.87 eV
acima da banda de valência. Deste modo, podemos, baseados em nossos resultados,
afirmar que esta superfície apresenta um certo caráter metálico. O estado V
2
está bem
localizado ao redor do ponto
K
e a 0.50 eV acima da banda de valência. O estado V
3
da
superfície corresponde à interação entre os átomos da segunda e terceira camadas da
superfície. No “gap” interno, observamos quatro estados de superfície, V
4,
V
5,
V
6
e V
7
com combinações do estado s localizadas nos átomos de As da superfície.
Na figura B.3, mostramos a densidade de carga para a superfície (100) de GaAs,
terminada em As, na reconstrução (2x1), vistas de lado e do topo. Nesta figura,
podemos observar que a densidade de carga está concentrada entre os átomos que
compõem o dímero de As e sobre os átomos de H, que saturam a superfície,
representada pelas linhas verdes.
74
Figura B.1: Estrutura eletrônica para a superfície (100) de GaAs, terminada em As, na reconstrução
2x1. A área hachurada representa a projeção da estrutura eletrônica do cristal. As letras denotam os
estados de superfície.
'J
J
K
Γ
[110]
[110]
Figura B.2: Zona de Brillouin para a superfície (100) de GaAs, terminada em As, na reconstrução 2x1.
(a) (b) (c)
Figura B.3: Densidade de carga para a superfície (100) de GaAs, terminada em As, na reconstrução
2x1, em (a) e (b) vista de lado e (c) vista do topo.
75
Apêndice C
Propriedades eletrônicas da superfície (100)
de GaAs, terminada em As, na reconstrução
(2x1) com um átomo de N adsorvido
Na figura C.1, mostramos os nossos resultados para as bandas de energia da
superfície (100) de GaAs, terminada em As, na reconstrução (2x1) com um átomo de N
adsorvido sobre esta superfície. A área hachurada representa a projeção da estrutura
eletrônica do cristal perfeito. Nesta superfície, o sítio de adsorção considerado foi o sítio
C, pois é o sítio mais estável para a adsorção do N.
Analisando a figura C.1, observamos que a superfície (100) de GaAs, terminada
em As, na reconstrução (2x1), de acordo com os nossos resultados, apresenta dez
estados de superfície. O primeiro e o segundo estados desocupados, C
1
e C
2
, estão
aproximadamente a 0.61 e 1.09 eV, respectivamente, abaixo da banda de condução, ao
redor do ponto
K
. O estado C
3
está cerca de 1.87 eV abaixo da banda de condução ao
longo da direção 'JK − . O estado C
4
, por sua vez, situa-se a 2.14 eV abaixo da banda de
condução ao longo da direção 'JK − e, neste caso, não cruza em nenhuma direção com
o estado V
1
, que está a 1.53 eV acima da banda de valência. Deste modo, podemos
afirmar, baseados em nossos resultados, que a adsorção de um átomo de N sobre esta
superfície suprime o caráter metálico desta superfície. O estado V
2
está bem localizado
ao redor do ponto
K
e a 0.48 eV acima da banda de valência. O estado V
3
e V
4
da
superfície correspondem à interação entre os átomos da segunda e terceira camadas da
superfície. No “gap” interno, observamos apenas dois estados de superfície,
V
5
e V
6
,
com combinações do estado s localizadas nos átomos de As da superfície.
Na figura C.2, mostramos a densidade de carga calculada para a superfície (100)
de GaAs, terminada em As, na reconstrução (2x1) com um átomo de N adsorvido no
sítio C. Nesta figura, é possível observar que a concentração da densidade de carga é
76
maior, comparada com a densidade de carga calculada para a superfície (100) de GaAs,
terminada em As, na reconstrução (2x1) antes da adsorção do átomo de N, sobre o
átomo de N adsorvido, pois o N é mais eletronegativo que os átomos de As e Ga e sobre
o átomo de H, que satura a superfície, representada pelas linhas verdes. É possível,
também, observar uma certa inclinação dos dímeros de As, devido à presença do átomo
de N.
Figura C.1: Estrutura eletrônica para a superfície (100) de GaAs, terminada em As, na reconstrução 2x1
com um átomo de N adsorvido sobre esta superfície. A área hachurada representa a projeção da
estrutura eletrônica do cristal. As letras denotam os estados de superfície.
(a) (b)
(c)
Figura C.2: Densidade de carga para a superfície (100) de GaAs, terminada em As, na reconstrução
2x1,após a adsorção de um átomo de N, em (a) e (b) vista de lado e (c) vista do topo.
77
Referências Bibliográficas
[1] TRAGER-COWAN, C., “Report on the evening rump session on InN – July 21
st
,
2004 at the 2004 International Workshop on Nitride Semiconductors”, phys. stat.
sol. (c) v. 2, pp. 2240-2245, 2005.
[2] JAIN, S. C.,WILLANDER, NARAYAN, J., et al. “III-nitrides: Growth,
characterization, and properties”, J. Appl. Phys. v. 87, pp.965-1006, 2000.
[3] KARCH K., BECHSTEDT, F. “Ab initio lattice dynamics of BN and AnN:
Covalente versus ionic forces”, Phys. Rev. B v.56, pp.7404-7415, 1997.
[4] TSAO, J. Y., Materials Fundamentals of Molecular Beam Epitaxy. New York,
Academic Press, 1996.
[5] DESJONQUÈRES, M. C., SPANJAARD, D., Concepts in Surface Science. Berlin,
Springer-Verlag, 1996.
[6] DAVIES, JOHN., The Physics of Low-dimensional Semiconductors. New York:
Cambridge University Press, 1998, Stradling and Klipstein, Eds. Growth and
Characterization of Semiconductors. New York: Adam Hilger, 1990.
[7] STRADLING and KLIPSTEIN, Eds. Growth and Characterization of
Semiconductors. New York: Adam Hilger, 1990.
[8] POLITI, P. et al., “Instabilities in crystal growth by atomic or molecular beams”,
Physics Reports, v. 324, pp. 271-404, 2000.
[9] DUPUIS, R. R., DAPKUS, P. D., HOLONYAK, N. et al., “Room-temperature
laser operation of quantum-well Ga
(1-x)
Al
x
As-GaAs laser diodes grown by
metalorganic chemical vapor deposition” Appl. Phys. Lett. v. 32, pp.295-297, 1978.
[10] NAKAMURA, S., MUKAI, T., SENOH, M. “Candela-class high-brightness
InGaN/AlGaN double-heterostructure blue-light-emitting diodes”, Appl. Phys. Lett.
v. 64, pp. 1687 -1689, 1994.
[11] NAKAMURA, S., SENOH, M., NAGAHAMA, S., et al. ”Room-temperature
continuous-wave operation of InGaN multi-quantum-well-structure laser diodes
with a long lifetime” Appl. Phys. Lett., v.70, pp. 868 - 870 , 1997.
78
[12] YANG, H., ZHENG, L. X., LI, J. B., et al., “Cubic-phase GaN light-emitting
diodes” Appl. Phys. Lett. v. 74, pp. 2498-2500, 1999
[13] AS, D. J., RICHTER, A., BUSCH, J., et al., “Electroluminescence of a cubic
GaN/GaAs(001) p-n junction” Apll. Phys. Lett. v. 76, pp. 13-15, 2000
[14] GAMEZ-CUATZIN, H., TARDY, J. ROJO-ROMEO, P., PHILIPPE, A. et al.,
phys. stat. sol (a) v. 176, pp. 131, 1999.
[15] TANIYASU, Y., SUZUKI, K., LIM, D. H., et al. phys. stat. sol. (a) v. 180, pp.
241, 2000.
[16] ABERNATHY, C. R., MACKENZIE, J.D., BHARATAN, S. R., et al., “Growth of
In
x
Ga
1-x
N and In
x
Al
1-x
N on GaAs metalorganic molecular beam epitaxy” J. Vac.
Sci. Technol. A v. 13, pp. 716 - 718, 1995.
[17] ABERNATHY, C. R., MACKENZIE, J.D., DONOVAN, S. M., et al., “Dry
patterning of InGaN and InAlN” Appl. Phys. Lett. v. 64, pp. 3643-3645, 1994.
[18] VARTULI, C. B., PEARTON, S. J., LEE, J. W., et al., “lCl/Ar electron cyclotron
resonance plasma etching of III-V nitrides” Appl. Phys. Lett. v. 69, pp.1426 - 1428,
1996.
[19] LISCHKA, K. phys.stat. sol. (a) v. 177, pp. 135, 2000, AS, D. J., LISCHKA, K.,
phys. stat. sol. (a) v. 176, pp. 475, 1999.
[20] SMITH, A. R., FEENSTRA, R. M., GREVE, D. W., et al. “Reconstructions of
GaN (0001-bar) surface”, Phys. Rev. Lett. v. 79, pp. 3934 - 3937, 1997.
[21] SILVA, C. C., LEITE ALVES, H. W.,SCOLFARO, L. M. R., et al. “Pressure-
induced phase transitions and polytype structures in III-Nitrides”, phys. stat. sol.
(c), v. 2, pp. 2468 – 2471, 2005.
[22] BRANDT, O., YANG, H., JENICHEN, B., et al., “Surface reconstructions of zinc-
blende GaN/GaAs(001) in plasma-assisted molecular-beam epitaxy”, Phys. Rev. B
v. 52, pp. R2253 - 2256, 1995.
[23] YANG, H., BRANDT, O., WASSERMEIER, M., et al., “Evaluation of the surface
stoichiometry during molecular beam epitaxy”, Appl. Phys. Lett. v. 68, pp. 244 -
246, 1996.
[24] SCHIKORA, S., HANKELN, M., AS, D. J., et al., “Epitaxial growth and optical
transitions of cubic GaN films”, Phys. Rev. B v. 54, pp. R8381 – 8384 1996.
[25] FEUILLET, G., HAMAGUCHI, H., OHTA, K., et al., “Arsenic mediated
reconstructions on cubic (001) GaN”, Appl. Phys. Lett. v. 70, pp. 1025 - 1027,
1997.
79
[26] BERGER, A., TROOST, D., MÖNCH, W., “Adsorption of atomic nitrogen at
GaAs (110) surfaces”, Vacuum v. 41, pp. 669-671, 1990.
[27] TROOST, D., BAIER, H., -U., BERGER, A., et al., “Nitride layers on GaAs (110)
surfaces”, Surf. Sci. v.242, pp. 324 – 330, 1991.
[28] BRANDT, O., YANG, H., PLOOG, K., “Surface kinetics of zinc-blende (001)
GaN, Phys. Rev. B v. 54, pp. 4432 – 4435, 1996.
[29] NEUGEBAUER, J., ZYWIETZ, SCHEFFLER, M., et al., ”Clean and As- covered
zinc-blende GaN (001) surfaces: Novel surface structures and surfactant
behaviour”, Phys. Rev. Lett. v. 80, pp. 3097 – 3100, 1998; NEUGEBAUER, J.,
“Ab initio analysis of surface structure and adatom kinetics of group-III nitrides”,
phys. sta. sol.(b) v. 227, pp. 93 – 114, 2001; NEUGEBAUER, J.,”Surfactants and
antisurfactants on group-III-nitride surfaces, phys. stat. sol. (c) v. 0, pp. 1651 –
1667, 2003.
[30] GROSSE, F., GYURE, M. F., “Ab initio based modeling of III-V semiconductor
surfaces: Thermodynamic equilibrium and growth kinetic on atomic scales”, Phys.
Rev. B v. 66, pp. 075320: 1-13, 2002; e referências aí contidas.
[31] KRATZER, P., PENEV, E., SCHEFFLER, M., “First-principles studies of
kinetics in epitaxial growth of III-V semiconductor”, Appl. Phys. A v. 75, pp. 79 –
88, 2002; e referências aí contidas.
[32] KRATZER, P., PENEV, E., SCHEFFLER, M., “Reaction- limited island
nucleation in Molecular Beam Epitaxy of compound semiconductors”, Phys. Rev .
Lett. v. 88, pp. 036102: 1-4, 2002.
[33] ZANGWILL, A., Physics at Surfaces. UK, Cambridge University Press, 1992.
[34] NEWMAN, M. E. J., BARKEMA, G. T., Monte Carlo Methods in Statistical
Physics. UK, Clarendon Press, 2001.
[35] BINDER, K., HEERMANN, D. M., Monte Carlo Simulation in Statistical Physics.
Berlin, Springer-Verlag, 1997.
[36] FUCHS, M., SCHEFFLER, M., “Ab initio pseudopotentials for electronic
structure calculations of poly-atomic systems using density- functional theory”,
Comput. Phys. Commun. v. 119, pp. 67 – 98, 1999.
[37] GONZE, X., BEUKEN, J. –M., CARACAS, R., et al.,”First-principles
computation of material properties: the ABINIT software project”, Comput. Mater.
Sci. v. 25, pp. 478 – 492, 2002; e referências aí contidas.
80
[38] LABELLA, V. P., DING, Z., BULLOCK, D. W., EMERY, C., “A union of the
real-space and reciprocal-space view of the GaAs (001) surface”, Int. J. Mod. Phys.
B v. 15, pp. 2301 – 2333, 2001; e referências aí contidas.
[39] HOHENBERG, P., KOHN, W., “Inhomogeneous electron gas”, Phys. Rev. v. 136,
pp. B864 – B871, 1964.
[40] KOHN, W.,SHAM, L. J., “Self-consistent equations including exchange and
correlations effects”, Phys. Rev. v. 140, pp. A1133 – A1138, 1965.
[41] PAYNE, M.C., TETER, M. P., ALLAN, D. C., et al., ”Iterative minimization
techniques for ab initio total energy calculations: molecular dynamics and
conjugate gradients”, Rev. Mod. Phys. v. 64, 1045 -1097, 1992.
[42] PERDEW, J. P., ZUNGER, A., “Self-interation correction to density-functional
approximations for many-electrons”, Phys. Rev. B v. 23, pp.5048 - 5079, 1981.
[43] CEPERLEY, D. M., ALDER,B. J., “Ground state of the electron ga by a stochastic
method”, Phys. Rev.Letters v. 45, pp. 566 - 569, 1980.
[44] STAMPLF, C., VAN DE WALLE, C. G., “Density-functional calculations for III-
nitrides using the local-density approximation and generalized gradient
approximation”, Phys. Rev. B , v. 59, pp. 5521 - 5535 , 1999
[45] FUCHS, M., DA SILVA, J. L. F., NEUGEBAUER, J., et al., “Cohesive properties
of group-III nitrides: A comparative study of all-electron and pseudopotential
calculations using the generalized gradient approximation”, Phys. Rev. B , v. 65,
pp. 245212 - 245224, 2002
[46] ASHCROFT, N. W., MERMIN, N. D., “Solid State Physics”, Florida, Harcourt
College Publishers, 1976.
[47] TROULLIER, N., MARTINS, J. L., “Efficient pseudopotentials for plane-wave
calculations”, Phys. Rev. B v. 43, 1993 - 2006, 1991.
[48] STICK, I, PARRINELLO, R. CAR., BARONI, S., “Conjugate gradient
minimization of the energy functional: A new method for electronic structure
calculation”, Phys. Rev. B v. 39, 4997 -5004, 1989.
[49] BALERESCHI, A. “Mean-value point in the Brillouin zone”, Phys. Rev. B v.7, pp.
5212 - 5215, 1973.
[50] CHADI, D. J., COHEN, M. L., “Special points in the Brillouin Zone” Phys. Rev. B
v.8, pp. 5747 - 5753, 1973
[51] MONKHORST, H. J., PACK, J. D.,”Special points for Brillouin zone
integrations”, Phys. Rev. B v.13, pp. 5188 - 5192, 1976.
81
[52] HELLMANN, H., “Einführung in die quantenchemic” (Denticke, Leipzig), 1937.
[53] FEYNMAN, R. P., “Force in molecules”, Phys. Rev. v. 56, 340 - 343, 1939.
[54] KOONIN, S. E., Computational Physics. USA, Addison-Wesley Pub. Co, Inc,
1986.
[55] METROPOLIS, N., ROSENBLUTH, M., TELLER, A., et al., “Equation of state
calculations by fast computing machines”, J. Chem. Phys. v. 21, pp. 1087 - 1092,
1953.
[56] KORTLA, M., “Numerical simulations in the theory of crystal growth”, Comp.
Phys. Comm. V. 97, pp. 82 – 100, 1996.
[57] BORTZ, A. B., KALOS, M. H., LEBOWITZ, J. L., “A new algorithm for Monte
Carlo simulation of Ising spin systems”, J. Comp. Phys. v. 17, pp. 10 - 18, 1975.
[58] XIAO, R. F., ALEXANDER, J. D., ROSENBERGER, F., “Growth morphologies
of crystal surfaces”, Phys. Rev. A, v. 43, pp. 2977 - 2992, 1991.
[59] KRATZER, P., “Atomistic simulations of processes at surfaces”, In: Predictive
simulation of semiconductor processing status and challengs springer series in
materials science, v. 72 (Eds.) J. Dabrowski, E. R. Weber. Springer Verlag, 2004,
pp. 39 – 70. ISBN: 3 – 540 – 20481-4.
[60] RATSCH, C., RUGGERONE, P., SCHEFFLER, M., “Density functional theory of
surface diffusion and epitaxial growth of metals”, in Surface Diffusion: Atomistic
and collective processes, NATO ASI Series B, volume 360, Plenum Press, New
York, 1997.
[61] FICHTHORN, K. A., WEINBERG, “Theoretical foundations of dynamical Monte
Carlo Simulations”, J. Chemical Physics, v. 95, n
o
2, pp. 1090-1096, 1991.
[62] KREUZER, H. J., PAYNE, S. H., “Equilíbrio and dynamics of gas adsorption on
heterogeneous solid surfaces”, v. 104, pp. 153, Elsevier, Amsterdam, 1997.
[63] FICHTHORN, K. A., SCHEFFLER, M., “Island nucleation in thin-film epitaxy: A
first-principles investigation”, Phys. Rev. Lett., v.84, pp.5371 - 5374, 2000.
[64] FICHTHORN, K. A., MERRICK, M. L., SCHEFFLER, “A kinetic Monte Carlo
investigation of island nucleation and growth in thin-film epitaxy in the presence of
substrate-mediated interactions”, to be published in Applied Physics A.
[65] MURNAGHAN, F. D., “The compressibility of media under extreme pressures”,
Proc. Nat. Acad. Sci. USA v. 30, pp. 244 – 247 , 1944.
82
[66] FRITSCH, J., ECKL, C., PAVONE, P., et al., “Ab initio calculation of the
structure and dynamics of III-V semiconductor surfaces”, Festkörperprobleme,
Adv. In. Solid State Phys., v.36, pp. 135-162, 1997.
[67] TRUHLAR, D. G., GARRET, B. C., KLIPENSTEIN, S. J., “Current status of
transition-state theory”, J. Phys. Chem. v. 100, pp. 12771 - 12800, 1996.
[68] HÄNGGI, P., TALKNER, P., BORKOVEC, M., “Reaction-rate theory: fifty years
after Kramers”, Rev. Mod. Phys. v. 62, pp. 251 - 341, 1990.
[69] SRIVASTAVA, G. P., “Theory of semiconductor surface reconsturuction”, Rev.
Prog. Phys. v. 60, pp. 561 – 613, 1997.
[70] ALVES, J. L.A., HEBENSTREIT, J., SCHEFFLER, M., “Calculated atomic
structures and electronic properties of GaP, InP, GaAs and InAs (110) surfaces”,
Phys.Rev. B v. 44, pp. 6188 - 6198 , 1991.
[71] PASHLEY, M. D., “Electron counting model and its application to island
structures on molecular-beam epitaxy grown GaAs (001) and ZnSe (001)”, Phys.
Rev. B v. 40, pp. 10481 – 10487, 1989.
[72] NORTHRUP, J. E., FROYEN, S., “Structure de GaAs (001) surfaces: The role of
electrostatic interactions”, Phys. Rev. B v. 50, pp. 2015 – 2018, 1994.
[73] HASHIZUME, T., XUE, Q. K., ZHOU, Z., et al., “Structures of As-rich GaAs-
rich GaAs(001) – (2x4) Reconstructions”, Phys. . Lett. v. 73, pp. 2208 – 2211,
1994.
[74] MORIARTY, P., BETON, P. H., HENINI, M., “Sb-induced GaAs (111)B surface
reconstructions: success and failure of the electron-counting rule”, Surf. Sci. v. 365,
pp. L663 – L668, 1996.
[75] WHITMAN, L. J., THIBADO, P. M., ERWIN, S. C., et al., “Metallic III-V (001)
Surfaces: Violations of the electron counting model”, Phys. Rev. Lett. v. 79, 693 –
696, 1997
[76] BIEGELSEN, D. K., BRINGANS, R. D., NORTHRUP, J. E., et al., “Surface
reconstructions of GaAs(100) observed by scanning tunneling microscopy”, Phys.
Rev. B, .v. 41, pp. 5701 -5706 , 1990.
[77] LANDOLT-BÖRNSTEIN, New Series: Numerical Data and Functional
Relationships in Science and technology, 22a, ed. O. Madelung, Springer-Verlag,
Berlin, 1989.
[78] SCHMIDT, W., G. BECHSTEDT, ”Geometry and electronic structure of GaAs
(001) (2x4) reconstructions”, Phys. Rev. B, v.54, pp. 16742 - 16748, 1996.
83
[79] SCHMIDT, W. G., BECHSTEDT, F., “Atomic structures of GaAs(100)-(2x4)
reconstructions” Surf. Sci. v. 360, pp. L473, 1996.
[80] PAIVA, A. de, ALVES, J.L., LEITE ALVES, H., W. “Initial growth of GaN and
InN over GaAs (110) substrates”, braz. J. Phys., v. 34, pp. 608 – 610, 2004.
[81] ZHANG,S. B., NORTHRUP, J. E., “Chemical potential dependence of defect
formation energies in GaAs: application to Ga self-diffusion”, Phys. Rev. Letters,
v. 67, pp. 2339 – 2342, 1991.
Livros Grátis
( http://www.livrosgratis.com.br )
Milhares de Livros para Download:
Baixar livros de Administração
Baixar livros de Agronomia
Baixar livros de Arquitetura
Baixar livros de Artes
Baixar livros de Astronomia
Baixar livros de Biologia Geral
Baixar livros de Ciência da Computação
Baixar livros de Ciência da Informação
Baixar livros de Ciência Política
Baixar livros de Ciências da Saúde
Baixar livros de Comunicação
Baixar livros do Conselho Nacional de Educação - CNE
Baixar livros de Defesa civil
Baixar livros de Direito
Baixar livros de Direitos humanos
Baixar livros de Economia
Baixar livros de Economia Doméstica
Baixar livros de Educação
Baixar livros de Educação - Trânsito
Baixar livros de Educação Física
Baixar livros de Engenharia Aeroespacial
Baixar livros de Farmácia
Baixar livros de Filosofia
Baixar livros de Física
Baixar livros de Geociências
Baixar livros de Geografia
Baixar livros de História
Baixar livros de Línguas
Baixar livros de Literatura
Baixar livros de Literatura de Cordel
Baixar livros de Literatura Infantil
Baixar livros de Matemática
Baixar livros de Medicina
Baixar livros de Medicina Veterinária
Baixar livros de Meio Ambiente
Baixar livros de Meteorologia
Baixar Monografias e TCC
Baixar livros Multidisciplinar
Baixar livros de Música
Baixar livros de Psicologia
Baixar livros de Química
Baixar livros de Saúde Coletiva
Baixar livros de Serviço Social
Baixar livros de Sociologia
Baixar livros de Teologia
Baixar livros de Trabalho
Baixar livros de Turismo
