Transcript Modelli cinetici di trasporto

4
METODO MONTE CARLO
4
60
Metodo Monte Carlo
L’equazione di Boltzmann descrive il trasporto di cariche elettriche nei semiconduttori. A causa della complessità degli operatori di collisione e del fatto
che l’equazione di Boltzmann è accoppiata in modo non lineare all’equazione
di Poisson, un approccio numerico basato sul metodo delle differenze finite
o degli elementi finiti è molto complicato, anche se negli ultimi anni notevoli
passi avanti sono stato fatti in questa direzione.
Uno degli approcci numerici alternativi che si è dimostrato di maggiore
efficacia nell’ affrontare il problema della risoluzione dell’equazione di Boltzmann è, senza dubbio, il metodo Monte Carlo che consiste, essenzialmente,
nel calcolare stocasticamente il tempo di volo balistico dei portatori di carica ed il loro stato finale dopo una collisione. Questo metodo non presenta
eccessive difficoltà di calcolo ed è di immediata interpretazione fisica. Gli
svantaggi del metodo Monte Carlo sono i tempi di calcolo piuttosto lunghi ed
i problemi di rumore statistico. Questi svantaggi lo rendono un metodo poco
pratico ai fini della progettazione e dell’analisi dei dispositivi a semiconduttore, tuttavia esso rappresenta un termine di confronto imprescindibile ogni
qual volta si voglia valutare la bontà di un nuovo modello di trasporto.
4.1
Metodo Monte Carlo a singola particella
Il trasporto di cariche in un semiconduttore è un problema a molti corpi.
Tuttavia, quando le cariche si possono considerare indipendenti, è possibile
usare un metodo approssimato che descrive l’ ‘ensemble’ di cariche seguendo
la storia temporale di un singolo portatore che subisce molti processi di
scattering: è questo il metodo Monte Carlo a particella singola.
Esso è utile nello studio delle proprietà di trasporto delle cariche in un
blocco omogeneo di semiconduttore.
Nel seguito consideremo il caso in cui il portatore di carica sia un elettrone. Il metodo a singola particella simula il moto di un elettrone nello
spazio dei vettori d’onda: ‘spazio degli impulsi’. Il vettore d’onda iniziale
viene scelto in maniera stocastica, dopo di che l’impulso del’lelettrone, supponendo che siano soddisfatte condizioni tali che sia valida la descrizione
semiclassica, evolve secondo l’equazione
d k(t)
q
= ∇x V (x(t)),
dt
~
da cui
Z
k(t + τ ) − k(t) =
t
t+τ
q
∇x V (x(s)) ds
~
(1)
(2)
4
METODO MONTE CARLO
61
che nel caso in cui il campo elettrico sia un campo elettrico E uniforme
applicato dall’esterno diventa
q
k(t + τ ) − k(t) = − Eτ.
~
(3)
Fino a questo momento abbiamo ragionato come se il moto di un elettrone
all’interno di un cristallo fosse libero: volo libero. Ma sappiamo che un
elettrone è continuamente soggetto a scattering che ne cambiano istantaneamente il vettore d’onda. Il volo libero si ha quindi solo tra due scattering
consecutivi.
Questo volo libero avrà quindi una durata τ stocastica che dipende dalla
probabilità totale per unità di tempo che un elettrone sia scatterato, P (k),
definita da
Z
P (k) =
P (k0 , k) dk0 .
(4)
B
Chiamiamo W (t) la probabilità che un elettrone non subisca urti nell’intervallo di tempo [0, t]. Gli urti sono processi Markoviani, il che significa che la
probabilità che un elettrone subisca un urto ad un certo istante dipende solo
dallo stato dell’elettrone in quell’istante e non dalla sua storia precedente.
Pertanto, essendo P (k(t)) dt la probabilità che un elettrone che ha vettore d’onda k(t) all’istante t subisca un urto nell’intervallino di tempo [t, t +
dt], si ha
W (t + dt) = W (t) (1 − P (k(t)) dt) ,
da cui
W (t + dt) − W (t) = −W (t)P (k(t)) dt.
Dividendo per dt e passando al limite, otteniamo l’equazione differenziale
ordinaria
d
W (t) = −W (t)P (k(t)).
(5)
dt
Questa equazione può essere integrata immediatamente, con l’ovvia condizione iniziale W (0) = 1. La soluzione è
Z t
W (t) = exp −
P (k(s)) ds .
(6)
0
Usando questa formula possiamo esprime la durata di volo libero. Si procede
nel modo seguente:
1. si sceglie in maniera random un numero r1 uniformemente distribuito
tra 0 ed 1;
4
METODO MONTE CARLO
62
2. si risolve rispetto a τ l’equazione
Z
r1 = exp −
τ
0
0
P (k(t )) dt
.
(7)
0
L’inversione dell’equazione (7) è molto costosa in termini di calcolo,
perché richiede la conoscenza della traiettoria k(t). Per superare questa
difficoltà si procede in un altro modo (metodo della costante Γ).
Sappiamo che la probabilità P (k) è la somma di tanti termini, rappresentanti le probabilità per ogni fenomeno di scattering. Se supponiamo di
tener conto di N fenomeni di scattering, possiamo scrivere
P (k) =
N
X
Pj (k).
(8)
j=1
Per calcolare τ da (7), introduciamo un ulteriore scattering virtuale P0 (k),
detto self-scattering, durante il quale l’elettrone conserva inalterato il suo
stato. La probabilità di self-scattering è definita in modo tale che
N
X
Pj (k) = Γ,
(9)
j=0
per una certa costante Γ. Questa costante viene scelta in modo da soddisfare
la disuguaglianza
N
X
Γ ≥ sup
Pj (k),
(10)
k
j=1
la quale assicura che
P0 (k) = Γ −
N
X
Pj (k) ≥ Γ − sup
k
j=1
N
X
Pj (k) ≥ 0.
j=1
Inoltre Γ deve essere quanto più piccola possibile, in modo da minimizzare
il numero di self-scattering. Osserviamo che per soddisfare (10) è necessario
conoscere le probabilità totali di scattering, ma non la traiettoria di un dato
elettrone.
Con questo artificio, usando (9), τ risulta individuato dall’equazione
r1 = e−Γτ ,
4
METODO MONTE CARLO
63
da cui si ricava la semplice formula
τ =−
ln r1
.
Γ
(11)
Dopo un tempo di volo libero determinato dalla formula precedente,
la particella subirà un certo processo di scattering, che deve essere scelte
anch’esso in modo stocastico.
Per determinare il processo di scattering, si procede nel modo seguente:
1. si definiscono le quantità
Λn (k) =
n
X
Pj (k
j=1
Γ
,
n = 1, . . . N,
(12)
che sono le somme successive delle probabilità di scattering normalizzate a Γ, e si pone Λ0 (k) = 0;
2. si genera un numero random r2 uniformemente distribuito tra 0 e 1, e si
sceglie come meccanismo di scattering quello corrispondente all’indice
n tale che
Λn−1 (k) ≤ r2 < Λn (k);
(13)
3. se risulta
ΛN (k) ≤ r2 ,
allora si sceglie il self-scattering, cioè non succede niente e si genera
un altro numero random.
Si noti che in questo modo il principio di esclusione di Pauli non viene preso
in considerazione.
Il problema successivo è quello di individuare il vettore d’onda dell’elettrone dopo lo scattering.
Il modulo |k0 | del vettore d’onda k0 dopo l’urto è noto per la conservazione dell’energia1 , quindi bisogna determinare solo la sua direzione. Per
far ciò, è sufficiente determinare due angoli θ e ϕ rispetto ad un sistema
di assi (x̂, ŷ, ẑ) solidale con il cristallo. Più precisamente, indichiamo con θ
l’angolo formato dal versore k̂0 = k0 /|k0 | con il versore ẑ,
cos θ = k̂0 · ẑ,
1
0 ≤ θ ≤ π,
Per l’energia della banda di conduzione si usa l’approssimazione parabolica o di Kane.
4
METODO MONTE CARLO
64
e con ϕ l’angolo formato dalla proiezione sul piano x̂ŷ di k̂0 con il versore x̂,
cos ϕ =
[k̂0 − (k̂0 · ẑ)ẑ] · x̂
|k̂0 − (k̂0 · ẑ)ẑ|
=
k̂0 · x̂
,
| sin θ|
0 ≤ ϕ ≤ 2π.
Se lo scattering è isotropo, cioè se dopo l’urto la probabilità di tutte
le direzioni è la stessa, la probabilità che l’elettrone abbia una direzione
individuata da angoli compresi tra ϕ e ϕ + dϕ, θ e θ + dθ è proporzionale
all’elemento di angolo solido
dΩ = sin θ dθ dϕ ≡ −d cos θ dϕ.
Quindi possiamo considerare cos θ e ϕ come variabili stocastiche omogenee.
Per quanto visto, gli angoli si individuano nel seguente modo:
1. si scelgono in maniera random due numeri r3 ed r4 uniformemente
distribuiti tra 0 e 1;
2. a partire da r3 , si prende ϕ come determinazione di una variabile
random uniformemente distribuita tra 0 e 2π, ovvero
ϕ = 2πr3 ;
3. a partire da r4 si prende θ in modo che cos θ sia la determinazione di
una variabile random uniformemente distribuita tra −1 e 1, ovvero
cos θ = 1 − 2r4 .
In conclusione, il numero d’onda di un elettrone dopo uno scattering
isotropo avrà componenti
 0
k = |k| sin θ cos ϕ,


 x
ky0 = |k| sin θ sin ϕ,
(14)


 0
kz = |k| cos θ.
Più complicata è la situazione nel caso in cui lo scattering sia anisotropo.
In tal caso, la probabilità di transizione per unità di tempo dipenderà
dall’angolo θ0 tra le direzioni dell’elettrone prima e dopo l’urto.
Ad esempio, nel caso di urto con le impurezze, la probabilità di scattering
Pimp (k, k0 ) dipende da |k − k0 |, e lo scattering è quindi anisotropo. Per calcolare la probabilità totale di scattering con le impurezze, Pimp (k), occorre
effettuare un’integrazione rispetto a k0 . Per calcolare questo integrale usiamo
4
METODO MONTE CARLO
65
l’approssimazione parabolica, per cui la zona di Brillouin viene rimpiazzata
da tutto R3 , e introduciamo coordinate sferiche (|k0 |, ϕ0 , θ0 ) rispetto ad un
sistema cartesiano (x̂0 , ŷ0 , ẑ0 ) solidale con k, in cui l’asse ẑ0 coincide con k̂.
Gli angoli ϕ0 e θ0 , con 0 ≤ ϕ0 ≤ 2π, 0 ≤ θ0 ≤ π, sono definiti in modo analogo
a ϕ e θ. Otteniamo cosı̀
Z
Pimp (k) =
Pimp (k, k0 ) dk0
3
R
Z ∞ Z 2π Z π
Pimp (k, k0 )|k0 |2 sin θ0 dθ0 dϕ0 d|k0 |.
=
0
0
0
Adesso siamo in grado di definire la probabilità (non normalizzata) che l’angolo tra k e k0 sia compreso tra 0 e θ0 , semplicemente cambiando gli estremi
d’integrazione rispetto alla variabile θ0 :
Z ∞ Z 2π Z θ0
0
Pimp (k, k0 )|k0 |2 sin θ0 dθ0 dϕ0 d|k0 |.
Pimp (k, θ ) =
0
0
0
La probabilità (normalizzata) che l’angolo tra k e k0 sia compreso tra 0 e θ0 ,
è data da
Pimp (k, θ0 )
.
Pimp (k)
Questa funzione, che è una funzione stocastica omogenea, può essere calcolata esplicitamente, ottenendo
2 2|k|
0
(1 − cos θ ) 1 + Kimp
Pimp (k, θ0 )
=
2 .
Pimp (k)
2|k|
0
2 + (1 − cos θ ) Kimp
In definitiva, gli angoli ϕ0 , θ0 si individuano nel seguente modo:
1. si scelgono in maniera random due numeri r5 ed r6 uniformemente
distribuiti tra 0 e 1;
2. a partire da r5 , si prende ϕ0 come determinazione di una variabile
random uniformemente distribuita tra 0 e 2π, ovvero
ϕ = 2πr5 ;
3. a partire da r6 si prende θ0 in modo che Pimp (k, θ0 )/Pimp (k) sia la
determinazione di una variabile random uniformemente distribuita tra
0 e 1, ovvero
Pimp (k, θ0 )
= r6 ,
Pimp (k)
4
METODO MONTE CARLO
66
da cui si ricava
2r6
cos θ0 = 1 −
1 + (1 − r6 )
2|k|
Kimp
2 .
Gli angoli ϕ0 e θ0 ci permettono di determinare le coordinate di k0 in una
terna che ha k̂ come asse ẑ0 . È conveniente scegliere come asse x̂0 la linea
dei nodi, cioè l’intersezione tra il piano x̂ŷ di riferimento e il piano passante
per l’origine ed ortogonale a k. Allora, usando il formalismo degli angoli di
Eulero, è possibile eprimere le componenti di k0 nella terna di riferimento
mediante due opportune rotazioni.
Usando il metodo descritto fin qui, possiamo determinare il moto di un
elettrone fino ad un tempo massimo fissato, tmax . Per ergodicità si può
assumere che se il tempo di volo è sufficientemente lungo (dell’ordine di
2 · 10−9 s), le medie temporali sono equivalenti alle medie d’insieme. Quindi,
calcolando le medie temporali delle grandezze fisiche si hanno informazioni
sui valori medi di queste grandezze per gli elettroni.
La velocità istantanea dell’elettrone è
v(k) =
1
∇k E(k),
~
quindi la velocità media in un volo libero è
hviτ =
1 ∆E
,
~ ∆k
(15)
dove il simbolo ∆ indica incrementi della grandezza nel tempo τ e l’espressione s’intende valida componente per componente. Ora, essendo
∆k = −
qe
Eτ,
~
dove E è il campo elettrico, supposto uniforme, si ricava
hviτ = −
∆E
.
qe Eτ
(16)
Quest’ultima equazione si intende valida componente per componente.
Per quanto visto, la velocità media durante il tempo di simulazione totale
tmax è
1 X
1 X
∆E
hvi =
hviτ τ =
−
τ,
tmax
tmax
qe Eτ
v.l.
v.l.
4
METODO MONTE CARLO
67
dove la sommatoria è su tutti i voli liberi dell’elettrone. In conclusione si
ricava la formula
X
1
hvi = −
∆E.
(17)
qe Etmax
v.l.
Se introduciamo l’energia dell’elettrone all’inizio di un volo libero, Ei ,
e l’energia dell’elettrone alla fine del volo libero, Ef , la formula precedente
diventa
X
1
hvi = −
(Ef − Ei ).
(18)
qe Etmax
v.l.
Infine, approssimando l’energia media durante un volo con
Ei + Ef
,
2
(19)
1 X Ei + Ef
τ.
tmax
2
(20)
hEiτ =
si ha
hEi =
v.l.