Równanie Schr¨odingera w polu magnetycznym

Jeżeli cząstka o ładunku q znajduje się w polu magnetycznym B(r) = ∇ × A(r), gdzie A — potencjał wektorowy pola, to jej hamiltonian wyraża się poprzez uogólniony pęd P = −i¯h∇ = p+(q/c)A. Nierelatywistyczne równanie Schr¨odingera dla elektronu (cząstki o spinie 1/2 i ładunku q = e = −|e|) zapisujemy jako równanie Pauliego [163–166]

" elektronów. W tym przypadku dogodnie jest obrać jako jednostki długości, czasu i energii wielkości związane z polem magnetycznym:

l_u=

gdzie B_u — jednostkowe pole magnetyczne (np. wartość średnia lub charakterystyczna dla układu). W tych jednostkach otrzymujemy parametr skali α = 1, a magneton Bohra µ_B = |e|¯h/(2mc) redukuje się do jedności.

Odpowiednikiem stacjonarnego bezwymiarowego równania masy efektywnej (3.17) jest



[−i∇ + A(r)] 1

m(r)[−i∇ + A(r)] + V (r) + g(r)B · σ



ψ(r) = αεψ(r), (3.21) przy czym w niejednorodnym materiale efektywny czynnik Land´ego g może być funkcją położenia.

Zmianie musi ulec także definicja wektora gęstości prądu prawdopodobieństwa (3.18), gdyż jest on związany z operatorem prędkości cząstki v = (P + A)/m:

j(r, t) = 2

αm(r)<{Ψ^∗(r, t)[−i∇ + A(r)]Ψ(r, t)}. (3.22) Dalej zajmować się będziemy szczególnym przypadkiem równania (3.21), kiedy poten-cjał V , pole B, masa efektywna m i czynnik g zależą od jednej tylko współrzędnej (x),

a kierunek pola magnetycznego jest stały w przestrzeni (B = (0, 0, B(x)); stąd automa-tycznie ∇ · B = 0). Równanie to opisuje wówczas układ kwazijednowymiarowy, którego lokalne właściwości zmieniają się wzdłuż jednego wyróżnionego kierunku (osi X). Wpro-wadzając potencjał wektorowy w postaci

A= [0, A(x), 0], A(x) = względu na translacyjną niezmienniczość układu w płaszczyźnie YZ można napisać

ψ_σ(r) = ψ_σ(x)e^i(kyy+kzz),

a wyrazy pochodzące od pola magnetycznego i ruchu w płaszczyźnie YZ potraktować jako dodatkowy efektywny potencjał zależny tylko od x:

V_ef(x) = [k_y + A(x)]²+ k_z²

m(x) ± g(x)B(x). (3.24)

Składowe prądu prawdopodobieństwa dane są teraz formułami j_x,σ(x) = 2

Przedstawione i wyprowadzone powyżej równania posłużą w kolejnym rozdziale do sfor-mułowania algorytmów numerycznych pozwalających obliczać przewodnictwo rozważanych układów niskowymiarowych.

4

Metody obliczeń

Rozdział ten przedstawia dwie metody stosowane przy rozwiązywaniu równania Schr¨o-dingera i równania masy efektywnej: metodę różnic skończonych i metodę elementów skoń-czonych. Opisuje również analityczny formalizm map śladów–antyśladów wyprowadzony z tej ostatniej dla aperiodycznych supersieci. Opisane jest także zastosowanie metody róż-nic skończonych do wyznaczania rozkładu potencjału elektrostatycznego w strukturze, do której przyłożono zewnętrzne napięcie.

4.1. Metoda różnic skończonych

Podrozdział ten opisuje metodę rozwiązywania jednowymiarowego stacjonarnego rów-nania Schr¨odingera i rówrów-nania masy efektywnej za pomocą różnic skończonych (metodę siatkową) oraz jej zastosowanie do stanów rozproszeniowych.

4.1.1. Dyskretne jednowymiarowe równanie Schr¨odingera

Zagadnienia kwantowomechaniczne, dla których równanie Schr¨odingera (lub równa-nie masy efektywnej) posiada dokładne rozwiązania analityczne, są stosunkowo równa-nieliczne.

W większości przypadków konieczne jest stosowanie metod przybliżonych, w tym metod numerycznych.

Szeroką i stosunkowo uniwersalną klasę metod numerycznych rozwiązywania równania Schr¨odingera stanowią metody różnic skończonych (metody siatkowe) [175–178]. Najprost-szą postać przyjmują one dla zagadnień jednowymiarowych.

Niech potencjał V (x) będzie zadany na przedziale ha, bi. Zdefiniujmy na tym przedziale siatkę N + 1 równoodległych punktów

x_i = a + i b − a

N + 1 = a + is (i = 0, . . . , N), (4.1)

gdzie s — krok siatki. W punktach tych określone będą wartości potencjału V_i = V (x_i) oraz funkcji falowej ψ_i = ψ(x_i).

Stosując do jednowymiarowego stacjonarnego równania Schr¨odingera

"

− d²

dx² + αV (x)

#

ψ(x) = αε(x), (4.2)

37

