Categorias

Modelação Ab Initio da Difusão do Boro em Silício

Resumo

O presente trabalho tem como principal objectivo a exposição dos principais elementos inerentes à modelação ab initio de defeitos em estruturas cristalinas
de silício. Para tal foi usado o pacote de software AIMPRO, largamente usado na comunidade científica, incluindo no Departamento de Física da Universidade de Aveiro, cujo código foi desenvolvido por R. Jones e P.R. Briddon e tem sofrido várias contribuições de outros colaboradores.

Como aplicação foram recriadas várias conclusões teorico-experimentais relativas à difusão de átomos de boro em estruturas de silício cristalino e ligas SiGe, no mecanismo de difusão proposto recentemente por W. Windl et al. em 2001.

Foram confirmadas as conclusões de Cowern et al. (1994) e J. Bang et al. (2007) de que a presença de átomos de Ge na vizinhança dos átomos de B contribui para a retardação da difusão destes.

Introdução

Em 1854, o silício cristalino foi pela primeira vez sintetizado por Henri Deville e desde então tem se tornado dos materiais mais estudados pela comunidade científica. Tratando-se do segundo elemento mais abundante à face da Terra (a seguir ao oxigénio) permite que este seja bastante barato, comparativamente a outros semicondutores não tão abundantes.

As propriedades eléctricas dos materiais semicondutores são altamente influenciadas pela existência de defeitos e impurezas. Torna-se assim necessário compreender e controlar essas propriedades, de forma a adaptar os semicondutores às necessidades da indústria. Maioritariamente, esse controlo é feito por meio de dopagem com átomos de outros elementos diferentes dos átomos da rede.

O silício, elemento do grupo IV da tabela periódica, quando dopado com átomos de um elemento do grupo V, resulta num excesso de electrões livres, que se comportam como transportadores de carga, tornando o material tipo-n.

Aplicando a mesma analogia, quando dopado com átomos de um elemento do grupo III, resulta numa deficiência de electrões livres (o equivalente à formação de buracos) e assim torna o material tipo-p. A fronteira entre os materiais dos dois tipos é a conhecida junção p-n, o elemento crucial em todos os transístores e, inerentemente, na maioria dos dispositivos electrónicos. O boro é de longe o dopante tipo-p mais usado.

A última geração de chips informáticos apresenta transístores com apenas  e planos para  já estão em marcha. As regiões tipo-p e tipo-n são obtidas por métodos litográficos e implantação iónica. Contudo estes métodos, para além das limitações no que respeita à dimensão, provocam igualmente muitos danos na rede cristalina (intersticiais e aglomerados).

O boro substitucional difunde-se por meio de um mecanismo auxiliado por auto-intersticiais (Windl et al., 1999; Sadigh et al., 1999; Hakala et al., 2000; Jeong e Oshiyama, 2001). O boro intersticial é extremamente móvel, difundindo-se com uma energia de activação de  (Watkins, 1975) e assim torna-se móvel muito abaixo da temperatura ambiente (Tipping e Newman, 1987). Este fenómeno, em combinação com a geração significativa de intersticiais, resulta numa difusão significativa do boro substitucional. Este fenómeno é vulgarmente denominado transient enhanced diffusion (TED). Durante a TED o boro difunde-se ao longo de grandes distâncias. Stolk et al. (1995) observou a difusão de boro ao longo de distâncias na ordem dos , em meios saturados com intersticiais. Por difusão térmica (ou seja, na ausência de intersticiais) o boro difundiu-se apenas por . Assim, no processo de fabrico de , a ocorrência de tal fenómeno como a TED é desastrosa.

São conhecidas muitas técnicas que demonstram a redução eficaz da TED. Por exemplo a implantação de iões de baixa energia (Agardwal et al., 1997), co-implantação de boro e flúor (Downney et al., 1998), como também recozimentos térmicos rápidos (Mannino et al., 2001), entre outros. Normalmente uma ou mais combinações dos métodos anteriores são usados.

Contudo a barreira nanométrica dos processos de fabricos de chips não têm-se revelado a mais preocupante. Actualmente, as grandes empresas fabricantes de chips têm virado a sua atenção para o fabrico de chips cada vez mais rápidos e energeticamente eficientes. Em 2005 os primeiros chips de baseados em ligas SiGe foram produzidos em massa. Estes chips permitem um menor consumo de energia e velocidades de operação superiores a 200 GHz.

