################################################################################
Teoria otwartych układów kwantowych
Heinz Peter Breuer, Francesco Petruccione
Tytuł oryginału : „The theory of open quantum systems”
Oxford University Press 2002, 2003
************************************************************************************************
Tłumaczenie : R. Waligóra Pierwsze tłumaczenie : 2014
Ostatnia modyfikacja : 2018-10-20 Tłumaczenie całości książki ( w dwóch częściach )
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
CZĘŚC II Skróty.
M-C – Monte -Carlo
************************************************************************************************
Rozdział 7 Stochastyczne metody modelowania.
Najczęściej spotykane procesy stochastyczne w nauce nie mogą być analizowane z użyciem metod analitycznych.
Dlatego też należy polegać na metodach numerycznych w celu prognozowania w oparciu o zadana dynamikę
stochastyczną. Algorytmy numeryczne mogą być wykorzystywane w celu rozwiązywania równania ruchu, kierującego rozkładem probabilistycznym danej wielkości stochastycznej, równania podstawowego lub np. równania Fokkera- Plancka. Inne podejście polega na wprowadzeniu układu deterministycznych rr dla momentów danego procesu i
rozwiązania ich z użyciem całkowania numerycznego. Jednakże metody takie nie mogą być wielokrotnie zrealizowane w praktyce w wyniku wysokiego stopnia złożoności, leżącego u podstaw danego zagadnienia. Rozkład probabilistyczny może być przedstawiony w sposób pełny w przestrzeni wielowymiarowej (nawet nieskończenie wymiarowej ), podczas gdy dynamika momentów może być sprowadzona do nieograniczonej hierarchii sprzężonych równań nieliniowych. W obu tych przypadkach całkowanie numeryczne z użyciem deterministycznych algorytmów jest skrajnie trudne.
Alternatywnym i bardzo mocnym instrumentem analizy numerycznej procesów stochastycznych są tzw. metody Monte- Carlo (* M-C *) lub stochastyczne metody modelowania. Podstawowa idea takich metod polega na wykorzystaniu algorytmów numerycznego generowania niezależnych realizacji procesu stochastycznego i ocenie wszystkich średnich matematycznych takich realizacji z pomocą uśredniania statystycznego. Modelowanie stochastyczne, jest zatem równoważne realizacji eksperymentu z użyciem komputera. To zaś prowadzi do wyników dla pojedynczych wyników pomiaru z odpowiednim prawdopodobieństwem dając jako dodatek do wartości średnich oceny błędów statystycznych interesujących nas wielkości.
W niniejszym rozdziale rozpatrujemy w pewnych detalach zastosowanie metod Monte-Carlo do procesów
stochastycznych w przestrzeni Hilberta, które to analizowaliśmy w rozdziale 6. Niezależnie od tego, czy dany proces jest procesem kawałkami zdeterminowanym lub tez procesem dyfuzyjnym w przestrzeni Hilberta, odpowiednie oddzielne jego realizacje składają się z interwałów deterministycznych ewolucyjnych okresów, przeplatanych natychmiastowymi skokami lub tez nigdzie nieróżniczkowalnymi trajektoriami. Odpowiednie metody modelowania obu takich typów procesów opisane będą w podrozdziałach 7.1 i 7.2 ich ilustracje na specjalnych przykładach podamy w podrozdziale 7.3
Aby numerycznie określić macierz gęstości dla otwartego układu kwantowego, można bezpośrednio scałkować równanie dla macierzy gęstości lub, inaczej zamodelować proces dla stochastycznej funkcji falowej i ocenić macierz kowariancji.
W ogólności, równanie dla macierzy gęstości prowadzi do układu równań liniowych, zawierających N2 zmiennych zespolonych, gdzie N oznacza efektywny wymiar przestrzeni Hilberta. W przeciwieństwie do tego, modelowanie stochastyczne wymaga jedynie rozpatrzenia N zmiennych zespolonych, charakteryzujących wektor stanu.
Dla dużego N, tj. dla wysokiego wymiaru przestrzeni Hilberta, oczekujemy iż modelowanie z użyciem metod M-C będzie bardziej efektywne, niż całkowanie równania dla macierzy gęstości, ponieważ w tym przypadku zapewniona jest wymagana liczebność reprezentacji, która to nie rośnie zbyt szybko wraz ze wzrostem N.
Taka własność będzie zilustrowana w podrozdziale 7.4, w którym to podajemy szczegóły porównania analizy
7.1 Numeryczne algorytmy modelowania PDP.
W kilku miejscach w rozdziale 6 pokazaliśmy już wyniki, otrzymane przy modelowaniu numerycznym PDP w
przestrzeni Hilberta. Obecnie omówimy szczegółowo metody modelowania i szczególne własności odpowiednich metod numerycznych. Podstawowa część zastosowań metod M-C w fizyce statystycznej i fizyce materii skondensowanej może być znaleziona w zbirze prac pod redakcją Binder’a (Binder 1995 ) oraz w książce Landaua i Binder’a
( Landau, Binder 2000 )
7.1.1 Ocena średnich matematycznych.
Stochastyczny algorytm modelowania służy nam dla generacji próbki niezależnych realizacji danego procesu
stochastycznego ψ(t) dla funkcji falowej. Oznaczmy takie realizacje jako ψr(t), r = 1,2, ..., R gdzie R jest liczbą realizacji w próbce. Cel modelowania polega na ocenie średnich matematycznych :
Mt ≡ E( F[ψ, t] ) = ∫ P[ψ, t]F[ψ, t]DψDψ* (7.1) lub funkcjonałów rzeczywistych F[ψ, t] losowego wektora stanu, który może zależeć jawnie od czasu. Obiektywna i konsekwentna ocena dla średniej matematycznej Mt zadana jest przez średnią po próbkach :
Teraz i dalej daszek wykorzystywany jest w celu oznaczenia oszacowania średniego. Szczególny interes reprezentują dla nas średnie matematyczne postaci :
obserwabli B układu otwartego, dla której ocena przyjmuje postać :
Jest jasne, że ocena średniej matematycznej podlega określonym błędom statystycznym.
Naturalną miarą dla fluktuacji statystycznych wielkości losowej F[ψ, t] jest dyspersja :
W przypadku F = < ψ | B | ψ > (7.5) zadaje nic innego jak dyspersje :
σt2 = Var2(B) (7.6)
która już jest zawarta w równaniach (5.14) i (5.68). Odpowiednia ocena dla błędów statystycznych w definicji Mt z ograniczonej próbki o objętości R zadana jest przez wyrażenie :
Wielkość σt^ znana jest jako błąd standardowy wartości średniej Mt. Dla oceny wartości średniej obserwabli B zadana jest ona przez wyrażenie :
Jeśli realizacje w próbce są statystycznie niezależne, błąd standardowy σt^ ma postać pierwiastka kwadratowego z objętości próbki R :
σt^ ~ 1/√R (7.9)
Jest to zależność pomiędzy błędem standardowym wartości średniej i liczba realizacji próbki.
7.1.2 Generowanie realizacji procesu.
Rozpatrzmy QND w przestrzeni Hilberta o postaci analizowanej w podrozdziale 6.1.1, który może być przedstawiony poprzez srr (zobacz (6.25)) :
Inkrementy poissonowskie dNi(t) spełniają równania (6.26) i (6.27), jednocześnie niehermitowski hamiltonian H6 zadany jest przez wyrażenie (zobacz (6.13)) :
H^ = H – ½ i Σ γi A†
i Ai (7.11)
i
Z podrozdziału 6.1.1 wiadomo również, że rozkład kumulatywny czasu oczekiwania (6.21) dla skoków kwantowych zadany jest przez wyrażenie :
F[ψ~, t] = 1 – || exp( –iH^τ) ψ~ ||2 (7.12)
gdzie τ reprezentuje losowy czas oczekiwania pomiędzy dwoma skokami.
Z analiz przeprowadzonych w podrozdziałach 1.5.2 i 6.1.1 widzimy, że przykład realizacji ψr(t) procesu, określony przez równanie (7.10) w interwale czasu [0, tf ], może być wygenerowany przez następujący algorytm :
1) Przyjmujemy, że unormowany stan ψr(t) został osiągnięty poprzez skok w chwili t i przyjmujemy, że ψr(t) = ψ~.
Jeśli t jest czasem początkowym (t = 0 ), to ψ~ jest stanem początkowym procesu, który powinien być wybrany z rozkładu początkowego P[ψ~,t = 0 ].
2) Definiujemy losowy czas oczekiwania τ zgodnie z funkcją rozkładu (7.12). Można to zrobić np. poprzez losowanie liczby η, która jest rozłożona równomiernie w interwale [0, 1] i określone poprzez τ z równania :
η = 1 – F[ψ~, τ ] = || exp( –iH^τ) ψ~ ||2 (7.13) Dla η > q istnieje jednoznaczne rozwiązanie; q – jest defektem rozkładu czasu oczekiwania, określonym w równaniu (6.22). Dla η ≤ q ustanawiamy τ = ∞, to oznacza, że dalszych skoków nie będzie.
W granicach interwału czasu [t, t + τ ] realizacja zgodnie z ewolucją czasową, zadawaną przez (6.15) :
3) W chwili t + τ (jeśli τ jest określony i t + τ < tf ) następuje jeden z możliwych skoków, oznaczanych przez indeks i w równaniu (7.10). Wybieramy określony typ skoku i z prawdopodobieństwem :
i zamieniamy :
4) Powtarzamy kroki 1 – 3 do chwili póki, nie osiągniemy żądanego czasu tf , który prowadzi do realizacji ψr(t) na całym interwale czasu [0, tf ].
5) Wynikiem będzie wygenerowana próbka realizacji ψr(t) ; r = 1, 2, ... ,R zgodnie z powyższym algorytmem. Dowolna wielkość statystyczna może być oceniona poprzez odpowiednie średnie ansamblu, tak jak to opisano w poprzednim podrozdziale.
Dalej przeanalizujemy oddzielne części algorytmu dokładniej.
7.1.3 Określenie czasu oczekiwania.
Specyficzna częścią algorytmu modelowania jest określenie losowego czasu oczekiwania.
7.1.3.1 Eksponencjalne rozkłady czasów oczekiwania.
Rozpatrzmy teraz najprostszy przypadek, ekspotencjalnego czasu oczekiwania :
P[ψ~,τ ] = 1 – exp(–Γτ ) (7.17)
Taki rozkład czasu oczekiwania pojawia się zawsze, kiedy sumaryczna intensywność skoku nie zależy od czasu tj.
Γ[ ψ~,τ ] = Γ = const. > 0
Jeśli η jest równomiernie rozłożoną wielkością losową na interwale [0, 1], to czas oczekiwania τ łatwo określić z pomocą równania (7.13), które daje :
τ = – (1/Γ ) ln(η ) (7.18)
Realizacja odpowiednich generatorów liczb losowych omawiana jest np. w Press, Flannery, Teukolsky, Vetterling (1992) oraz Knuth (1981)
7.1.3.2 Multieksponencjalne rozkłady czasów oczekiwania.
Rozkład czasów oczekiwania, który jest suma funkcji ekspotencjalnych, pojawia się jeśli operatory skoku Ai są operatorami własnymi hamiltonianu H, co prowadzi do (zobacz (3.125)) :
Taki przypadek miał miejsce np. w słabo powiązanym równaniu podstawowym (3.140) bez zewnętrznych pól wzbudzenia. Stąd wynika, że H oraz dodatni operator :
Σ γi A† i Ai
mają ogólna wspólna bazę { | α > } tj. :
z rzeczywistymi wartościami własnymi Eα i Γα ≥ 0. Jednakże H^ w takiej bazie jest diagonalny i możemy zapisać :
H^ | α > = ( Eα – ½ Γα ) | α > (7.22)
Stąd wynika, że rozkład czasu oczekiwania jest suma funkcji ekspotencjalnych, których eksponenty zadane są przez urojone wartości własnych H^. Zatem, czas oczekiwania τ jest określony przez równanie :
η = || exp( –iH^τ) ψ~ ||2 = Σ | < α | ψ~ > |2 exp( –Γατ ) (7.23) Jeśli ψ~
jest proporcjonalna do jednego z | α >, to otrzymujemy ekspotencjalny czas oczekiwania (7.17). W ogólnym przypadku równanie (7.23) nie może być rozwiązane analitycznie względem τ. Jednakże może być ono rozwiązane numerycznie przy warunku, że | α > i odpowiednie wartości własne Γα są znane. Równanie (7.23) pokazuje również, że mamy defekt jeśli Γα = 0 dla pewnego α z < α | ψ~
> ≠ 0.
7.1.3.3 Podstawowe rozkłady dla czasów oczekiwania.
Drugi krok algorytmu modelowania, opisanego w podrozdziale 7.1.2, oparty jest na założeniu, że rozkład czasów oczekiwania F[ψ~, t] jest znany jawnie, tak że równanie (7.13) można rozwiązać względem τ zarówno analitycznie jak i numerycznie. Jeśli F[ψ~, t] jest znane, to opisany algorytm jest zazwyczaj bardziej efektywny, w szczególności, jeśli losowy krok czasu τ staje się duży i jeśli wyrażenie analityczne dla exp(–iH^τ ) jest również znane.
Jednakże w przypadku ogólnym ani ewolucja deterministyczna, ani rozkład czasu oczekiwania nie są znane jawnie. Czas oczekiwania τ powinien być określony wraz z numerycznym rozwiązaniem ewolucji deterministycznej. Aby to osiągnąć, rozpoczniemy od unormowanego stanu początkowego ψ~ i rozwiązujemy numerycznie równanie :
d/dt ψ(t) = –iH^ψ(t) (7.24)
które daje nieunormowany wektor stanu ψ(t) = exp(–iH^t )ψ~, którego norma zmniejsza się monotonicznie. Zgodnie z równaniem (7.13) czas oczekiwania τ może być określony przez sprawdzenie po każdym kroku całkowania kwadratu normy, zmniejszającego się do wartości η.
Zatem drugi krok algorytmu, może być zamieniony na następującą procedurę :
2) Bierzemy losową wartość η, rozłożoną równomiernie na odcinku [0, 1] i określamy nieunormowane rozwiązanie ψ(t) dla liniowego równania Schrödingera (7.24), odpowiadające wartości początkowej ψ(0) = ψ~. Czas oczekiwania τ jest wtedy określony przez warunek :
|| ψ(τ) ||2 = η (7.25)
Deterministyczna ewolucja pomiędzy skokami przyjmuje postać :
ψr(t + s ) = ψ(s)/ || ψ(s) || ; 0 ≤ s ≤ τ (7.26)
Taka forma algorytmu może być użyteczna przy realizacji modelowania w przypadku ogólnym. Dla rozwiązania numerycznego równania (7.24) pomiędzy skokami często dogodnie jest rozłożyć funkcje falową względem bazę własna hamiltonianu H^ i sprowadzić (7.24) do układu liniowego rr dla składowych s. Całowanie może być zrealizowane z wykorzystaniem standardowych procedur numerycznych , odpowiadających deterministycznym rr. Na zakończenie zauważymy, że podany powyżej algorytm, a w szczególności warunek (7.25) i równanie (7.26) mogą być również stosowane, jeśli generator w równaniu (7.24) zależy od czasu, tj. jeśli H^ = H^(t).
7.1.4 Wybór skoków.
Zgodnie z krokiem 3 opisanego algorytmu, powinniśmy wybrać skok kwantowy typu i z prawdopodobieństwem pi , zadanym przez równanie (7.15) i powinniśmy zrealizować odpowiednie przejście (7.16). Indeks i, znaczący skok, jest wielkością losową o wartościach w pewnym zbiorze indeksowym I = { 1, 2, .... , K } odpowiadającym rozkładowi pi.
W istocie istnieją dwa sposoby numerycznego przedstawienia wyboru skoku. Jeśli zbiór indeksowy I jest mały, to i może być wyznaczony z użyciem metody poszukiwania liniowego ( Gillepie 1976, 1992 ).
W tym celu interwał [0, 1] dzielimy na K podinterwałów Ji o długości pi i losujemy równomierną wielkość losową x ∈[0, 1]. Indeks i znajdujemy wtedy poprzez poszukiwanie interwału Ji, który zawiera x. Ważnym jest zauważyć, że prawdopodobieństwa pi powinny być większe niż zdolność rozdzielcza generatora liczb losowych stosowanego w takim algorytmie ( Hanusse, Blanche1981 )
Jeśli wartości prawdopodobieństw pi są wzajemnie równe, to może być bardziej efektywne wykorzystanie następującej metody. (Press, Flannery, Teukolsky, Vetterling 1992 )
Bierzemy dwie niezależne wielkości losowe x, y ∈ [0, 1]. Pierwsza z nich standardowo wykorzystywana jest dla wyboru indeksu i = int(Kx) + 1 ze zbioru indeksowego I. Jeśli druga wielkość losowa spełnia nierówność y < pi /pmax to taki indeks jest przyjmowany (pmax oznacza maksimum pi ). W przypadku przeciwnym i odrzucamy powtarzając procedurę.
Przewagą metody odrzucania nad poszukiwaniem liniowym jest to, że jest ona numerycznie bardziej efektywna, nawet jeśli liczba K elementów w zbiorze wskaźnikowym staje się duża. Odpowiedni warunek polega na tym, że iloczyn pmax i K nie powinien być większy niż 1.
Poszukiwanie liniowe i metoda odrzucania są nieefektywne, jeśli zbiór indeksowy I jest bardzo duży i jeśli rozkład pi jest bardzo niejednorodny. W tych przypadkach może być użyteczne zgrupowanie możliwych skoków w logarytmiczne klasy (Fricke, Schnakenberg 1991 ) lub wykorzystanie poszukiwania drzewiastego (Maksym 1988 ; Blue, Beichl, Sullivan 1995 )Jednakże czas, który jest zachowany dzięki ulepszonemu algorytmowi poszukiwania, może być znowu utracony przy aktualizacji struktur danych (Breuer, Huber, Petruccione 1996 )
7.2 Algorytmy dla stochastycznych równań Schrödingera
W niniejszym podrozdziale przeanalizujemy algorytmy modelowania stochastycznych równań Schrödingera (srs) typu dyfuzyjnego. Takie równania spotykaliśmy w podrozdziale 6.1.3 jako dyfuzyjne przybliżenia odpowiednich PDP w przestrzeni Hilberta. Pojawiają się one również w kontekście fundamentalnych modeli kolapsu dynamicznego (zobacz podrozdział 6.5.2 )
Dalej omówimy cztery różne algorytmy rozwiązywania srs : dobrze znany schemat Eulera, adaptowanane dla potrzeb srs schematy Heun’a i Runge’go-Kutty, oraz schemat drugiego rzędu małości dla równań różniczkowych zaproponowany przez Platen’a (Kloeden, Platen 1992 )
Podstawowy akcent kładziony jest na rząd zbieżności takich schematów używając przykład stochastycznego równania typu Schrödingera, spotykanego w podrozdziale 6.4.2 (Breuer, Dorner, Petruccione 2000 ).
Rozpatrzone metody mogą być oczywiście zastosowane dla analizy innych typów srs. W szczególności mogą być one łatwo uogólnione dla modelowania srs, przyjmując zespolone inkrementy w miejsce inkrementów rzeczywistych typu Wienera.
7.2.1 Ogólne uwagi o zbieżności.
Rozpatrzmy srs o następującej uogólnionej postaci Ito :
dψ(t) = D1(ψ(t)) dt + D2(ψ(t)) dW(t) (7.27)
o składowych : dryftowej :
i dyfuzyjnej :
gdzie wykorzystaliśmy umowne oznaczenie :
< A >ψ ≡ < ψ | A | ψ > (7.30)
gdzie dW(t) jest inkrementem rzeczywistym.
Równanie postaci (7.27) otrzymaliśmy już w podrozdziale 6.4.2 biorąc granicę dyfuzyjna dla detekcji homodynowej (zobacz równanie (6.181)). Dla uproszczenia rozpatrzymy przypadek pojedynczego operatora Lindblada A. Odpowiednia macierz gęstości zadana jest przez równanie :
Zauważmy, że równanie (7.27) jest srr z rosnącym szumem, co oznacza, że człon dyfuzyjny D2(ψ(t)), który zwiększa inkrement wienerowski, zależny jest od zmiennej stochastycznej ψ(t). Szum nazywa się szumem addytywnym, jeśli człon dyfuzyjny nie zależy od zmiennej stochastycznej.
Należy podkreślić, że ψ(t) jest procesem w przestrzeni Hilberta ℵS otwartego układu kwantowego. Zatem w ogólności reprezentuje on proces w nieskończenie wymiarowej przestrzeni Wektorowej. Jednakże w wielu interesujących przypadkach dyssypatywne zachowanie równań ruchu gwarantuje, że dynamika ogranicza się do skończonej
podprzestrzeni ℵS. Uzasadnia się to fizycznie koniecznością ograniczania się do przestrzeni skończenie wymiarowych tj. możemy rozpatrywać równanie (7.27) jako proces zachodzący w przestrzeni skończenie wymiarowej.
Teraz rozpatrzymy schemat numeryczny dla całkowania równania (7.27) na interwale czasu [ t0 ,t0 + T ].
Taki schemat prowadzi do dyskretnemu przybliżeniu po czasie ψk dla właściwego procesu ψk(t) przy czasach tk ≡ t0 + k∆t , k = 0, ... , n = T/∆t
Dyskretyzacja nie koniecznie musi być równo odcinkowa, jednakże dla uproszczenia tak właśnie przyjmiemy.
W niniejszym rozdziale ψk oznacza dokładny proces, określony przez równanie (7.27). analogicznie przez :
ρk = E[ | ψk > < ψk | ] (7.32)
- oznaczymy dyskretno- czasową aproksymacje macierzy gęstości, jednocześnie :
ρ(t) = E[ | ψ(t) > < ψ(t) | ] (7.33)
- jest dokładną macierzą gęstości, spełniająca równanie Lindblada.
Aby scharakteryzować zachowanie zbieżności schematów numerycznych, omawianych dalej, porównamy przybliżenie dyskretne ρk z rozkładem w szereg Taylora dokładnej macierzy gęstości ρ(t), który dany jest przez wyrażenie :
ρ(t + ∆t) = ρ(t) + £ρ(t) ∆t + ½ £2ρ(t)∆t2 + O(∆t3 ) (7.34) Bez utraty ogólności zawsze będziemy przyjmowali, że zadano zdeterminowany (tj. sztywny ) stan początkowy
ψ0 ≡ ψ(t0). Jednokrokowy błąd schematu numerycznego może być wyrażony przez różnicę :
ρ1 – ρ(t1) = O(∆tβ’ ) (7.35)
co oznacza, że schemat odtwarza rozkład w szereg Taylora ρ(t) włączając człony rzędu β’ – 1 do ∆t.
Łatwo jest pokazać, że całkowanie po ograniczonym interwale czasu [t0 , t0 + T ] zmniejsza rząd zbieżności o jeden, a ponieważ potrzebujemy n = T/∆t kroków czasowych w celu obliczenia gęstości w chwil czasu t0 + T = tn tj. :
ρn – ρ(tn ) = O(∆tβ ) (7.36)
z β = β’ – 1.
Jeśli równanie to spełnia schemat numeryczny, to będziemy mówili, że rozpatrywany schemat jest rzędu β.
Równanie (7.36) jest szczególnym przypadkiem słabej zbieżności rzędu b. Stopień przybliżenia stochastycznego równania różniczkowego może być również scharakteryzowany z użyciem pojęcia silnej zbieżności (Kloeden, Platen 1992 ). Jeśli algorytm numeryczny prowadzi do aproksymacji ψn dla ψ(tn) spełniającej zależność :
E[ || ψn – ψ(tn) || ] ≤ O(∆tβ ) (7.37)
Dla wystarczająco małego ∆t, to mówimy, że schemat jest zbieżny silnie w rzędzie β. Taka definicja nakłada warunek na bliskość wielkości losowych ψn i ψ(tn) na końcu interwału całkowania. W przeciwieństwie do tego, słaba zbieżność wymaga tylko tego iż rozkłady probabilistyczne ψn i ψ(tn) są wzajemnie bliskie, co jest najsłabszym kryterium.
W praktyce często interesujemy się tylko tą słaba postacią zbieżności, np. przy analizie przybliżenia funkcjonałów dla wielkości stochastycznej.
7.2.2 Schemat Eulera.
Schemat Eulera, jest prawdopodobnie najprostszym algorytmem rekursywnym dla obliczenia przybliżenia wektora stanu.
Ma on postać :
ψk+1 = ψk + D1(ψk)∆t + D2(ψk )∆Wk (7.38)
Wienerowskie inkrementy ∆Wk , k = 0, 1, 2, ... , n – 1sa niezależnymi losowymi wielkościami gaussowskimi o zerowej średniej matematycznej i dyspersją ∆t, indeks k pokazuje różne realizacje na każdym kroku. Zatem, możemy zapisać :
∆Wk = √∆t ξk (7.39)
gdzie ξk jest wielkością losową typu gaussowskiego o zerowej średniej matematycznej i jednostkowej dyspersji, tj. jest wielkością losową spełniającą rozkład standardowy.
Należałoby wspomnieć, że wykorzystanie wielkości losowych typu gaussowskiego nie jest konieczne, jeśli rzecz idzie o zbieżność równania w sensie (7.36).
Można wykorzystać i inne wielkości losowe, które pokrywają się tylko w pierwszym i drugim momencie z rozkładem standardowym.
Łatwo pokazać, że równanie (7.36) jest spełnione przy β = 1. Schemat Eulera, daje zatem słaby schemat zbieżności rzędu 1. Mimo słabego stopnia zbieżności, algorytm Eulera jest użyteczny, ponieważ jest on bardzo łatwy w
implementacji numerycznej i często daje wyniki o rozsądnym stopniu przybliżenia. W przypadku ogólnym, schemat ten nie zachowuje normy wektora stanu, w przeciwieństwie do srr (7.27). Dlatego też wektor stanu powinien być
normowany po każdej iteracji.
7.2.3 Schemat Heun’a.
Kolejna metoda, która będziemy nazywali schematem Heun’a jest uogólnieniem metody Heun’a, znanej z całkowania numerycznego deterministycznych rr. Jest on określony przez zależność rekursywną o postaci :
gdzie ∆Wk ponownie przyjmuje postać równania (7.39). W przypadku nie występowania dyfuzji metoda Heun’a staje się metoda dla deterministycznych rr. To w pewien sposób oznacza, że algorytm Heun’a prowadzi do zależności (7.36) z β = 1. Taki schemat rzędu 1 jest podobny do metody Eulera. To kontrastuje z srr z szumem addytywny, gdzie metoda Heun’a rzeczywiście zapewnia drugi rząd zbieżności.
7.2.4 Schemat Rungego-Kutty czwartego rzędu.
Dla deterministycznych rr istnieje cały szereg różnorodnych metod numerycznych i łatwo jest się pomylić przy heurystycznym ich wykorzystaniu w celu całkowania srr. Następujący przykład pokazuje, że ostrożność jest tutaj w istocie konieczna.
Dobrze znany numeryczny schemat dla deterministycznych rr, to metoda Rungego-Kutty czwartego rzędu.
Heurystyczna modyfikacja dla przybliżonego rozwiązania srr może być otrzymana, jeśli na każdym kroku po czasie dryft całkujemy metoda Rungego-Kutty, a składowa dyfuzyjna srr zgodnie ze schematem Eulera. To prowadzi do
następującego schematu :
Rząd zbieżności ponownie może być określony poprzez obliczenie realizacji w równaniu (7.36). Zadziwiającym
wynikiem jest, to że ten schemat jest rzędu β = 1.
Stochastyczne uogólnienie metody Rungego-Kutty, jest rzędu 1 w przeciwieństwie do metody Rungego-Kutty, wykorzystywanej dla deterministycznych równań falowych, gdzie rząd zbieżności wynosi 4.
Z powyższego przykładu widać, że ogólne proste uogólnienia schematów Rungego-Kutty wysokiego rzędu
niekoniecznie prowadzą do wysokich rzędów dla srr. Podana powyżej metoda jest oczywiście zbieżna, ale w przypadku ogólnym nie będzie ona posiadała lepszych cech niż np. metoda Eulera.
7.2.5 Schemat drugiego rzędu małości.
Ostatni algorytm, który rozpatrzymy został zaproponowany przez Platena (Kloeden, Platen 1992 )W naszej formie zapisu taki schemat może być przedstawiony poprzez zależność rekursywną o postaci :
gdzie inkrementy ∆Wk powinny spełniać określone warunki.
Warunki takie będą spełnione, jeśli spełnione będzie równanie 97.39). Dla przypadku nie występowania szumu schemat sprowadza się do schematu Heun’a dla deterministycznych rr.
Ponownie realizując analizę błędów, możemy znaleźć iż schemat Platena posiada rząd zbieżności β =1.
Zatem, w przeciwieństwie, które przedstawiliśmy wcześniej, w/w schemat rzeczywiście jest schematem wysokiego rzędu w sensie określonym powyżej. Fakt ten zilustrowano w przedstawionych dalej przykładach.
7. Przykłady.
W niniejszym podrozdziale zastosujemy algorytmy modelowania dla PDP i stochastycznych równań Schrödingera ku pewnym wcześniej sygnalizowanych przykładów.
7.3.1 Tłumiony oscylator harmoniczny.
7.3.1.1 Modelowanie PDP.
Przeanalizujemy teraz następujący proces dla tłumionego oscylatora harmonicznego (zobacz podrozdział 6.7 ) :
Zamodelujemy to równanie dla najprostszego warunku początkowego, tj. dla stanu o zadanej liczbie cząstek ψ(0) = | n >.
W takim przypadku deterministyczne okresy ewolucji są trywialne ponieważ :
exp(–iH^t ) | n > = exp(–γnt/2 ) | n > (7.51)
Również skoki trywialnie zadane są przez warunek :
| n > → a | n > / || a | n > || = | n – 1 > (7.52)
Ponieważ rozkład czasu oczekiwania ma postać :
to znajdujemy, że losowy czas oczekiwania może być otrzymany z równania :
τ = –(1/γn) ln(η) (7.54)
gdzie η jest równomiernie rozłożoną na interwale [0, 1] wielkością losową.
Rysunek 7.1 pokazuje realizacje procesu dla stanu początkowego ψ(0) = | n0 > = | 9 >. W charakterze przykładu obliczenia średnich matematycznych pokazujemy na rysunku 7.2 średnie operatora liczby cząstek, określone z próbki z R = 100 realizacji (zobacz (7.3)) :
Analityczne wyrażenie dla takiej wielkości ma postać :
< a†a(t) > = n0 exp(–γt ) (7.56)
Zgodnie z równaniem (7.4) błąd standardowy dla takiej wartości średniej może być przedstawiony przez wyrażenie :
Wartości dla błędu standardowego pokazano na rysunku 7.2 przez koordynaty błędów.
Rys. 7.1 Jedna realizacja procesu (7.50) dla tłumionego oscylatora harmonicznego z warunkiem początkowym
| n0 > = | 9 > (linia ciągła ). Linia przerywana reprezentuje analityczne wyrażenie (7.56) dla średniej matematycznej operatora liczby cząstek.
Rys. 7.2 Uśrednienie po 100 realizacjach dla tłumionego oscylatora harmonicznego z warunkiem początkowym
| n0 > = | 9 >. Linia przerywana reprezentuje wynik analityczny.
7.3.1.2 Modelowanie stochastycznego równania Schrödingera.
Stochastyczne równanie Schrödingera (7.27) dla tłumionego oscylatora harmonicznego ma postać :
W modelowaniu numerycznym wektor stanu reprezentuje się w fokowskiej bazie stanów | n >. Ponownie stan
początkowy bierzemy w postaci ψ(0) = | n0 = 9 >, a przestrzeń Hilberta obcinamy do nmax = 12, tj. modelowanie zrealizowano w podprzestrzeni o wymiarze N = 13.
Oddzielna realizacja procesu przedstawiona jest na rysunku 7.3. Linia ciągła pokazuje krzywą analityczna (7.56). Linia cienka pokazuje rozwiązanie otrzymane z uwzględnienie drugiego rzędu małości w schemacie Platena, linia gruba to rozwiązania otrzymane z użyciem innych metod, które nie są wymienione w niniejszym przykładzie. Wielkość kroku czasowego ∆t = 0,01. Zauważmy, że jeden i ten sam ciąg liczb pseudolosowych został wykorzystany dla każdego przebiegu, po to aby podkreślić różnica lub podobieństwa pomiędzy omawianymi metodami.
Aby przeanalizować charakter zbieżności, obliczamy wartość średnią (7.55) operatora liczby cząstek przy ustalonym czasie T dla R = 10 realizacji i dla różnych kroków czasowych ∆t. Wyniki takiej analizy pokazano na rysunku 7.4.
Koordynaty błędów nie są narysowane na tym rysunku, ponieważ miałyby one wielkości porównywalne z wymiarami symboli graficznych. Oczywiście, platenowski schemat drugiego rzędu małości jest zbieżny bardzo dobrze już dla całkiem dużych rozmiarów kroku (porównując go do innych schematów ).
Heurystyczna metoda Rungego-Kutty pokazuje najgorsza charakterystykę zbieżności.
Rys. 7.3 Oddzielna realizacja średniej dla operatora liczby cząstek dla tłumionego oscylatora harmonicznego. Linia cienka odpowiada słabemu schematowi drugiego rzędu, linie grube schematom Eulera, Heun’a i Rungego-Kutty ( nie rozróżnialne w skali rysunku ). Linia gładka reprezentuje wynika analityczny (7.56).
Rys. 7.4 Wykres średniej liczby fotonów M^ analitycznie obliczonej wartości < n(T) > przy γT = 1,2 w zależności od rozmiaru kroku ∆t dla różnych metod. Linie ciągłe odpowiadają liniowemu i kwadratowemu dopasowaniu.
Rys. 7.5 Czas CPU (unormowany do jedności ) w zależności od błędu absolutnego dla danych z rysunku 7.4.
Mając wiadomości o rzędzie zbieżności, można dobrać funkcje liniowe ( w przypadku metod Eulera, Heune’a Rungego- Kutty ) lub funkcje kwadratowe (w przypadku metod platenowskich ) dla dopasowania punktów danych. Odpowiednie obiekty (pokazane również na rysunku 7.4) są dobrze zgodne z punktami danych i potwierdzają wyniki, otrzymane w podrozdziale 7.2.1
Innym kryterium oceny algorytmów numerycznych jest czas CPU. Rysunek 7.5 pokazuje czas CPU( unormowany do jedności ) w porównaniu z błędami dla danych pokazanych na rysunku 7.4. widzimy, że schemat Platena demonstruje najlepsze własności, po nim jest schemat Eulera.
Rys. 7.6 Jedna realizacja PDP (7.59) dla wzbudzanego dwupoziomowego atomu. Linia ciągła pokazuje ρ11 obliczone z realizacji. Linia przerywana reprezentuje rozwiązanie analityczne dla ρ11, zgodnie z równaniem (3.289).
Parametry : Ω = γ = 0,4.
Rys. 7.7 Uśrednienie po 1000 realizacjach PDP (7.59) dla wzbudzanego dwuatomowego atomu.
Linia ciągła pokazuje ρ11, otrzymane z uśrednienia po realizacjach. Linia przerywana reprezentuje rozwiązanie analityczne dla ρ11 odpowiadające równaniu (3.289). Błąd standardowy wartości średniej wskazywany jest przez znaczniki błędów. Parametry – te same co dla rysunku 7.6.
7.3.2 Wzbudzany dwupoziomowy układ
W charakterze drugiego przykładu rozpatrzmy PDP zadany przez równanie :
Jak pokazano w podrozdziale 6.3 równanie to opisuje bezpośrednią fotodetekcje wzbudzanego układu (zobacz (6.129)), jednocześnie stochastyczne równanie Schrödingera :
odpowiada detekcji homodynowej (zobacz podrozdział 6.4)
W modelowaniu numerycznym rozpoczynamy od atomu w stanie podstawowym | g > i przyjmujemy prawdopodobieństwo :
ρ11(t) = < e | ρ(t) | e > (7.61)
znalezienia atomu w stanie wzbudzonym | e >. Z próbki realizacji takiego prawdopodobieństwa oceniamy determinowana średnią :
Błąd standardowy σ^t zadany jest z zależności ;
Analizę algorytmu dla generacji realizacji procesu kawałkami –deterministycznego można skrócić, ponieważ zebraliśmy już wszystkie odnoszące się ku niemu zależności (rozkład czasu oczekiwania i skoków ) w podrozdziale 6.3.
Zgodnie z równaniami (7.13) i (6.133) losowy czas oczekiwania τ może być określony równaniem :
dla τ. Dla znalezienia pierwiastka tego równania wykorzystaliśmy standardową procedurę numeryczną.
Na rysunku 7.6 pokazaliśmy | e | ψ(t) > |2 dla jednej realizacji PDP. Odpowiednie średnie po ansamblu z 1000 realizacji pokazano na rysunku 7.7, który to jasno pokazuje zbieżność użytego algorytmu.
Teraz zwrócimy uwagę na modelowanie stochastycznego równania Schrödingera (7.600. Oddzielne realizacje obliczone na podstawie różnych metod numerycznych, omówionych w podrozdziale 7.2 pokazano na rysunku 7.8.
Na rysunku wstawkowym, widać że słaby schemat rzędu 2 różni się od innych schematów, które z trudem rozróżniamy na większej skali.
Wyniki modelowania z rożnymi rozmiarami kroku dla R = 5 • 105 realizacji na punkt pokazano na rysunku 7.9. Dla jasności pokazujemy interwały błędów tylko w tych punktach, które nie są bardzo bliskie siebie. Zostały one policzone zgodnie z równaniem (7.63). W dowolnym z czterech przypadków wyniki numeryczne są dobrze zgodne z analitycznymi przewidywaniami. Metoda Rungego-Kutty wydaje się lepsza co do zbieżności, niż schematy Eulera lub Huena dla dużych rozmiarów kroku. Jednakże, jeśli przyjąć do wiadomości odpowiednie czasy CPU (tak jak to zrobiono na rysunku 7.10), staje się oczywiste, że metody te posiadają mniej lub więcej równy stosunek prędkość zbieżności – dokładność, jednakże schemat platenowski ponownie przejawia najlepsze własności.
Sytuacja zmienia się, jeśli wybieramy parametry takie, że γ << Ω, tj. jeśli deterministyczna część stochastycznego równania Schrödingera ma przewagę nad częścią szumową.
Metoda Rungego-Kutty pokazuje najlepsza charakterystykę, jednocześnie schemat Platena ma takie same własności jak schemat Huena. Oczywiście, tak jak tego oczekiwaliśmy dla γ << Ω, deterministyczna część równania ruchu dominuje i dlatego metoda Rungego-Kutty czwartego rzędu staje się bardziej efektywna.
Rys. 7.8 Oddzielna realizacja dla ρ11, otrzymane z równania (7.60), zgodnie z innymi schematami z podrozdziału 7.2 z krokiem ∆t = 0,01. Wyniki nie są rozróżnialne w skali rysunku. Linia gładka daje rozwiązanie analityczne (3.289).
Wstawka pokazuje powiększenie małej części rysunku (linia cienka – schemat Platena, linia gruba – schemat Eulera, Huena i Rungego- Kutty )
Parametry : Ω = γ = 0,4.
Rys. 7.9 Różnica pomiędzy oceną średniej | < e | ψr(T) > |2 i dokładnymi wartościami dla różnych wielkości kroków i
Rys. 7.10 Czas CPU (unormowany do jedności ) w zależności od błędu absolutnego dla danych punktów z rysunku 7.9.
7.4 Analiza realizacji numerycznej.
Algorytm M-C jest często sposobem analizy układów o wysokim wymiarze, dla których obliczenia deterministyczne oparte na numerycznym całkowaniu równań ruchu, znajdują się daleko poza możliwościami dowolnego komputera.
Z drugiej strony, metody deterministyczne mogą być użyteczne dla układów niskowymiarowych.
W niniejszym podrozdziale chcemy przedstawić systematyczną analizę numerycznej efektywności metody funkcji falowych M-C i porównać ja z całkowaniem odpowiedniego równania dla macierzy gęstości.
Najbardziej interesującą jest zależność czasu rozwiązania danego zagadnienia od wymiaru układu N. Pokażemy, że czas potrzeby CPU dla obu w/w podejść może być wyrażony w języku prostych praw potęgowych.
7.4.1 Efektywność numeryczna i prawa skalowania.
Wyniki modelowania z użyciem metody M-C niosą ze sobą zawsze błędy statystyczne. Takie błędy w określony sposób są związane z liczba realizacji, które ś generowane w modelowaniu. Jednakże sumaryczny czas CPU, wymagany dla modelowania M-C, zależy krytycznie od żądanej dokładności. Typowym zagadnieniem analizy numerycznej może być np. „ocena średniej matematycznej energii z błędem względnym 1%” lub też „obliczenie macierzy gęstości z
dokładnością lepsza niż 10–4 w każdym z jej elementów”.
Wtedy to nieuchronnie pojawia się pytanie, które wymaga odpowiedzi : przy jakich okolicznościach modelowanie stochastyczne jest preferowane ?
Wielkością krytyczną dla względnej charakterystyki dokładności jest liczba realizacji R. Zależność błędów modelowania od wymiaru zagadnienia zasługuje na szczegółowe omówienie. W tym celu oznaczymy N jako liczbę zmiennych zespolonych, które są wykorzystywane dla numerycznego przedstawienia funkcji falowej, tj. liczbę stanów bazowych.
Zatem, liczba zmiennych zespolonych, które są wymagane dla przedstawienia macierzy gęstości wynosi ½ N Zależność (7.9) pozwala nam zapisać następujące równanie :
σ^t2 = λB(N)/R (7.65)
w którym czynnik λB(N) uwzględnia zależność błędu statystycznego od obserwabli B, oraz od rozmiaru układu N, ale nie zależy od rozmiaru próbki R. Przy wykorzystaniu dostatecznie dużej próbki realizacji wielkość λB(N) może być określona przez równanie (7.65) dla danych modelowych. Wtedy to równanie (7.65) może być rozwiązane dla R :
R ≡ R(N) = λB(N)/ σ^t2 (7.66)
Jest to liczba realizacji, która jest wymagana dla osiągnięcia danej dokładności σ^t dla obserwabli B i rozmiaru układu N. Jeśli λB(N) zmienia się jak funkcja potęgowa od rozmiaru układu :
λB(N) ~ N–x
to można przyjąć następująca klasyfikacje (Ferrenberg, Landau, Binder 1991 ) :
1) jeśli x = 1, to obserwabla B jest silnie samouśredniona 2) jeśli 0 < x < 1, to obserwabla B jest samouśredniona 3) jeśli x = 0, to obserwabla B nie jest samouśredniona
W ramach takiego podejścia możemy przeanalizować charakterystykę skalowania czasu CPU, wymagana przez deterministyczne i stochastyczne podejście dla numerycznej oceny średnich matematycznych z zadaną dokładnością.
W przypadku podejścia w oparciu o równanie dla macierzy gęstości część komend numerycznego całkowania, które dominują przy zapotrzebowaniu czasu CPU, związanych jest z obliczaniem generatora £ρ równania Lindblada. Jedno takie obliczenie dla dostatecznie dużego N wymaga nakładu czasu CPU, proporcjonalnego do potęgi N, a czas CPU wymagany dla całkowania macierzy gęstości na zadanym fizycznym interwale czasu w głównym rzędzie po N ma postać :
TDME = k1s1(N) Nβ (7.67)
Gdzie s1(N) – jest liczbą razy, która powinna być obliczana wielkość £ρ ; k1i β - zależne są od typu w/w zagadnienia, ale nie od N.
Oprócz tego, k1 w szczególności, zależy od formy przedstawienia danych na komputerze.
Analogicznie w wielu przypadkach z sztywne ograniczoną po czasie częścią modelowania stochastycznego jest obliczenie generatora H^ψ, a wymagany czas CPU dla takiego modelowania ma postać :
TStS = k2R(N)s2(N) Nα (7.68)
Gdzie R(N) – jest liczbą realizacji procesu, która jest generowana przy obliczeniach dla układu o wymiarze N, s2(N) – jest liczbą obliczeń H^ψ dla danej realizacji, k2 jest wielkością analogiczną do k1.
W wielu sytuacjach s1(N) i s2(N) są w przybliżeniu równe. Przy warunku, że wykorzystywana jest podobna procedura całkowania numerycznego, będzie to przypadek, kiedy najmniejsza skala czasu dla dynamiki funkcji falowej jest w przybliżeniu równa takiej że skali dla macierzy gęstości.
Ponieważ chcemy oddzielić efekty rozmiaru układu od zjawisk dynamicznych, będzie to najbardziej interesujący nas przypadek, a przykład podany w podrozdziale 7.4.2 ilustruje ten właśnie przypadek.
Tylko zasygnalizujemy, że istnieją również takie sytuacje, w których s1(N) i s2(N) są bardzo różne. W ogólności ich stosunek może zależeć od N, a wielkości te nie koniecznie rosną wraz z czasem fizycznym w ten sposób jak maiło to miejsce w analizowanym przez nas przypadku.
Dalej rozpatrzymy przypadek w którym skala czasu dynamiki oddzielnej realizacji procesu stochastycznego zmienia się w czasie ewolucji czasowej w szerokich granicach np. ma to miejsce w chłodzeniu laserowym (podrozdział 8.3) Modelowanie jednej realizacji będzie zawierało siatki z bardzo długimi krokami czasu, przerywane fazami z szybkimi fazami ewolucji i krótkimi krokami czasu. Z drugiej strony, integrator równania macierzy gęstości, który opisuje dynamikę całego anasamblu, powinien być zawsze aproksymowany do skali małego czasu. Jest jasne, ze w takich przypadkach metoda stochastycznej funkcji falowej jest dogodniejsza.
Liczba operacji zmiennoprzecinkowych wymaganych dla obliczenia £ρ i obliczenia H^ψ różnią się czynnikiem N i jak tego oczekiwaliśmy :
β ≈ α + 1 (7.69)
Zależność ta będzie sprawdzona w dalszym przykładzie.
Wnioskując na podstawie powyższych rozważań, możemy zapisać równania (7.67) i (7.68) w krótkiej formie :
TDME = k1Nα+1 (7.70)
TStS = k2 Nα–x (7.71)
Otrzymujemy tutaj to, że liczba kroków s1 i s2 jest w przybliżeniu równa i może być ona włączona do stałych α, k1 i k2.
Reprezentacja metody stochastycznej funkcji falowej w porównaniu z całkowaniem macierzy gęstości może być oceniona przez pomiar różnicy pomiędzy wykładnikami potęg, różnica ta wynosi 1 w przypadku, kiedy nie ma samouśrednienia (x = 0 ), oraz 2 w przypadku silnego samouśrednienia ( x = 1 ).
7.4.2 Tłumiony pobudzany oscyaltor Morse’a (* The damped driven Morse oscillator *)
Prawa skalowania (7.70) i (7.71) zilustrujemy dla czasowych zapotrzebowań CPU. W tym celu rozpatrzymy nietrywialny przykład, a mianowicie tłumiony oscylator Morse’a, który jest pobudzany zależną od czasu siłą.
W takim przykładzie przekonamy się, że prawa skalowania są bardzo dobrze wypełnione w takim nieliniowym zagadnieniu.
7.4.2.1 Opis modelu.
Hamiltonowski opis cząstki koherentnej ma postać :
HS(t) = HM + HL(t) (77.2)
gdzie
HM = (1/2m)p2 + V(q) (77.3)
a
V(q) = D[ 1 – exp(–bq )]2 (7.74)
jest potencjałem Morse’a.
Hamiltonian Morse’a (7.73) może być wykorzystany dla modeli np. molekularnego stopnia swobody w granicach oddzielnej elektronowej powierzchni energii potencjalnej. Przy odpowiednim wyborze parametrów w jego ramach otrzymujemy przejrzysty realistyczny opis np. wibracyjnej dynamiki lokalnego wiązania O-H w molekule wody.
Zewnętrzne zaburzenie opisujemy członem oddziaływania :
HL(t) = µqF0s(t) sin(ωLt ) (7.75)
gdzie µq – jest odpowiednia składowa molekularnego momentu dipolowego, F0 – maksymalna siła pola, s(t) – obwiednia pędu wzbudzenia.
Dla konkretyzacji przyjmiemy :
s(t) = sin2(πt/tp ) (7.76)
gdzie tp reprezentuje czas trwania impulsu.
Fizyczne omówienie takiego modelu można znaleźć w zbiorze artykułów pod redakcja Manza i Woestea (1995)
Hamiltonian Morse’a HM posiada ograniczona liczbę N stanów związanych, która zależy od parametrów potencjału Morse’a V(q). Nasza analiza będzie ograniczona do przypadku, kiedy dynamika jest dobrze określona przez ograniczony wektor stanu oscylatora Morse’a. Z tego powodu modelowania będzie przeprowadzone w energetycznej bazie własnej | j > ; j = 0, 1, 2, ... , N – 1 dla HM z odpowiadającymi tej bazie wartościami własnymi energii Ej (znanych analitycznie ).
Operatory skoku maja postać :
Ajk = | j > < k | ; j, k = 0, 1, 2, .... N – 1 (7.77)
z odpowiadającymi im intensywnościami γjk. Takie intensywności przyjmujemy jako zadane dla słabo sprzężonego równania podstawowego, opisującego dysypatywną dynamikę w rezerwuarze termicznym. Operatory Ajk są operatorami własnymi HM, odpowiadającymi wartościom własnym ( Ej – Ek ). Efekt działania operatora skoku Ajk na funkcje falową układu może być zatem interpretowany jako opis przejścia promienistego lub absorpcji drgań energii kwantowej
| Ej – Ek |. W reprezentacji energetycznej podstawowe równanie Markowa dla elementów macierzowych ρjk = < j | ρS | k > można zapisać następująco :
Odpowiedni PDP kierowany jest przez poniższe równanie dla składowych ψj = < j | ψ > unormowanego wektora stanu ψ :
W powyższych równaniach f(t) jest niezależną od czasu zewnętrzną siłą ( porównaj z (7.75)) :
f(t) = µF0s(t) sin(ωLt) (7.80)
Qjk = < j | q | k > - są to elementy macierzowe operatora dipolowego.
Macierz (Qjk ) jest macierzą rzeczywista i symetryczną. Γj – jest sumaryczna intensywnością wszystkich skoków z | j > : N–1
Γj = Σγmj (7.81)
m=0
Na koniec – poissonowskie inkrementy spełniają warunki :
7.4.2.2 Wyniki modelowania.
Dalej rozwiązujemy równania 97.78) i (7.79) i porównujemy charakterystyki numeryczne dla obu tych równań. Za stan początkowy zawsze przyjmuje się stan podstawowy oscylatora Morse’a. Aby przeanalizować wpływ rozmiaru układu na zapotrzebowanie czasu programów numerycznych, analizuje się szereg podobnych oscylatorów z różna liczbą N stanów związanych. Rozmiar układu N zmienia się w granicach N = 12, ... , 18. Aby nie przesadzić, liczba wymaganych stanów faktycznej dynamiki, zawsze znajduje się w granicach tejże liczby stanów, odpowiednio do tego należy również przeskalować pole wzbudzenia. W dalszej analizie uczestniczą dwie kombinacje parametrów : jedna odpowiadająca słabemu tłumieniu, druga – silnemu. Zauważmy również, że dla scałkowania równania dla macierzy gęstości i dla deterministycznych kawałków PDP wykorzystuje się ten sam schemat Rungego-Kutty.
Szczegóły wykorzystywanych parametrów fizycznych i reprezentacje numeryczna opisano przez Breuera, Hubert’a i Petruccione (1997)
W pierwszej kolejności należy zwrócić uwagę na wykładniki α, β, które wprowadzono w równaniach (7.67) i (7.68). Na rysunku 7.11 pokazano czas CPU na jeden krok czasu całkowania numerycznego, jako funkcje rozmiaru układu N.
Pomiar nachylenia linii daje szacunek :
β = 3,0 ± 1 ; α = 2, 0 ± 0,1 (7.84)
co potwierdza sensowność równania (7.69). Takie wartości można łatwo zrozumieć – w przypadku modelowania stochastycznego największą część zapotrzebowania czasu, pochłania przemnożenie ψ z macierzą dipolowa Q (porównaj z (7.79)), które wymaga O(N2 ) operacji zmiennoprzecinkowych. Analogicznie dla całkowania macierzy gęstości obliczenia prawej części równania (7.78) dla wszystkich j i k zawierają O(N3 ) operacji zmiennoprzecinkowych.
rys. 7.11 Czas CPU na jeden krok przy całkowaniu równia dla macierzy gęstości (∆) i dla propagacji wektora falowego ψ( ). Linie przerywane pokazują jak zostały obliczone wykładniki potęg α i β, które zostały wyliczone z współczynnika nachylenia. Linie ciągłe łącza po prostu punkty danych.
Rysunek 7.12 pokazuje liczbę kroków całkowania s1 i s2, które są konieczne dla obliczenia pędu całkowitego. Jak można zobaczyć s1 i s2 rosną wraz ze wzrostem N, ale pozostają w przybliżeniu równe sobie. Ich wzrost ma miejsce dzięki szczególnym własnością omawianego przykładu. Oba układy rr, które powinny być rozwiązane dla obliczenia macierzy gęstości i dla modelowania stochastycznego, staja się bardziej stabilne wraz ze wzrostem N.
Rys. 7.12 Trójkąty ∆ pokazują s1liczbę kroków integratora dla obliczenia jednego pędu, z wykorzystaniem macierzy gęstości. Romby pokazują s2 średnią wartość kroków dla obliczenia jednej realizacji procesu stochastycznego.
Rys. 7.13 Funkcja λH(N) mierzy własność samouśrednienia średniej energii oscylatora Morse’a HH ( porównaj z (7.65)). Trójkąty odpowiadają słabej dyssypacji, romby – silnej dyssypacji.
Aby zbadać charakterystykę R(N), liczby realizacji M-C, które powinny być wygenerowane dla analizy układu o rozmiarze N, należy rozpatrzyć błąd standardowy obserwabli układu.
Rozpatrzmy w tym celu np. energię HM oscylatora Morse’a, której błąd standardowy oznaczymy jako σ^H. Zgodnie z równaniem (7.66) otrzymujemy :
R(N) = λH(N)/ σ^H2 (7.85)
Na rysunku 7.13 pokazano zachowanie funkcji λH(N) dla różnych zbiorów parametrów otoczenia. Z dokładnością do fluktuacji statystycznych mamy λ(N) i dlatego R(N) jest nierosnącą funkcją od N. Podobne wyniki otrzymujemy dla innych obserwabli i elementów macierzowych ρ. To oznacza, że aby osiągnąć stały błąd statystyczny w wynikach modelowania, kiedy rozmiar układu N zwiększa się, liczba realizacji nie powinna się zwiększać. Stad wynika, że w ostatecznym rozrachunku modelowanie stochastyczne dla układów o dużych rozmiarach będzie szybsze niż rozwiązanie równania dla macierzy gęstości.
Wniosek ten zilustrowano przykładem na rysunku 7.14. Wykresy pokazują wymagany czas CPU dla całkowania
równania dla macierzy gęstości i dla generowania takiej liczby realizacji procesu stochastycznego, jaka jest potrzebna dla otrzymania błędu standardowego σ^H = 4 • 10–3 dla średniej energii oscylatora. Liczba R realizacji została wyznaczona z pomocą równania (7.85). Zgodnie z rysunkiem 7.13 wybieramy λH(N) = 10–3 niezależnie od N. Krzywe odpowiadają różnym zależnościom potęgowym i przecinają się w pewnym punkcie N0. W przedstawionym przykładzie :
N0 ≈ 35 dla przypadku słabej dyssypacji i N0 ≈ 55 dla silnej dysypacji. Przy dużym N0 modelowanie stochastyczne przebiega szybciej.
Podstawowym wynikiem modelowania numerycznego jest, to że czas CPU dla całkowania macierzy gęstości TDME i czas dla modelowania stochastycznego funkcji falowej TStS skaluje się z rozmiarem układu N tak :
TDME ~ N3 (7.86)
TStS ~ N3 (7.87)
Chociaż analiza numeryczna została przeprowadzona na szczególnym przykładzie, to wyniki osiągnięte w podrozdziale 7.4.1 są znacznie ogólniejsze. W szczególności, podczas gdy wartości absolutne wykładników w równaniach (7.86) i (7.87) zależą od konkretnych własności analizowanego układu, to ich różnica ma ogólniejsze znaczenie.
Przy odpowiednio ogólnych warunkach wykładnik w wyrażeniu dla TDME przewyższa analogiczna wartość wykładnika TStS o wielkość leżąca w przedziale od 1 do 2.
Wielkości TDME i TStS które analizowaliśmy powyżej, reprezentują czas w przeciągu którego programy pracują na oddzielnym procesorze. Przy porównaniu efektywności kodów numerycznych, które zostały przygotowane dla pracy na komputerach równoległych (* obecnie procesorach wielordzeniowych *) istnieją inne ważne kryteria :
Przyspieszenie i skalowalność. Przyspieszenie definiuje się jako stosunek miedzy czasami, wymaganymi dla
zrealizowania pracy na oddzielnym procesorze i na komputerze równoległym. Pod skalowalnością mamy na myśli to, że w szerokim zakresie liczby procesorów przyspieszenie jest bliskie liczbie procesorów maszyny równoległej. To oznacza, że czas tracony na współpracę pomiędzy procesorami i nakłady na ich synchronizację są małe.
Zauważmy, że efektywne skalowane równoległe procedury całkowania macierzy gęstości są zagadnieniem bardzo złożonym. Jednocześnie, metoda stochastycznych funkcji falowych jest w swej naturze algorytmem równoległym i skalowalnym. A ponieważ oddzielne realizacje próbek są generowane niezależnie synchronizacja algorytmu jest wymagana tylko dla wynikowego uśrednienia lub archiwizacji i dla celów kontrolnych.
Metody M-C, są zatem idealne dla algorytmów równoległych.
Rys. 7.14 CPU- czas wymagany dla całkowania równania macierzy gęstości (trójkąty ) i dla wygenerowania takiej liczby realizacji procesu stochastycznego (romby ), która jest wymagana dla otrzymania błędu standardowego
σ^H = 4 • 10–3 dla średniej energii oscylatora. Wykresy po prawej stronie obejmują pełna wariacje rozmiaru układu, wykresy po lewej – dają powiększony obraz przy mniejszym N.
Literatura.
************************************************************************************************
Rozdział 8 Stochastyczne metody modelowania.
W rozdziale 6 zaprezentowaliśmy szereg prostych przykładów dla procesów stochastycznych w przestrzeni Hilberta, jak również ich rozwinięcie w ramach teorii pomiaru ciągłego. Niniejszy rozdział poświęcony jest omówieniu teorii w oparciu o mikroskopowe równania elektrodynamiki kwantowej QED. Akcent kładziemy na pochodzenie dowolnych operacji kwantowych i odpowiednich procesach stochastycznych dla realizacji schematów detekcji, które to wynikają bezpośrednio z hamiltonianu, który opisuje oddziaływanie stopni swobody materii z kwantowym polem EM. We wszystkich przypadkach rozdział ten będzie rozpatrywał stopnie swobody materii w sposób nierelatywistyczny, przyjmując do wiadomości zastosowania w optyce kwantowej.
Po krótkim omówieniu zagadnienia kwantowania pola EM, rozpoczynamy od wyprowadzenia podstawowego wyrażenia dla QED- operacji, która to opisuje odwrotne działanie na stopnie swobody materii, indukowane przez pomiar
zmiennych polowych. Z takiego podstawowego wyrażenia ustanowimy reprezentacje stochastyczną dla promieniowania multipolowego danej materii. Stochastyczna dynamika będzie sformułowana w języku stochastycznego równania ruchu dla zredukowanej macierzy gęstości źródła.
W przypadku pełnego pomiaru liczb kwantowych, emitowanych fotonów, otrzymujemy stochastyczna ewolucje dla wektora stanu źródła. Będziemy również omawiali pomiary niepełne. Prowadza one do dynamiki w postaci równania stochastycznego dla macierzy gęstości, które nie zachowuje czystości stanów kwantowych.
Niniejszy rozdział zawiera omówienie pewnych ważnych zagadnień fizycznych, dla których wzajemne powiązanie koherentnej ewolucji kwantowej i procesów dysypatywnych odgrywa ważną rolę. Na początku przeanalizujemy pojawienie się stanów ciemnych przy oddziaływaniu atomów z koherentnymi polami laserowymi. Oprócz tego, przeanalizujemy mechanizm koherentnego wychwytu przy ochładzaniu atomów w reżimie subodbicia w polach laserowych. Wyprowadzimy przy tym odpowiedni proces kawałkami deterministyczny dla kwantowego ruchu atomów.
(* In addition, we study the mechanism of coherent population trapping in the sub-recoil cooling of atoms
in laser fields. *) To da nam interesujący przykład pojawienia się długo zakresowego rozkładu Levy’ego, wynikłego z efektu kwantowej interferencji.
Na koniec, przeanalizujemy zjawiska dyssypatywne w dynamice układów w silnych periodycznych polach wzbudzenia.
Odpowiednia dla nich stochastyczna dynamika funkcji falowej może być wprowadzona z wykorzystaniem reprezentacji wektora stanu zredukowanego układu w bazie składającej się ze stanów Floquet’a.
8.1 Pomiary ciągłe w QED.
8.1.1 Konstrukcja mikroskopowego hamiltonianu.
W pierwszej kolejności przypomnimy podstawowe zasady kwantowania swobodnego pola promieniowania
(Bjorken, Drell 1965 ), które może być przedstawione przez swobodny operator polowy A(x, t ) dla wektora potencjału.
Nakładając warunki cechowania Coulomba :
∇ • A(x, t ) = 0 (8.1)
operatory pola elektrycznego i magnetycznego zadane są następująco ( W niniejszym podrozdziale jawnie zapisujemy wszystkie stałe fizyczne takie, jak stała Plancka h, prędkość światła c, ładunek elektronu e i jego masa m. W dalszej kolejności wykorzystywać będziemy układ jednostek Gaussa, tak więc stała struktury subtelnej ma postać
α = e2/hc ≈ 1/137 )
E(x, t ) = – (1/c) ∂/∂t A(x, t ) (8.2)
B(x, t ) = ∇× A(x, t ) (8.3)
Operatory polowe są to operatory Heisenberga, spełniające równanie falowe :
(( 1/c2) ∂2/∂t2 –∇ )A(x, t ) = 0 (8.4)
Dla dalszych analiz użytecznym będzie rozłożyć pole promieniowania w wielomian w zbiorze zupełnym funkcji Uλ(x ), numerowanych przez pewien indeks λ. Spełniają one równanie Helmholtza :
( ∆ + ωλ2/c2 ) Uλ(x ) = 0 (8.5)
z odpowiednimi warunkami granicznymi i warunkiem normalizacji, takim że :
Wielkość ωλ oznacza częstość modu λ.
W dalszej kolejności nałożymy taki warunek iż mody są poprzeczne :
∇ • Uλ(x) = 0 (8.7)
Relacja zupełności dla modów przyjmuje postać :
gdzie δTij(x, x’ ) – jest poprzeczną δ -funkcją, która przy zastosowaniu do dowolnego pola wektorowego rzutuje go na składowa poprzeczną pola.
Przy wykorzystaniu funkcji modowych możemy przedstawić potencjał wektorowy następująco :
jednocześnie operator pola elektrycznego przyjmuje postać :
Operatory bλ†, bλ są operatorami, odpowiednio kreacji i anihilacji dla fotonów modu λ. Operatory te spełniają następujące relacje komutacji :
W szczególności przez :
| λ > = bλ† | 0 > (8.13)
oznaczymy stany jednofotonowe; | 0 > - jest próżniowym stanem pola.
Wykorzystując takie relacje komutacyjne, otrzymujemy jednoczasowe relacje komutacyjne dla operatorów pola np. :
Rzutowanie na składową poprzeczna, wbudowane w poprzeczną δ -funkcje z prawej strony powyższego równania gwarantują, że relacje komutacyjne są zgodne z prawem Gaussa ∇ • E = 0 dla pola swobodnego. W wyniku tego w cechowaniu radiacyjnym swobodny hamiltonian dla poprzecznych stopni swobody pola może być zapisany następująco :
Stopnie swobody pola EM związane są z elektronowymi stopniami swobody układu kwantowego. Przykładowo ze stopniami swobody atomu, którego hamiltonian oznaczymy jako HS. Hamiltonian HI opisujący oddziaływanie pomiędzy elektronami układu i polem promieniowania, może być zapisany następująco :
Rozpatrujemy tylko istotne stopnie swobody w przybliżeniu nierelatywistycznym i zaniedbujemy człon diamagnetyczny, który jest proporcjonalny do kwadratu potencjału wektorowego. Wielkość j(x, t) jest paramagentyczna elektronową gęstością prądu :