będącego szczególnym przypadkiem równania (3.23), trójpunktowe przybliżenie drugiej pochodnej (z dokładnością do wyrazów małych rzędu s²) [179]

d²

otrzymujemy dyskretną postać stacjonarnego jednowymiarowego równania Schr¨odingera

−ψi−1+ (2 + αs²V_i)ψ_i− ψi+1 = αs²εψ_i (i = 1, . . . , N). (4.3) Równanie to jest podstawą różnych podejść numerycznych, takich jak metoda algebraiczna i metoda strzałów [176, 178].

Jeżeli poszukujemy stanów związanych, których funkcje falowe muszą zanikać przy x → ±∞, wówczas możemy przyjąć warunki brzegowe

ψ(a) = ψ₀ = ψ(b) = ψ_{N +1} = 0. (4.4)

Warunki te są w ogólnym przypadku tylko przybliżone, jednak dokładność tego przybliże-nia jest wystarczająca, o ile przedział ha, bi jest wybrany odpowiednio. Uwzględprzybliże-niając je, możemy zapisać układ równań (4.3) dla i = 1, . . . , N w postaci algebraicznego zagadnienia własnego z symetryczną macierzą trójdiagonalną:

Hψ = αs²εψ, (4.5)

Rozpatrzmy teraz stacjonarne równanie masy efektywnej (3.17). Wprowadźmy na siatce pomocnicze punkty o współrzędnych

x_i+1/2 = a + (i + ¹₂)s, (4.6)

w których określimy wartości masy efektywnej m_i±1/2 = m(x_i±1/2). Stosując dwukrotnie trójpunktowe przybliżenie pierwszej pochodnej otrzymujemy kolejno [177, 178]

 d

Ostatecznie dyskretne jednowymiarowe równanie masy efektywnej zapisujemy jako

− 1

Dla stanów związanych, po uwzględnieniu warunków brzegowych (4.4), otrzymujemy algebraiczne zagadnienie własne (4.5) z symetryczną macierzą trójdiagonalną

H=

Wyraźmy składową x-ową gęstości prądu prawdopodobieństwa j_i+1/2 = j(x_i+1/2), daną wzorem (3.25), poprzez wartości funkcji falowej na siatce. Przyjmijmy

ψ_i+1/2= ψ_i+ ψ_i+1

korzystamy z formuły trójpunktowej dla drugiej pochodnej. W ten sposób dochodzimy do j_i+1/2= 2

Dla kombinacji fal płaskich postaci

ψ(x) = Ae^+ikx+ Be^−ikx (4.10)

ze wzoru analitycznego (3.25) otrzymujemy j = 2k

αmA^∗A − 2k

αmB^∗B = j_→ − j←, (4.11)

gdzie j_→ — prąd cząstek poruszających się w prawo (w kierunku dodatnim), j_← — prąd cząstek poruszających się w lewo. Dla wartości funkcji falowej (4.10) na siatce ze wzoru (3.25) wynika, że

j_i+1/2= 2

αm_i+1/2s=[(Ae^+ik(i+1)s+ Be^−ik(i+1)s)(Ae^−ikis+ Be^+ikis)]

= 2 sin ks

αm_i+1/2s(A^∗A − B^∗B).

W granicy s → 0 wynik ten zgadza się ze wzorem (4.11).

Potrzebny jest jeszcze sposób reprezentacji na siatce potencjałów osobliwych typu V (x) = Ω
(x − x0). Podstawowa równość wynikająca z definicji delty Diraca

będzie spełniona, jeśli wybierzemy V_i = Ω
_i,l/s, gdzie x_l jest punktem siatki najbliż-szym x₀. Wstawiając tak określone wartości potencjału V_i do równania (4.8) i mnożąc je przez s otrzymujemy przy s → 0 skok pochodnej funkcji falowej zgodny z (3.19).

4.1.2. Metoda różnic skończonych dla stanów rozproszeniowych

Jedną z metod rozwiązywania dyskretnego stacjonarnego jednowymiarowego równania Schr¨odingera, dogodną dla stanów niezwiązanych (rozproszeniowych, zdelokalizowanych, tj. opisanych nienormowalną funkcją falową: R+∞

−∞ |ψ(x)|²dx = ∞), jest opisany poniżej wariant metody macierzy przejścia.

Wyrażając w równaniu (4.3) ψ_i+1 jako funkcję ψ_i−1 i ψ_i, możemy napisać macierzowe równanie [175]

gdzie P⁽ⁱ⁾ jest unimodularną macierzą przejścia dla i-tego punktu siatki; łatwo zauważyć, że det P⁽ⁱ⁾ = 1. Warto zwrócić uwagę, że wyrazy tej macierzy są liczbami rzeczywistymi, co dla dowolnych potencjałów czyni metodę siatkową dogodniejszą do obliczeń numerycznych od metody elementów skończonych.

Dyskretne stacjonarne równanie masy efektywnej (4.8) również możemy wyrazić za pomocą macierzy przejścia:

Macierz przejścia P dla całego przedziału ha, bi, zajmowanego przez badany układ fizyczny, jest iloczynem macierzy P⁽ⁱ⁾:

 a jej wyznacznik jest równy stosunkowi mas efektywnych po prawej i lewej stronie rozpa-trywanego przedziału:

det P = P₁₁P₂₂− P12P₂₁ = m_{N +1/2} m_1/2 .

Przyjmując funkcję falową cząstki na zewnątrz przedziału ha, bi w postaci kombinacji fal płaskich: padającej i odbitej po lewej stronie oraz przechodzącej po prawej, otrzymujemy

ψ(x) =

gdzie ρ i τ — odpowiednio fazowe współczynniki odbicia i przejścia, a k i k⁰ — wektory falowe po lewej i prawej stronie.

Po wstawieniu powyższych wartości ψ₋₁, ψ₀, ψ_N i ψ_{N +1} do równania (4.12) powstaje układ równań z niewiadomymi ρ i τ:

τ e^+ik0(b+s) = e^+ika[P₁₁(1 + ρe^−2ika) + P₁₂(e^−iks+ ρe^−2ikae^+iks)]

Mnożąc drugie równanie przez e^+ik0s i porównując z pierwszym, eliminujemy τ. Otrzymu-jemy

ρ = e^2ika−P11− P12e^−iks+ P₂₁e^+ik0s+ P₂₂e^i(k0−k)s

P₁₁+ P₁₂e^+iks− P21e^+ik0s− P22e^i(k0+k)s = e^2ikaZ_n Z_d.

Podstawiając powyższy wynik do drugiego równania w (4.13), dochodzimy do Z_dτ e^i(k0b−ka)= P₂₁(Z_d+ Z_n) + P₂₂(e^−iksZ_d+ e^+iksZ_n)

Gęstości prądów prawdopodobieństwa dla cząstek padających, odbitych i przechodzą-cych wynoszą odpowiednio

jⁱ = sin ks

m , j^r = sin ks

m ρ^∗ρ, j^t= sin k⁰s m⁰ τ^∗τ.

Energetyczne współczynniki odbicia R i przejścia T , zdefiniowane wzorami R = j^r

jⁱ, T = j^t jⁱ,

wyrażają się poprzez elementy macierzy P jako R = ρ^∗ρ = Z_n^∗Z_n

Prawidłowość powyższych formuł można zbadać, obliczając

Wartości energetycznych współczynników odbicia i przejścia nie zależą od kierunku padania cząstek. Jeżeli rozważymy cząstki padające z prawej strony, mamy

ψ(x) = i dalej, po podstawieniu do zależności (4.12),

ρ⁰e^2ik0be^+ik0s = (P₁₁+ P₁₂e^+iks)e^i(k0b−ka)τ⁰− e^−ik0s) ρ⁰e^2ik0b = (P₂₁+ P₂₂e^+iks)e^i(k0b−ka)τ⁰− 1.



 (4.14)

Mnożąc drugie równanie przez e^+ik0s i porównując z pierwszym, eliminujemy ρ⁰. Otrzymu-jemy

a podstawiając ten wynik do drugiego równania w (4.14), dochodzimy do

ρ⁰= e^−2ik0b−P11 − P12e^+iks+ P₂₁e^−ik0s+ P₂₂e^i(k−k0)s

P₁₁+ P₁₂e^+iks− P21e^+ik0s− P22e^i(k0+k)s = e^−2ik0bZ_n^∗ Z_d.

Wynikają stąd równości dla energetycznych współczynników odbicia i załamania:

R⁰ = ρ^0∗ρ⁰ = Z_n^∗Z_n

4.2. Metoda elementów skończonych

W niniejszym podrozdziale przedstawiona jest inna metoda, znajdująca zastosowanie do rozwiązywania równania Schr¨odingera w układach kwazijednowymiarowych, określana zazwyczaj mianem metody macierzy przejścia [165].

4.2.1. Macierz przejścia dla pojedynczej bariery potencjału

Niech jednowymiarowy potencjał V (x) będzie funkcją przedziałami stałą, tak że V (x) = V_i dla x_i < x < x_i+1: umieśćmy dodatkowo w punktach nieciągłości x_i potencjały osobliwe o amplitudach Ω_i proporcjonalne do delty Diraca (rys. 4.1):

V (x) =

Rysunek 4.1. Potencjał przedziałami stały z osobliwościami Ω_i
(x−xi) w punktach nieciągłości x_i

Funkcję falową cząstki poruszającej się w potencjale (4.16) można przedziałami przed-stawić jako kombinację fal płaskich

ψ(x) =

Na granicach przedziałów, tj. w punktach x_i, spełnione są warunki ciągłości (3.19) funkcji falowej i jej pochodnej. Dla funkcji (4.17) możemy je zapisać w zsymetryzowanej postaci

Rozwiązywanie stacjonarnego równania Schr¨odingera (4.2) lub masy efektywnej (3.17) na przedziałach hxi, x_i+1i i zszywanie rozwiązań przy pomocy warunków (4.18) należy do klasy metod zwanych metodami elementów skończonych. Opisany poniżej wariant nazy-wany jest metodą macierzy przejścia; przedstawiamy go dla najbardziej ogólnego przy-padku, tj. równania masy efektywnej z potencjałem osobliwym (4.16). Autor nie spotkał się dotychczas z takim podejściem w literaturze, a potencjał tego typu pojawia się w ma-gnetycznym modelu Kroniga–Penneya [135] po uwzględnieniu spinu elektronu.