Com este mercado emergente, torna-se claro que é uma boa aposta investir no estudo destas ligas de forma a conhecer as suas características e propriedades. Assim, foi objectivo principal deste projecto estudar a difusão do boro em ligas de SiGe. Diversas experiências demonstram que a difusibilidade do boro em ligas SiGe decresce rapidamente, com o aumento da concentração de Ge (Kuo et al., 1993; Lever et al., 1998; Wang et al., 2004).

Capítulo I: O Método

1.1. O problema dos n-corpos.

Toda informação de um sistema mecânico quântico está contida na função de onda que obedece à equação de Shrödinger
(1.1)
onde é a energia total do sistema e o operador ℋ é o conhecido hamiltoniano. Este contém os termos cinéticos e potenciais de electrões e núcleos,
, (1.2)
onde e representam o número de electrões e núcleos, respectivamente, e , e representam a massa, carga e localização do núcleo α e as coordenadas do electrão i. Assim, a função de onda do sistema
, (1.3)
sendo o spin do electrão. Constata-se que Ψ é função de variáveis escalares, o que implica que a resolução dos problemas, mesmo os mais simples, é impraticável, mesmo com os computadores mais rápidos actualmente disponíveis. Torna-se assim necessário desenvolver um método, por meio de aproximações, tendo sempre presente a necessidade de manter a fiabilidade dos resultados.

1.2. Aproximação adiabática de Born-Oppenheimer


Como a massa do electrão é várias ordens de grandeza menor que a massa dos nucleões é razoável assumir que os electrões reagem instantaneamente ao movimento do núcleo. Assim, Born e Oppenheimer (1925) propuseram que a função de onda total seja dada por
(1.4)
onde e representam as funções de onda dos electrões e do núcleo, respectivamente, e o índice indica que a função de onda electrónica é determinada para uma dada configuração fixa de núcleos. 
Substituindo (1.4) em (1.1) obtém-se
(1.5)
, (1.6)
em que
(1.7)
As equações (1.5) e (1.6) apresentam todos os termos cinéticos ( ) e potenciais ( ) do operador hamiltoniano de (1.2). Os índices n e e denotam os termos nucleares e electrónicos.
No ponto de vista matemático, a aproximação adiabática considera o termo nulo. O estado estacionário do problema dos n-corpos é normalmente definido depois de encontrada a função para um conjunto de posições atómicas fixas.

1.3. Princípio variacional


Existem dois métodos principais para obter uma solução estacionária do hamiltoniano – métodos de elementos finitos e métodos variacionais. O primeiro, como o nome indica, integra a função de Schrödinger usando métodos discretos, e o segundo usa um princípio variacional. Neste método, é escolhido a dedo um subespaço do espaço de Hillbert para obter uma função aproximada da função de onda , da forma


. (1.8)
O principio variacional estabelece que o valor esperado para a energia total , funcional da função da função de onda tem valor mínimo,
. (1.9)
O valor aproximado da energia total será
, (1.10)
com e . Para o estado fundamental, a derivada de com deve-se anular, implicando que,
, para . (1.11)

1.4. Teoria de Hartree-Fock

Para se ter em atenção o Princípio de Exclusão de Pauli, o produto de funções em (1.8) deve ser anti-simétrico. Convenientemente, expressa-se este produto na forma de um determinante de Slater (Slater, 1929),
. (1.12)
O determinante garante a anti-simetria pois a troca de duas das spin-orbitas do electrão irão alterar a função num factor de enquanto a presença de duas spin-orbitas iguais resultam em .
Considerando como uma solução ao hamiltoniano do sistema (1.5) e (1.6), a minimização do valor da energia total, sujeita à condição de ortogonalidade , dá origem às equações de Hartree-Fock (Mc Weeney, 1989; Thijssen, 1999),
. (1.13)
O operador de Fock é dado por,
(1.14)
em que inclui a energia cinética dos electrões e a energia potencial de interacção electrões-núcleo. A energia potencial de interacção entre os electrões é traduzida pelos termos e ,
(1.15)
e

 .  (1.16)

O termo traduz o campo médio de Hartree e traduz o potencial de troca, que surge devido à anti-simetria das funções de onda.
A matriz na equação (1.13) é composta pelos multiplicadores de Lagrange resultantes do processo de minimização. Os elementos da diagonal são dados por,
(1.17)

que, após comparação com a equação (1.13), resulta em,
(1.18)
Normalmente, numa implementação computacional do método de Hartree-Fock, as spin-orbitas são expressas como uma combinação linear das orbitais atómicas. Após escolha cuidadosa do conjunto de funções , define-se
(1.19)
A equação de Fock toma a forma matricial, obtendo-se um problema de valores próprios generalizado, denominado equação de Roothaan (Roothaan, 1951),
(1.20)
onde é a matriz de dimensão definida por . A equação anterior é resolvida até a função de onda convergir, a par com a energia total. Koopmans (1934) demonstrou o significado físico dos valores próprios .

1.5. Teoria dos Funcionais da Densidade (DFT)


A Teoria dos Funcionais da Densidade (DFT) (Hohenberg e Kohn, 1964; Kohn e Sham, 1965), considera a energia total como uma funcional da densidade electrónica (equação (1.21)), em oposição à teoria HF em que considera individualmente as funções de onda electrónicas.


(1.21)
1º Teorema de Hohenberg-Kohn O potencial externo é um funcional único da densidade electrónica (Hohenberg e Kohn, 1964).
A energia total do sistema é dada por,
(1.22)
onde é o potencial externo sentido pelos electrões devido, por exemplo, às interacções núcleo-electrão, e é um funcional independente do sistema que toma em consideração e energia cinética electrónica, a energia de Hartree, a energia de correlação electrónica e a energia de troca electrónica.
2º Teorema de Hohenberg-Kohn Para uma densidade electrónica , obedecendo a e , a energia do estado fundamental toma valor mínimo (Hohenberg e Kohn, 1964),
. (1.23)
Demonstra-se assim a significativa simplificação que advém do uso da densidade de carga, pois este método, ao contrário da teoria HF, não apresenta aproximações. Contudo a contribuição de troca e correlação para o funcional mantêm-se não-local e é geralmente desconhecida.

1.5.1. Equações Kohn-Sham

A implementação da DFT transforma a equação de Schrödinger de muitos electrões (equação (1.4)) num sistema de equações de partícula única, que são as equações de Kohn-Sham (Kohn e Sham, 1965),
, (1.24)
em que a densidade é obtida a partir de (1.21).
Os primeiros três termos da equação (1.24) traduzem a energia cinética, o potencial externo imposto pelos iões (núcleos) e a energia de Hartree, respectivamente. O quarto termo agrega todos os efeitos dos corpos do sistema electrónico, sob a forma de um funcional de troca e correlação. A consolidação da DFT vem do facto de que existe uma densidade funcional universal (que depende exclusivamente da densidade electrónica ) que permite obter a solução exacta da densidade de carga e da energia total de um sistema no seu estado fundamental.
A energia total funcional para um sistema com muitos electrões é dada por
(Parr e Yang, 1989)
(1.25)
Na equação (1.25), a energia de Hartree e o potencial de troca e correlação ( e , respectivamente) são dados por
(1.26)
. (1.27)
As equações de Kohn-Sham (1.24) são resolvidas iterativamente, até se verificar convergência. Primeiramente escolhe-se uma densidade de carga inicial e resolve-se (1.24). Os vectores próprios resultantes servem como ponto de partida para a equação (1.21), gerando uma nova densidade de carga, repetindo-se novamente o ciclo, até que se verifique que a densidade de carga não se altera significativamente, sendo assim um método auto-consistente.
A forma exacta de é contudo desconhecida, mas algumas aproximações são actualmente usadas. De facto, a procura deste funcional continua a ser um dos tópicos de investigação actuais (Doren et al., 2001; Tao et al., 2008).
O formalismo anterior negligencia a dependência com o spin, aparte do termo . Considerando a densidade de carga como a sobreposição de duas densidades spin-up e spin-down, , o spin é contabilizado (von Barth e Hedin, 1972; Rajogopal e Callaway, 1973).

1.5.2. O funcional de troca e correlação

Uma aproximação comum para contornar o problema de não se conhecer é a chamada Aproximação à Densidade Local (LDA) e quando se considera o spin , Aproximação à Densidade de Spin Local (LSDA) (Kohn e Sham, 1965; von Barth e Hedin, 1972; Perdew e Zunger, 1981). A energia de troca e correlação é assumida como local e normalmente distinguida nas duas contribuições, a de troca e a de correlação :
, (1.28)
usando a notação da LSDA. A forma analítica do funcional de troca de um gás de electrões homogéneo é (von Barth e Hedin, 1972),
(1.29)
O termo correlativo é mais complexo pois envolve para um regime de alta densidade a teoria das perturbações (Perdew e Zunger, 1981), enquanto que para um regime de baixa densidade, usa o método de Monte Carlo aplicado a uma função de Green (Ceperly, 1978; Ceperly e Alder, 1980). Vários funcionais na forma parameterizada têm sido apresentados por Perdew e Zunger (1981) (PZ), Vosko et al. (1980) (VWN) e Perdew e Wang (1992) (PW).
Uma estimativa mais precisa da energia de troca e correlação é obtida considerando a expansão de primeira ordem da na densidade de carga, incluindo assim os termos dependentes com o gradiente da densidade de carga (Perdew, 1991; Perdew et al. 1996a,b). Este método de obter é conhecido como a Aproximação Generalizada ao Gradiente (GGA).