Dla funkcji falowej postaci (4.17) warunki (4.18) można zapisać w formalizmie macie-rzowym:

przy czym Wyznacznik tej macierzy ma wartość

det M(k_i−1, m_i−1, k_i, m_i, x_i) = k_i−1

Rysunek 4.2. Pojedyncza prostokątna bariera potencjału Dla prostokątnej bariery potencjału (rys. 4.2)

V (x) =

możemy przyjąć funcję falową postaci

ψ_x =

×

Powyższe formuły pozostają słuszne dla cząstek o energii niższej niż energia bariery (ε < V_b), tj. kiedy k⁰ = i|k⁰|. Można wówczas wyeliminować urojone argumenty funk-cji trygonometrycznych, wprowadzając C⁰ = cosh[|k⁰|(b − a)], S⁰ = sinh[|k⁰|(b − a)], i w przypadku m_i = m_i−1 otrzymujemy

M_↑(x₀) =

Ten sam wynik można otrzymać, przechodząc do zera z szerokością prostokątnej bariery (a = x_i− ¹₂, b = x_i+ ¹₂) przy zachowaniu stałego iloczynu (b − a)Vb = Ω.

Dodając do prostokątnej bariery (4.23) potencjał osobliwy postaci V (x) = Ω[
(x−a)−

(x − b)], otrzymujemy macierz przejścia

M_↑u↓(k, m, k⁰, m⁰, Ω, a, b) = (4.27) taka szczególna postać potencjału osobliwego pojawia się w magnetycznym modelu Kro-niga–Penneya ze spinem. Dla stanów podbarierowych (ε < V_b) można ponownie wyelimi-nować urojone argumenty funkcji trygonometrycznych:

M_↑u↓(k, m, |k⁰|, m⁰, Ω, a, b) = (4.28) Nietrudno pokazać, że wyznaczniki macierzy (4.24)–(4.28) są równe jedności (są to macierze unimodularne), co wynika ze wzoru (4.22):

det M_u(k, m, k⁰, m⁰, Ω, a, b) =

Ponadto, o ile k jest liczbą rzeczywistą (tj. energia cząstek ε ­ 0), macierz ta ma szczególną postać

niezależnie od tego, czy k⁰ jest liczbą rzeczywistą, czy urojoną (tj. czy ε > V_b — cząstka przechodzi nad barierą, czy też ε < V_b — cząstka tuneluje).

4.2.2. Dowolny układ prostokątnych barier potencjału

Macierz przejścia M dla układu n prostokątnych barier (rys. 4.1) jest iloczynem po-szczególnych macierzy:

M = M^{(N )}_u M^{(N −1)}_u · · · M⁽²⁾u M⁽¹⁾_u , (4.30)

a zatem det M = 1. Warto też zauważyć, że struktura (4.29) macierzy zachowuje się przy ich mnożeniu: więc również macierz M jest postaci (4.29). Dodatkowo, ze względu na unimodularność, mamy

W sytuacji, gdy potencjał składa się z pewnej liczby jednakowych barier (lub barier kilku rodzajów) położonych w punktach x_i, można wykorzystać fakt, że przesunięcie ba-riery o ∆x powoduje tylko prostą transformację opisującej ją macierzy przejścia:

M_u(k, k⁰, a + ∆x, b + ∆x) = E(−k∆x)Mu(k, k⁰, a, b)E(k∆x).

Znając zatem macierze przejścia M⁽ⁱ⁾_u0 dla barier położonych w x = 0, możemy przepisać wzór (4.30) w postaci

M = E(−kxN)M⁽ⁿ⁾_u0E(kx_N)E(−kxN −1)M^{(N −1)}_u0 E(kx_{N −1})

× · · · E(−kx1)M⁽¹⁾_u0E(kx₁),

a podstawiając d_i = x_i+1− xi — długość segmentu zawierającego i-tą barierę, M = E(−kL)h z macierzy E(kd_i)M⁽ⁱ⁾_u0 zachowuje strukturę (4.29).

Macierz przejścia M dla dowolnego potencjału postaci (4.16) jest iloczynem macierzy postaci (4.24), o ile k₀ = k_N (V₀ = V_N). Dowód wystarczy przeprowadzić dla n = 3.

ponieważ E(a)E(b) = E(a + b), co kończy dowód.

Powyższy fakt umożliwia zastosowania metody macierzy przejścia do dowolnych po-tencjałów, które dają się przybliżać funkcjami postaci (4.16), przy zmniejszaniu odległości x_i+1 − xi. Odpowiedni algorytm numeryczny wymaga mnożenia dużej liczby macierzy zespolonych rozmiaru 2 × 2.

Przyjmując, że układ barier zajmuje przedział ha, bi, oraz że po prawej i po lewej jego stronie funkcja falowa cząstki jest kombinacją fal płaskich, dla cząstki padającej z lewej strony otrzymujemy

ψ(x) =