1.6. Pseudopotenciais

As propriedades químicas de um átomo são quase exclusivamente devidas aos electrões de valência e à sua interacção com os átomos, enquanto os estados do caroço são relativamente independentes do ambiente onde o átomo de encontra. Os electrões de valência em cada átomo interagem com o núcleo e os electrões do caroço via um potencial efectivo, chamado pseudopotencial.


O uso de pseudopotenciais tem a vantagem, comparativamente ao cálculo usando todos os electrões, de reduzir o número de electrões a tratar; uma vez que os electrões do caroço não são considerados.

1.7. Código Computacional AIMPRO

O pacote de software AIMPRO implementa alguns elementos da DFT que são importantes, contudo complexos. O entendimento do código não é o principal objectivo deste projecto, assim uma breve descrição de alguns métodos implementados pelo código AIMPRO são apresentadas nos apêndices desta dissertação.

1.8. Cálculo de observáveis


A DFT é uma teoria no estado fundamental do sistema. Assim, as propriedades dos estados excitados (incluindo níveis de energia não ocupados por electrões) estão fora do alcance dos cálculos da DFT. Contudo existem muitas propriedades de estado fundamentais que a DFT pode modelar com sucesso e cujos objectivos são prever ou confirmar observações experimentais. Esta secção explica como a maioria dessas propriedades importantes são derivadas.

1.8.1. Propriedades fundamentais

Estrutura local / Forças


Assim que as equações de Kohn-Sham são resolvidas para uma densidade de carga de uma configuração de átomos consistente é desejável saber que forças actuam em cada átomo. A força exercida num átomo é simplesmente igual a (Teorema de Hellman-Feynman) (Hellmann, 1937; Feymann, 1939) e assim, deslocando ligeiramente cada átomo da sua posição, e calculando , a força exercida em cada átomo pode ser calculada.
Variar a posição dos iões provoca uma variação da densidade de carga que se traduz de duas formas: uma variação em e uma alteração das funções de base . Assim que todas as contribuições de forem consideradas, sabe-se que, a partir da força que actua em cada átomo, os átomos poderão mover de acordo com essas forças. De uma forma iterativa, esse movimento é feito até que as forças e a variação na energia total em cada iteração sejam insignificantes. Terminado este processo, diz-se que a estrutura está relaxada. O AIMPRO faz uso de um algoritmo conjugado para o gradiente, o que significa que os átomos são movidos na direcção onde existia uma componente das forças obtidas na iteração anterior, bem como ao longo da componente das forças obtidas na iteração actual. A distância de que o átomo se desloca ao longo dessa direcção é escolhida, aproximando a uma função cúbica ou quadrática a função que expressa a variação da energia com a posição. Mas é importante salientar que tal aproximação simplista implica que, terminado o algoritmo, a estrutura relaxada estará localizada num mínimo local da superfície da energia total. Não garante assim que esta esteja localizada no mínimo global. Para garantir, na medida do possível, que a estrutura relaxada se encontra num mínimo global é necessário iniciar a optimização estrutural a partir de várias configurações iniciais diferentes.

Energia Total


Uma das propriedades fundamentais de um sistema é a energia total , equação (1.25). A energia total é a energia de uma super-célula ou agregado e é frequentemente usada para comparar a energia de sistemas semelhantes. Neste caso, um sistema similar contém o mesmo número de átomos de cada espécie (sendo cada espécie de átomo definida pelo pseudopotencial especifico) e a mesma carga global. Se estas condições não forem cumpridas é então necessário comparar as energias de formação. Já que a energia de formação tem em consideração a carga e o número de átomos de cada espécie diferente, através dos potenciais atómico e electrónico. A energia de formação é uma propriedade física bem definida.

1.8.2. Propriedades Derivadas

Energia de formação


O potencial químico de uma espécie de átomos é definida como a derivada da energia livre de Gibbs (Reif, 1965; Flynn, 1972),
. (1.30)

Considerando que se está perante o equilíbrio termodinâmico, é igual em toda a extensão do sistema, mesmo que se apresente diferenças de fase, assim é o equivalente à energia livre por partícula. Negligenciando o termo que é pequeno para as reacções no estado sólido e o termo que é uma aproximação apenas válida para baixas temperaturas, a energia de formação é dada por,
(1.31)

onde é a energia total do sistema, é o estado de carga do sistema, é o número de átomos da espécie e o potencial químicos dessa espécie. O potencial químico electrónico é usualmente,
(1.32)

onde é a energia de Fermi e é a energia da orbital ocupada de energia superior, habitualmente obtida a partir os níveis de Kohn-Sham. Alternativamente pode ser obtida, comparando e .
A energia de formação é uma quantidade que fornece muita informação pois, como tem em conta os potenciais químicos, é possível comparar super-células de diferentes tamanhos e comparar a estabilidade dos defeitos contendo diferentes números de espécies de átomos (útil para calcular energias de ligação, ver secção seguinte). Comparar a energia de formação de um dado defeito, para diferentes estados carregados, resulta na obtenção dos níveis eléctricos. A energia de formação pode também ser usada para calcular a solubilidade de um defeito.

Energias de ligação


Quando um defeito pode ser visto como a composição por dois ou mais defeitos primários, é normalmente útil saber a energia de ligação entre os constituintes. Poder-se-á calcular a energia de ligação de um complexo formado pelos dois constituintes e , comparando a energia de formação do complexo com as energias de formação dos constituintes,
(1.33)
onde e são as energias de formação dos dois constituintes que formam o complexo .
Um modo alternativo de calcular a energia de ligação é construir uma série de super-células em que os constituintes e estão inicialmente próximos (na forma do complexo ) e nas células seguintes estão gradualmente separados. Comparando as energias dessas super-células, podemos obter a energia como uma função da separação de e . Se a energia alcançar um limite assimptótico para a separação maior, a energia de ligação será a diferença entre a energia da super-célula contendo e separados e a energia da super-célula contendo o complexo . Este método é mais preciso que o método para determinar a energia de formação, especialmente quando as super-células estão igualmente carregadas, como é normalmente o caso, já que a tensão de primeira ordem e as interacções quadrupolares entre células vizinhas são iguais para cada separação e assim ocorrerá cancelamento das mesmas.

Energia de migração

O princípio fundamental para determinar a energia de migração de um átomo ou complexo é extremamente simples. O AIMPRO é capaz de implementar vários métodos para esse cálculo, contudo o método NEB (Nudged Elastic Band) foi o escolhido neste projecto.
O NEB permite identificar o caminho de energia mais baixo entre uma configuração inicial e uma configuração final. Esse caminho é muito importante no estudo de difusão dos defeitos e é normalmente denominado por caminho de energia mínima (MEP, do inglês minimum energy path).

Fig. 1 Gráfico de contornos de uma superfície potencial.

Observe-se a Fig. 1. Consideremos as configurações relaxadas e . Pretende-se determinar o MEP deste caso. Para iniciar o algoritmo que implementa o NEB no AIMPRO, é necessário definir-se configurações intermédias entre e , por exemplo, por simples interpolação linear – caminho a tracejado.
Este método considera que todas as configurações entre e , incluindo estas, estão ligadas por uma força com constante elástica . Iterativamente o NEB percorre todo o caminho inicial e avalia uma nova posição para as configurações intermédias (ou imagens) tendo em conta dois critérios: a) a força que actua perpendicularmente ao caminho entre as imagens e b) a força elástica que mantém coesas as configurações no caminho. Gradualmente estas operações levarão o método ao MEP desejado – caminho a traço cheio.

Fig. 2 Perfil do MEP.

A representação em perfil do MEP determinado pelo NEB (Fig. 2) permite observar uma barreira de difusão, sendo a altura dessa barreira a energia de migração do defeito. Contudo este método apresenta o mesmo problema, no que se refere ao mínimo local, como na secção onde se abordou a relaxação das configurações. Assim para garantir que o MEP adequado foi alcançado, deve-se repetir o processo para diferentes caminhos iniciais (ou seja, para diferentes imagens iniciais).
Estabilidade térmica
Muitas técnicas experimentais são capazes de observar a que temperatura um defeito desaparece, sendo assim desejável saber estimar a estabilidade térmica de um defeito. Assumindo que um complexo, quando sujeito ao tratamento térmico (annealing), se dissocia e não sofre qualquer reacção com outra espécie que se torna móvel à temperatura de annealing, torna-se simples estimar a estabilidade térmica. A taxa pela qual o complexo se dissocia é,
(1.34)
onde é uma dada frequência, é a constante de Boltzmann e é a temperatura. é a energia de activação do defeito que se dissocia, que é aproximadamente a soma da energia de formação do complexo e da energia de migração da espécie que se difunde, e que pode ser calculada usando o AIMPRO. A frequência em (1.34) é normalmente a frequência de Debye, outra propriedade que pode ser determinada por meio de modelação ab initio, contudo essa operação é complexa.