(e^+ikx+ ρe^−ikx (x ¬ a) τ e^+ik0x (x ­ b)

gdzie ρ i τ — odpowiednio fazowe współczynniki odbicia i przejścia (liczby zespolone), a k i k⁰ — wektory falowe na lewo i na prawo od układu barier.

Wybrane w ten sposób funkcje falowe związane są poprzez macierz przejścia układu barier:

otrzymujemy fazowe współczynniki odbicia i przejścia ρ = −M_1,2^∗

M_1,1^∗ , τ = M_1,1− M_1,2^∗ M_1,2

M_1,1^∗ = det M

M_1,1^∗ = 1 M_1,1^∗ .

Przy założeniu k⁰ = k, energetyczne współczynniki odbicia i przejścia oraz ich iloraz przyjmują wartości

Wynika stąd, że do obliczenia energetycznych współczynników odbicia i przejścia wystar-cza znajomość jednego z diagonalnych wyrazów macierzy przejścia. Ponadto pomnożenie macierzy M przez macierz E(u) postaci (4.20) zmienia tylko fazę ρ i τ, a nie zmienia w ogóle wartości R i T .

Dla cząstek padających na układ barier z prawej strony otrzymujemy

ψ(x) =

łatwo otrzymujemy (przy założeniu k = k⁰) ρ⁰= M_1,2

M_1,1^∗ = −ρ, τ⁰ = 1

M_1,1^∗ = τ, R⁰ = R, T⁰ = T .

4.2.3. Równoważność metody elementów skończonych i metody różnic skończonych Z układu równań (4.33) możemy znaleźć wyrazy macierzy przejścia M opisującej w me-todzie elementów skończonych ten sam układ barier potencjału, co macierz P w meme-todzie siatkowej. Otrzymujemy M_1,1 = 1/τ^∗, M_1,2= −ρ^∗/τ^∗ i stąd

Podobnie dla strumienia cząstek padających na układ z prawej strony: wyznaczona z układu równań (4.34) macierz przejścia w metodzie elementów skończonych

M =

jest równa macierzy (4.35), zgodnie z równościami (4.15).

Uzyskanie macierzy przejścia dla układu barier potencjału wymaga mnożenia macierzy zespolonych. Jeżeli potencjał V (x) i inne właściwości układu, jak np. m(x), zadane są funkcjami ciągłymi, wówczas ich dokładna aproksymacja funkcjami schodkowymi (4.16) wymaga użycia dużej liczby krótkich przedziałów hxi, x_i+1i i liczba tych macierzy jest duża.

Z kolei metoda siatkowa wymaga mnożenia dużej liczby macierzy (gdyż krok siatki s musi być dostatecznie mały) niezależnie od postaci potencjału, lecz są to macierze rzeczywiste o prostej postaci i operacje na nich są znacznie mniej kosztowne numerycznie. Można zatem stwierdzić, że metoda siatkowa jest bardziej uniwersalna, natomiast metoda elementów skończonych jest szczególnie efektywna dla układów, których lokalne charakterystyki są przedziałami stałe.

4.3. Formalizm map śladów i antyśladów

Poniżej przedstawiony jest formalizm umożliwiający efektywne obliczanie współczynni-ków przejścia i odbicia dla wybranych układów aperiodycznych — uogólnionych supersieci Fibonacciego. Wykorzystuje on szczególne właściwości macierzy przejścia występującej w metodzie elementów skończonych.

4.3.1. Właściwości wielomianów Czebyszewa i macierzy unimodularnych

Niech W będzie dowolną macierzą unimodularną rozmiaru 2×2, tj. det W = W11W₂₂− W₁₂W₂₁ = 1, i niech U_j(t) oznacza wielomiany Czebyszewa drugiego rodzaju, spełniające rekurencyjne związki [180]

U_j(t) ≡ 0 (j < 0), U₀(t) ≡ 1, U_j(t) = 2tU_j−1(t) − Uj−2(t) (j > 0) (4.36) oraz

1 + U_j+1(t)U_j−1(t) = U_j²(t). (4.37)

Drugiego związku dowodzi się indukcyjnie. Korzystając z równości 1 + U_j(t)U_j−2(t) = U_j−1² (t), 1 + U_j−1(t)U_j−3(t) = U_j−2² (t), słusznych — jak łatwo sprawdzić — dla j ¬ 2, otrzymujemy kolejno 1 + U_j+1(t)U_j−1(t) =

= 1 + [2zU_j(t) − Uj−1(t)][2zU_j−2(t) − Uj−3(t)]

= 1 + U_j−1(t)U_j−3(t) + 4z²U_j(t)U_j−2(t) − 2z[Uj(t)U_j−3(t) + U_j−1(t)U_j−2(t)]

= U_j−2² (t) + 4z²[U_j−1² (t) − 1]

− 2z{[2zUj−1(t) − Uj−2(t)]U_j−3(t) + [2zU_j−2(t) − Uj−3(t)]U_j−2(t)}

= U_j−2² (t) + [2zU_j−1(t)]²− 4z{z[1 + Uj−1(t)U_j−3(t) + U_j−2² (t)] − Uj−2(t)U_j−3(t)}

= U_j−2² (t) + [2zU_j−1(t)]²− 4z{Uj−2(t)[2zU_j−2(t) − Uj−3(t)]}

= [U_j−2(t) − 2zUj−1(t)]² = U_j²(t).

Potęgi dowolnej macierzy unimodularnej W rozmiaru 2 × 2 wyrażają się wzorem [181]

W^j = U_j−1(¹₂Tr W)W − Uj−2(¹₂Tr W)1. (4.38) Powyższego wzór, którego prawdziwość dla j = 1 jest oczywista, dowodzi się indukcyjnie.

Niech t = ¹₂Tr W. Mamy kolejno

W^j+1 = U_j(t)W − Uj−1(t)1 = [2tU_j−1(t) − Uj−2(t)]W − Uj−1(t)1

= U_j−1(t)[2tW − I] − Uj−2(t)W = U_j−1(t)W²− Uj−2(t)W = W^jW.

Dla dowolnej nieosobliwej macierzy P rozmiaru 2 ×2 możemy zastosować podstawienie W= (det P)^−1/2P. Wtedy det W = 1, Tr W = (det P)^−1/2Tr P oraz, ze wzoru (4.38),

P^j = [√

det P ]^j−1U_j−1(¹₂Tr P/√

det P )P − [√

det P ]^jU_j−2(¹₂ Tr P/√

det P )1.

4.3.2. Odwzorowania śladów dla uogólnionych supersieci typu Fibonacciego

Zgodnie z zasadą konstrukcji z elementów A i B uogólnionego łańcucha Fibonacciego S_l rzędu l, tj.

S₀ = B, S₁ = A, S_l = S_l−1ⁿ  Sl−2^m ,

dla supersieci zbudowanej z segmentów opisanych macierzami przejścia A= E(kd_A)M_A, B= E(kd_B)M_B,

gdzie M_i (i = A, B) jest jedną z macierzy (4.24)–(4.28) w metodzie elementów skończonych lub dowolną macierzą postaci (4.35) wyznaczoną przy użyciu metody siatkowej.

Zgodnie ze wzorem (4.32) mamy

S₀ = B, S₁ = A, S_l = S^m_l−2Sⁿ_l−1.

Macierz przejścia dla całej supersieci rzędu l o długości L ma postać M = E(−kL)Sl.

Każda z macierzy S_l zachowuje strukturę (4.29): i jest unimodularna, zatem ze wzoru (4.38) mamy

S_l = [U_m−1(t_l−2)S_l−2− Um−2(t_l−2)1][U_n−1(t_l−1)S_l−1− Un−2(t_l−1)1] Ze wzorów (4.39) i (4.41) otrzymujemy zależności

s_l = a_l(s_l−2s_l−1+ z_l−2z_l−1^∗ ) − bls_l−2− cls_l−1+ d_l, z_l = a_l(s_l−2z_l−1+ z_l−2s^∗_l−1) − blz_l−2− clz_l−1.

)