Níveis eléctricos


A energia de formação de defeitos carregados é uma função do potencial químico electrónico ou, equivalentemente, ao nível de Fermi. Assim, a energia de formação de um complexo ou defeito pode ser calculada como uma função de para diferentes estados carregados e, para um dado valor de , um dos estados apresentará um valor mais baixo, sendo assim o mais estável. O valor de para o qual dois estados carregados são iguais em energia, corresponde a um nível de ocupação. Por outras palavras, o valor de energia para o qual corresponde ao nível do defeito. Quando é maior que o seu valor crítico, o defeito estará no estado carregado e quando está no seu valor mais baixo, o defeito estará no estado carregado . Este método para calcular os níveis eléctricos é conhecido como o Método da Energia de Formação.
Contudo este método reduz frequentemente os níveis observados, por isso é habitual usar o Método do Marcador (Coutinho, 2002). Este método foi usado neste projecto e permite calcular os níveis eléctricos com uma melhor precisão, recorrendo à resolução da seguinte equação,
(1.35)
onde e são as energias do defeito no estados e , é a posição do nível do defeito. e são as energias do defeito marcador nos estados respectivos. Idealmente o defeito marcador deve ter propriedades eléctricas idênticas ao defeito em estudo e deve ter um nível cuja energia é conhecida. Conhecido o valor de e calculando , , e a equação (1.35) é resolvida e o valor de determinado. O Método do Marcador tem frequentemente um erro numérico associado que não costuma superar . Este método é contudo mais fiável devido ao cancelamento dos erros que advêm das interacções entre super-células. Calcular os níveis eléctricos é fundamental para identificar os defeitos electricamente activos.
Modos locais de vibração
Os modos de vibração de um cristal são dados a partir da matriz de dinâmica do sistema (Born e Huang, 1954) e resultam da resolução do seguinte problema de valores próprios,
, (1.36)
os valores próprios são as frequências quadradas associadas aos modos normais . Cada modo normal é um vector de dimensões que descreve a movimento de todos átomos desse modo. Os elementos da matriz são definidos por,
(1.37)
onde e são quaisquer coordenadas cartesianas e e são os desvios nessas direcções que os átomos de massa e sofreram, respectivamente.
A determinação das segundas derivadas na equação (1.37) é feita por cálculo numérico. A super-célula é primariamente relaxada, de forma que a força efectiva seja zero. O átomo é movido ligeiramente na proporção ( ) ao longo da direcção cartesiana . A densidade de carga desta nova configuração é então calculada e, a partir da nova densidade de carga calculada, as forças são daí derivadas. Note-se que a força efectiva já não é zero pois a estrutura foi perturbada. Considera-se como a componente da nova força que actua num átomo , na direcção cartesiana . O átomo é movido na mesma proporção mas na direcção oposta , resultando numa força . A 2ª derivada da energia é,
(1.38)
Esta forma de calcular a segunda derivada da energia inclui algumas contribuições não harmónicas. Por esta razão, as frequências obtidas deste modo são frequentemente referenciadas como frequências quasi-harmónicas (Jones et al., 1994).
Quando as impurezas, como no caso do boro em silício, são mais leves que os átomos da rede, os modos de vibração são superiores aos do cristal. Neste caso têm-se um modo local de vibração onde só a impureza e os seus primeiros vizinhos vibram. Por isso, só se calculam as segundas derivadas para a impureza e para os seus primeiros vizinhos.
Os modos de vibração de um defeito podem ser calculados e comparados com os resultados de espectroscopia do infra-vermelho ou fotoluminiscência. Como as frequências dos modos de vibração são extremamente sensíveis à estrutura de um defeito, se compararmos os resultados experimentais com os modelados, e se se verificar que são coerentes então há uma forte possibilidade de que a estrutura modelada e a observada sejam similares. Os avanços no estudo de impurezas isotópicas e efeitos de pressão têm conduzido a uma melhor identificação dos defeitos.