(4.42) Z pierwszego związku można wyeliminować wyrazy pozadiagonalne z_j. Podstawiając l 7→ l − 1, mamy

s^∗_l−1= a_l−1(s^∗_l−3s^∗_l−2+ z_l−3^∗ z_l−2) − bl−1s^∗_l−3− cl−1s^∗_l−2+ d_l−1, z_l−1^∗ = a_l−1(s^∗_l−3z_l−2^∗ + z_l−3^∗ s_l−2) − bl−1z_l−3^∗ − cl−1z_l−2^∗ .

)

(4.43) Wstawiając drugie z powyższych równań do pierwszej z formuł (4.42), dostajemy

s_l = a_l{sl−2s_l−1+ z_l−2[a_l−1(s^∗_l−3z_l−2^∗ + z^∗_l−3s_l−2) − bl−1z^∗_l−3− cl−1z_l−2^∗ ]} Zajmijmy się wyrażeniem zapisanym powyżej w nawiasach kwadratowych, stosując formuły (4.36) i (4.37):

(2a_l−1b_l−1t_l−2− a²l−1− b²l−1)s^∗_l−3− 2al−1d_l−1t_l−2+ b_l−1d_l−1+ a_l−1c_l−1− bl−1s^∗_l−1=

= U_m−1² (t_l−3)[2U_n−1(t_l−2)U_n−2(t_l−2)t_l−2− Un−1² (t_l−2) − Un−2² (t_l−2)]s^∗_l−3

Wstawiając ten wynik do równania (4.44), otrzymujemy ostatecznie dla l > 3 s_l = U_m−1(t_l−2)U_n−1(t_l−1) Powyższe zespolone odwzorowanie jest nowym wynikiem, pozwalającym na znacznie prostsze wyprowadzenie odwzorowań rzeczywistych użytych w pracach [78, 79].

Jak widać ze wzoru (4.45), do wyznaczenia wyrazów diagonalnych macierzy S_l (a co za tym idzie, współczynników przejścia T = 1/|sl|² i odbicia R = 1 − T dla uogólnionej supersieci typu Fibonacciego l-tego rzędu) wystarcza znajomość wyrazów diagonalnych trzech macierzy niższych rzędów: S_l−1, S_l−2, S_l−3.

Wzór (4.45) ma sens matematyczny pod warunkiem, że U_n−1(t_l−2) 6= 0, skąd wynika wcześniejsze założenie a_l−1 6= 0. Jeżeli Un−1(t_l−2) = 0, we wzorze tym musi występować nieoznaczoność typu 0/0, gdyż mnożenie macierzy o skończonych wyrazach nie może dać s_l = ∞. Nieoznaczoność wystąpi, jeśli tl−2jest miejscem zerowym wielomianu U_n−1(t). Nie stanowi to większego problemu w obliczeniach numerycznych, ponieważ trafienie dokładnie w miejsce zerowe jest mało prawdopodobne, a jeśli się zdarzy, można powtórzyć obliczenia, zmieniając nieznacznie wartość energii cząstki (i zarazem wyrazy macierzy przejścia).

Nietrudno zauważyć, że wyznaczając t_l = <(sl) = ¹₂(s_l+ s^∗_l) = ¹₂ Tr S_j ze wzoru (4.45), otrzymujemy formułę zawierającą tylko połowy śladów poszczególnych macierzy:

t_l = U_m−1(t_l−2) Formułę tę można określić terminem odwzorowanie śladów.

Podobnie wyznaczając połowy antyśladów u_l = =(sl) = _2i¹(s_l − s^∗l) otrzymujemy for-mułę

u_l = U_m−1(t_l−2)

"

U_n(t_l−1)u_l−2+ U_n−1(t_l−1)

U_n−1(t_l−2)[U_n−2(t_l−2)u_l−1+ U_m−1(t_l−3)u_l−3]

#

− Um−2(t_l−2)U_n−1(t_l−1)u_l−1. (4.47)

którą nazwać można odwzorowaniem śladów–antyśladów.

Używając liczb rzeczywistych t_l i u_l, zapisujemy współczynniki przejścia i odbicia l-tego pokolenia supersieci jako

T = 1

t²_l + u²_l , R = 1 − T .

Można pokazać, że liczba elementarnych operacji arytmetycznych potrzebnych do wy-znaczenia T przy użyciu odwzorowań śladów i antyśladów jest mniejsza niż przy bezpo-średnim mnożeniu potęg macierzy według formuły (4.40), nawet przy wykorzystaniu szcze-gólnej struktury macierzy S_l. Porównanie nakładu obliczeniowego obu metod przedstawia tabela 4.1. Przy szacowaniu efektywności należy wziąć pod uwagę fakt, że na typowych współczesnych komputerach czas wykonania operacji mnożenia i dodawania jest zbliżony, zaś dzielenie jest kilkakrotnie wolniejsze [182, 183].

Tabela 4.1. Liczba elementarnych operacji arytmetycznych przy wyznaczaniu współ-czynnika przejścia dla kolejnego pokolenia uogólnionej supersieci Fibonacciego dla me-tody bezpośredniego mnożenia macierzy i meme-tody odwzorowań śladów–antyśladów. Po-dane liczby zależą od szczegółów implementacji algorytmu numerycznego i mogą się nieznacznie zmienić przy jego modyfikacji

Operacja arytmetyczna Mnożenie macierzy Metoda odwzorowań

dodawanie/odejmowanie 14 8

mnożenie 24 14

dzielenie 0 1

Z drugiej strony, stosując odzwzorowania śladów–antyśladów w prezentowanej postaci, tracimy informację o wyrazach pozadiagonalnych macierzy przejścia, a co za tym idzie, o fa-zie współczynnika odbicia ρ. Możliwe jest co prawda zapisanie podobnych odwzorowań dla wyrazów pozadiagonalnych z_l (metoda ta stosowana jest np. w obliczeniach dotyczących transmisji światła przez aperiodczne układy wielowarstwowe [90]), jednak wykorzystanie ich wiąże się z dodatkowymi operacjami arytmetycznymi i przewaga metody odwzorowań znika. Zatem metoda odwzorowań śladów–antyśladów jest efektywna, ale nie uniwersalna.

4.4. Wyznaczanie rozkładu potencjału elektrostatycznego w układzie

Znalezienie rozkładu potencjału elektrostatycznego φ(r) w niejednorodnym przestrzen-nie układzie wymaga rozwiązania równania Poissona [184]

−∇^ (r)∇φ(r) = %(r),

w którym^ (r) — zależny od położenia tensor przenikalności elektrycznej, a %(r) — gęstość ładunku. W przypadku układu kwazijednowymiarowego zajmującego obszar a ¬ x ¬ b, do którego przyłożono wzdłuż osi X zewnętrzne napięcie U, równanie upraszcza się do postaci

− d

dxκ(x) d

dxφ(x) = %(x) (4.48)

z warunkami brzegowymi (z dokładnością do addytywnej stałej cechowania) φ(a) = 0, φ(b) = U.

Operator różniczkowy w równaniu (4.48) ma postać analogiczną do operatora energii ki-netycznej w równaniu masy efektywnej (3.17) i poddaje się podobnemu do (4.7) schematowi dyskretyzacji. Na siatce punktów (4.1) x₀, . . . , x_N określone są wartości gęstości ładunków

%_i = %(x_i) oraz szukane wartości potencjału φ_i = φ(x_i), a na siatce pomocniczej (4.6)

prowadzi do dyskretnej postaci jednowymiarowego równania Poissona:

−κi−1/2φ_i−1+ (κ_i−1/2+ κ_i+1/2)φ_i− κi+1/2φ_i+1 = ˜%_iφ_i (i = 1, . . . , N), (4.49) gdzie ˜%_i = s²%_i oznacza przeskalowaną gęstość ładunku. Przenosząc wyrazy b₀ = κ_1/2φ₀ oraz b_N = κ_{N +1/2}φ_N (określające warunki brzegowe) na prawą stronę, otrzymujemy układ równań o macierzowym zapisie

Jest to układ równań z dodatnio określoną symetryczną macierzą trójdiagonalną, który można efektywnie rozwiązywać kosztem O(n) operacji arytmetycznych za pomocą standar-dowych algorytmów numerycznej algebry liniowej, np. z wykorzystaniem rozkładu LDL^T macierzy [185, 186].

5

Wyniki numeryczne

Niniejszy rozdział przedstawia wyniki obliczeń przewodnictwa elektrycznego, otrzy-mane w ramach opisanego modelu kwantowego dla wybranych aperiodycznych super-sieci półprzewodnikowych i supersuper-sieci magnetycznych, których uporządkowanie odpowiada uogólnionym supersieciom Fibonacciego. Rezultaty obliczeń wykonanych w celach porów-nawczych dla wybranych supersieci periodycznych zawarte są w dodatku A.

5.1. Supersieci półprzewodnikowe

Rozważany w pracy model supersieci półprzewodnikowej został opisany w podroz-dziale 2.2. Przyjęte dla większości badanych układów parametry modelu odpowiadają su-persieciom GaAs-Al_xGa_1−xAs dla ułamka molowego x ≈ 0,3, których technologia wytwa-rzania jest dobrze opanowana. Rozmiary liniowe rozpatrywanych supersieci są rzędu kilku-nastu nanometrów. Dla wybranych aperiodycznych i periodycznych supersieci półprzewod-nikowych wyznaczono przewodność elektryczną oraz przewodność różniczkową w funkcji parametrów modelu oraz zewnętrznego napięcia U i poziomu Fermiego ε_F w kontaktach.

5.1.1. Schemat obliczeń numerycznych

Obliczenia numeryczne przewodnictwa supersieci półprzewodnikowych przeprowadzone zostały według przedstawionego poniżej schematu.

Danymi wejściowymi były:

typ uporządkowania układu, określający rozmieszczenie komórek typu A i B — dla uogólnionych supersieci Fibonacciego jest on zadany parametrami l, m i n w formu-łach (2.3) i (2.5);

szczegółowe postaci komórek A i B, zadane przez długości komórek L_i oraz wartości potencjału V_i, masy efektywnej m_i i przenikalności elektrycznej κ_i;

wartość zewnętrznego napięcia U przyłożonego do układu;

położenie poziomu Fermiego ε_F w kontaktach.

Kolejne kroki algorytmu numerycznego polegały na:

1. Ustaleniu porządku komórek elementarnych w supersieci — dla uogólnionych supersieci Fibonacciego stosowana była reguła konkatenacji (2.3);

55

2. Wyznaczeniu profilu potencjału, masy efektywnej i przenikalności elektrycznej na całej supersieci na podstawie danych opisujących komórki oraz porządku przestrzennego komórek;

3. Obliczeniu rozkładu potencjału elektrostatycznego φ(x) w układzie przy zadanym na-pięciu U poprzez rozwiązanie równania Poissona (4.50);

4. Obliczeniu wartości przewodności dla określonej energii ε elektronów i kąta ϑ zadają-cego składową poprzeczną k_⊥ wektora falowego: G = G(ε, ϑ); krok ten był wielokrotnie powtarzany dla różnych wartości θ i ε przy numerycznym obliczaniu całek postaci (3.8) i (3.10).

Na podstawie otrzymanych wartości G wyznaczane były dodatkowo:

charakterystyki I = f(U);

przewodność różniczkowa G_diff = dI/dU (poprzez numeryczne różniczkowanie).

Z uwagi na obecność zewnętrznego pola elektrycznego, stosowany był algorytm nume-ryczny wykorzystujący metodę różnic skończonych. Został on przetestowany przez porów-nanie wyników dla U = 0 z wynikami dawanymi przez algorytm wykorzystujący odwzoro-wania śladów i antyśladów.

5.1.2. Wyniki obliczeń

Poniżej, na rysunkach 5.1–5.35, przedstawione są wyniki numerycznych obliczeń prze-wodności G i przeprze-wodności różniczkowej G_diff otrzymane dla rozważanych supersieci pół-przewodnikowych. Są one rozszerzeniem na przypadek przewodnictwa wielokanałowego wcześniej opublikowanych rezultatów [81].

Przewodność elektryczna została wyznaczona dla pojedynczych kanałów przewodnic-twa oraz przy uśrednianiu po powierzchni Fermiego. Zbadany został także wpływ

Przewodność elektryczna została wyznaczona dla pojedynczych kanałów przewodnic-twa oraz przy uśrednianiu po powierzchni Fermiego. Zbadany został także wpływ

W dokumencie PRZEWODNICTWO ELEKTRYCZNE UKŁADÓW NISKOWYMIAROWYCH (Stron 35-121)