Capítulo II: Modelação Ab Initio

2.1 Introdução
O código AIMPRO começou a dar os primeiros passos à cerca de 20 anos nas mãos de R. Jones e P. R. Briddon, e desde então tem crescido e melhorado, muito em parte da sinergia entre os seus utilizadores. Com o auxílio deste software, muitos têm sido os estudos ab initio realizados. Muitos desses estudos contemplam o estudo de impurezas no silício cristalino, como por exemplo a dopagem de silício com boro. Neste domínio, destaca-se Adey (2004) pelo seu trabalho que também usou como ferramenta o AIMPRO.
Embora o estudo de defeitos de boro em silício já tenha sido realizado, poucos estudos da difusão deste foram publicados. O 1º mecanismo proposto, kick-out (Zhu et al., 1996) para a difusão do boro substitucional auxiliado por auto-intersticiais foi abandonado em favor de um mecanismo que sugere a difusão directa do boro (Windl et al., 1999). Considera-se que o mecanismo kick-out continua a ser válido, especialmente para temperaturas muito baixas, contudo o mecanismo de difusão directa proposto apresenta valores para a energia de migração muito próximos dos obtidos experimentalmente (Fahey et al., 1989; Haddara et al., 2000).
Experiências demonstram que em ligas SiGe a difusão do boro é retardada (Cowern et al., 1994). O objectivo é estudar a influência de átomos de Ge na difusão de boro em silício.
2.2 Método

Fig. 3 Super-célula cúbica de 64 átomos de Si.
O cristal foi simulado usando uma super-célula cúbica de 64 átomos de Si, que é suficiente para o estudo da difusão do boro (Windl et al., 1999). Utilizou-se o código computacional AIMPRO e todas as espécies químicas foram caracterizadas pelos respectivos pseudopodenticiais de Hartwigsen-Goedecker-Hutter (HGH) publicados por Hartwigsen et al. (1998).
2.3 Silício cristalino
Com carácter introdutório e de teste procedeu-se ao estudo do Si cristalino. Depois de relaxar a super-célula de silício e determinada a constante da rede , através do código AIMPRO, foram feitas algumas simulações de forma a determinar algumas propriedades importantes para o estudo do Si cristalino.
Equação de estado
O AIMPRO permite calcular os parâmetros que completam a equação de estado de
Birch-Murnaghan (Murnaghan, 1944; Birch, 1947),
. (2.1)
Na equação anterior, representa a volume da super-célula, o módulo de compressibilidade (bulk modulus) e é a sua primeira derivada com a pressão.
Tabela 1 Resultados dos parâmetros calculados para o Si cristalino. A constante de rede a0 é comparada com o valor experimental de James e Lord (1992). O módulo de compressibilidade é comparado com o valor experimental de Singh (1993) e a sua derivada com Beattie e Schirber (1970).
Parâmetros calculados experimental
(Å)
5.40 5.43
(GPa)
94.60 97.9
(u.a.)
3.94 4.09
Estrutura de bandas

Fig. 4 Estrutura de bandas do Si cristalino. Os pontos G, L1, L2, L3, L4, Y, X e Z correspondem aos pontos-k de coordenadas (0 0 0), (¼ ¼ ¼), (-¼ ¼ ¼), (-¼ -¼ ¼), (¼ -¼ ¼), (0 ½ ½), (½ 0 0) e (0, 0, ½), respectivamente.
Na Fig. 4 observa-se o gap indirecto do silício, de acordo com os dados experimentais, e a subestimação da energia do gap em 40%, inevitável no uso da LDA.
2.4 Boro substitucional

Fig. 5 Estrutura do boro substitucional Bs.
À super-célula usada nos cálculos para o Si cristalino, foi-lhe removida um átomo central e no seu lugar foi adicionado um átomo de 10B. Após relaxamento obteve-se a estrutura da Fig. 5. Os parâmetros mais importantes que definem a estrutura são apresentados na Tabela 1. Note-se que como o átomo de boro tem um raio menor, a relaxação da estrutura provocou uma ligeira compressão do centro.
Tabela 2 Parâmetros calculados importantes da configuração do boro substitucional nas diferentes configurações de carga.
carga -1 0 +1

2.2014 2.0584 2.0539

2.0513 2.0696 2.0766

111.062° 109.775° 110.107°

107.834° 109.170° 108.832°
Usando espectroscopia de infravermelho, observou-se modos locais de vibração do Bs a e para 10Bs e 11Bs, respectivamente (Smith e Angress, 1963; Angress et al., 1965; Balkanski e Nazarewicz, 1966; Bean et al., 1972). Das simulações levadas a cabo no AIMPRO, verificou-se os mesmos modos a e respectivamente, observando-se assim um desvio de dos valores experimentais. Os experimentais e os calculados apresentam o mesmo desvio isotópico de .
À semelhança no que se fez no caso do silício cristalino, procedeu-se à determinação da estrutura de bandas para a configuração Bs0.

Fig. 6 Estrutura de bandas para o boro substitucional. No estado neutro, as linhas a azul correspondem a níveis ocupados, os linhas a tracejado representam os níveis vazios e a linha a vermelho o nível aceitador. O caminho de pontos-k é idêntico ao usado para obter a Fig. 4.

2.5 Difusão do Boro em Silício
No mecanismo directo de difusão do boro positivo, proposto por Windl, está esquematicamente representado sequencialmente de a) a m), na Fig. 7.

a)
b)
c)

d)
e)
f)

g)
h)
i)

j)
k)
l)

m)
Fig. 7 Mecanismo directo de difusão do boro positivo em silício.
O método NEB aplicado às configurações da Fig. 7 permitiu obter o MEP deste mecanismo.
O perfil do MEP (Fig. 8) permite identificar a energia de migração deste defeito em , em concordância com o valor obtido por Windl (1999).

Fig. 8 Perfil do MEP da difusão do B positivo em Si.
2.6 Difusão do Boro em Ligas Silício-Germanio
O estudo da secção anterior foi aplicado a diversas super-células contendo 1 ou 2 átomos de Ge colocados em posições estratégicas na rede. Os resultados são apresentados nas Fig. 9 e Fig. 10.

barreira de migração = 1.3262
barreira de migração = 1.1857
barreira de migração = 1.1945
Fig. 9 Difusão do boro numa super-célula de Si decorada com um átomo de Ge. Os átomos a verde são os átomos de boro, a azul os de germânio e a laranja os de silício.

Nas restantes configurações onde o átomo de Ge ocupa uma posição distante do defeito (Fig. 9, ao centro) e em que o átomo de Ge é primeiro vizinho do átomo de boro (Fig. 9, à direita), nenhuma destas apresenta uma alteração significativa na barreira de migração.

barreira de migração = 1.3108
barreira de migração = 1.3238
barreira de migração = 1.1798

barreira de migração = 1.1640
barreira de migração = 1.2066
Fig. 10 Difusão do boro numa super-célula de Si decorada com 2 átomos de Ge. O mesmo código de cores da Fig. 9 que permite identificar as espécies de átomos é usado nesta figura.
Novamente, a barreira de migração é alterada significativamente quando o Ge ocupa a posição intersticial junto do boro. Verifica-se que também que a presença de átomos de Ge afastados do defeito não faz alterar a barreira de migração. No caso em que os átomos de Ge estão próximos do defeito, a barreira de migração sofre um aumento, contudo não tão significativo como no Ge intersticial.

Capítulo III: Conclusões

A realização deste projecto permitiu adquirir técnicas de modulação importantes para o estudo dos materiais cristalinos.

Calculou-se as propriedades estruturais e electrónicas do silício na fase de diamante, nomeadamente i) a constante de rede, o módulo de compressibilidade e a sua primeira derivada relativa à pressão, e ii) estrutura de bandas, que apresenta um gap indirecto, corroborando bem com os resultados experimentais, exceptuando para a largura do gap, que é subestimada, resultado da aproximação à densidade local (LDA).

O boro substitucional foi modelado. As suas frequências de vibração e o seu desvio isotópico concordam significativamente com os valores experimentais. A estrutura de bandas apresenta um nível aceitador, rasante ao máximo da banda de valência.

Aplicou-se o método NEB (Nudged Elastic Band) à difusão do boro positivo, assistida por um auto-intersticial de silício. Os valores da barreira de migração de  estão de acordo com os valores de Windl et al. (2001). Para alguns centros Bs-Sii, decorado com átomos de germânio, a barreira de migração aumentou de . Este centro poderá retardar a difusão do boro, conforme se verifica experimentalmente.

A difusão do boro em silício amorfo é muito diferente da difusão em silício cristalino (Mirabella et al., 2008). Poucos estudos sobre este tópico foram publicados. Este problema apresenta-se naturalmente num futuro projecto.