Lars Onsager i model Isinga, czyli fizyka statystyczna a przejścia fazowe

Jesienią 1945 roku także uczeni wracali do pokojowego życia. Hendrik Casimir, doktorant Ehrenfesta i asystent Wolfganga Pauliego w ETH w Zurychu, lata wojny spędził w okupowanej Holandii, pragnął się więc dowiedzieć od swego dawnego szefa, co wydarzyło się w fizyce po stronie alianckiej: w Wielkiej Brytanii i w Stanach Zjednoczonych. Pauli, który spędził ten czas w Princeton, stwierdził, że w gruncie rzeczy niewiele się wydarzyło, prowadzono wprawdzie wiele prac nad radarem czy bombą atomową, ale w oczach Pauliego niezbyt się te kwestie liczyły. Dla niego ważne były dokonania intelektualne, a nie techniczne zastosowania. Właśnie jako dokonanie tego rodzaju – „arcydzieło analizy matematycznej” wyróżnił Pauli pracę Larsa Onsagera nad modelem Isinga z roku 1941. Na pochwałę ze strony Pauliego wyjątkowo trudno było zasłużyć, słynął on z ostrych ocen wygłaszanych często wprost w oczy („to nawet nie jest źle”). Był też wirtuozem trudnych technik, to on pierwszy rozwiązał problem atomu wodoru w mechanice kwantowej, w jej wersji macierzowej, zanim jeszcze powstało równanie Schrödingera.

Norweg pracujący w Stanach Zjednoczonych, Lars Onsager należał do wielkich dziwaków nauki. Karierę zaczął od tego, że zgłosił się do Petera Debye’a w ETH, by mu powiedzieć, że jego teoria elektrolitów jest błędna. Szybko przeniósł się za ocean. Studenci nazywali prowadzony przez niego przedmiot „sadistical mechanics” – wykłady były trudne, matematyczne, wykładowca mówił z norweskim akcentem, a do tego zasłaniał swą dużą sylwetką tablicę. W Yale dopiero po zaoferowaniu mu posady postdoca zorientowano się, że Onsager, mimo dorobku naukowego wciąż nie ma doktoratu. Napisał więc doktorat o funkcjach Mathieu, z którym wydział chemii nie wiedział, co zrobić. W tej sytuacji matematycy zaproponowali, że mogą tę pracę uznać za doktorat na ich wydziale. Ostatecznie przyznano mu doktorat z chemii. Onsager w latach czterdziestych wykazał, że dwuwymiarowy model Isinga wykazuje przejście fazowe. Całości bardzo długiej pracy nigdy zresztą nie opublikował, lubił podsycać zainteresowanie kolegów na konferencjach, pisząc np. na tablicy postać uzyskanego przez siebie ścisłego wyniku. Konkurowali pod tym względem z Feynmanem, który też lubił nagle wtrącić w dyskusji jakiś niepublikowany dotąd wynik. Przez pewien czas obaj zajmowali się nadciekłością helu i nabrali do siebie wzajemnego respektu.

Przez ostatnie kilkadziesiąt lat podano wiele rozwiązań problemu Isinga, jednak choć krótsze niż oryginalna praca Onsagera, nadal wymagają one sporo pracy i dość zaawansowanych technik, toteż ograniczymy się poniżej do zarysowania kontekstu, w którym ta praca się pojawiła.

Model Isinga to wyprany z wszelkich zbędnych szczegółów model ferromagnetyka, czyli materiału takiego jak np. żelazo, wykazującego namagnesowanie. Każdy atom stanowi dla nas strzałkę, która może być skierowana do góry albo na dół, czyli przeciwnie do wektora pola magnetycznego \vec{B} albo zgodnie z nim. Nasze strzałki są skrajnie uproszczoną wersją igły kompasu: mogą mieć tylko dwa zwroty. Gdy strzałka skierowana jest zgodnie z polem, ma niższą energię, gdy przeciwnie – wyższą.

Energie równe są odpowiednio \pm \mu B, gdzie \mu jest tzw. momentem magnetycznym (np. elektron ma ściśle określony moment magnetyczny). Na razie mamy do czynienia z paramagnetykiem, bo nasze strzałki zwracają się chętniej równolegle do wektora pola niż antyrównolegle. Gdy jednak pole wyłączymy, prawdopodobieństwa obu orientacji staną się równe.

Model Isinga opisuje styuację, gdy rozmieszczone w sieci krystalicznej spiny-strzałki położone najbliżej siebie wolą ustawiać się zgodnie. Równoległe ustawienie najbliższych sąsiadów ma energię -J, antyrównoległe J. Zauważmy, że teraz do energii dają wkład wszystkie pary najbliższych sąsiadów, czyli całkowita energia będzie sumą po linkach między sąsiadami (linki te zaznaczone są na czerwono). W przypadku dwywymiarowym zaznaczyliśmy energie dla tylko jednego spinu i jego sąsiadów, żeby nie zaśmiecać rysunku.

Ponieważ sąsiednie spiny chętnie ustawiają się równolegle, mamy w takim układzie do czynienia z bliskim porządkiem: nasi sąsiedzi mają te same poglądy co my, a przynajmniej korzystniejsze energetycznie jest, żeby mieli takie same poglądy. Pytanie podstawowe dla takiego układu brzmi: w jakich sytuacjach ten bliski porządek rozciągnie się na całą wielką sieć, dając zgodne uporządkowanie większości spinów – daleki porządek („prawie wszyscy mają takie same poglądy”). Mówimy tu o poglądach, bo model Isinga można stosować do opisu każdej sytuacji, gdy bliski porządek może wytworzyć porządek daleki. Stosuje się pewne warianty modelu Isinga do badania rozpowszechniania się plotek albo aktywności neuronów w mózgu. Rzecz więc nie musi dotyczyć tylko naszych strzałek-spinów i fizyki. My ograniczymy się tutaj do fizyki, ale warto sobie zdawać sprawę, że wiele zjawisk zbiorowych, kolektywnych można opisywać metodami fizyki.

Wracając do modelu Isinga: jego zachowanie będzie zależeć od temperatury, a ściślej mówiąc od porównania dwóch charakterystycznych energii: energii oddziaływania J z energią termiczną kT, gdzie k to stała Boltzmanna (inaczej mówiąc kT to temperatura wyrażona nie w stopniach, lecz w jednostkach energii). W niskich temperaturach dominować powinno uporządkowanie, w wysokich nieuporządkowanie. Gdzieś pomiędzy tymi dwoma obszarami następuje przejście fazowe ferromagnetyk-paramagnetyk (ferromagnetyk jest uporządkowany, ferrum to żelazo). Na symulacjach komputerowych sieci 400×400 atomów wygląda to tak.

kT=2,0JkT=2,27J

konfiguracja całkiem chaotyczna, bez bliskiego porządku

kT=2,5J

(Obrazki uzyskane za pomocą programu Dana Schroedera)

Przed drugą wojną światową nie można było oczywiście zrobić takiej symulacji komputerowej. Poza tym istotne jest udowodnienie, czy rzeczywiście model Isinga wykazuje przejście fazowe, a jeśli tak to w jakiej temperaturze, co dzieje się w jej pobliżu itp. itd.

Zacznijmy od spinów nieoddziałujących, czyli pierwszego obrazka u góry. Podstawowe prawo fizyki statystycznej mówi, że prawdopodobieństwo danego stanu układu zależy od energii tego stanu:

p=C\exp{\left( -\frac{E}{kT}\right)},

gdzie C jest stałą proporcjonalności. Jest to rozkład Gibbsa albo Boltzmanna-Gibbsa, choć można by go też nazywać rozkładem Boltzmanna-Gibbsa-Einsteina, ponieważ Einstein, pracownik Urzędu Patentowego, rozwinął tę technikę w wolnych od pracy chwilach. Boltzmann był tu prekursorem, ale zajmował się wyłącznie przypadkiem gazu. Gibbs uogólnił jego podejście i opublikował o tym książkę w Stanach Zjednoczonych, Einstein poznał ją po kilku latach i nawet stwierdził, że gdyby znał ją wcześniej, nie ogłosiłby trzech swoich prac z lat 1902-1904.

Dla spinu w polu magnetycznym mamy tylko dwa przypadki:

p_{\pm}=C\exp{\left(\mp \frac{\mu B}{kT}\right)}\Rightarrow C=\dfrac{1}{Z},\, \mbox{gdzie }\, Z=\cosh \left({\frac{\mu B}{kT}}\right).

Średnia wartość spinu w kierunku pola równa jest

M=(+1)p_{+}+(-1)p_{-}= \mbox{tgh}\left(\frac{\mu B}{kT}\right).

Dla układu N spinów należy po prostu tę wartość przemnożyć przez liczbę spinów. Gdy wyrazimy pole w jednostkach \frac{kT}{\mu}, a wartość spinu jako ułamek wartości maksymalnej M_0, otrzymamy po prostu wykres tangensa hiperbolicznego.

Gdy nie ma pola magnetycznego B, wypadkowy kierunek spinu jest równy M=0. Przy niewielkich wartościach pola M (magnetyzacja) jest proporcjonalna do B. Przy dużych wartościach osiągamy nasycenie – praktycznie wszystkie spiny ułożone są wówczas w jednym kierunku. (Tak się składa, że dla prawdziwego elektronu w polu magnetycznym wynik jest ten sam, choć spin elektronu różni się technicznie od naszej strzałki. Ale to tylko nawiasem. Pozostajemy przy strzałkach).

Uwzględnienie oddziaływań między spinami bardzo komplikuje problem, gdyż nie możemy już traktować spinów jako niezależne statystycznie. Na symulacjach u góry widać, że w różnych temperaturach wyniki są odległe od całkiem przypadkowego ułożenia, mamy do czynienia z bliskim porządkiem. Rozkład Gibbsa daje nam wtedy prawdopodobieństwa z osobna dla każdej konfiguracji spinów – jest ich 2^{N}. W dodatku, żeby uzyskać wiarygodne wyniki musimy uwzględnić dużo spinów, w skończonych próbkach przejścia fazowe się rozmywają. Jeśli chcemy coś udowodnić, trzeba umieć obliczyć granicę przy N dążącym do nieskończoności (co było główną trudnością Onsagera przy rozwiązywaniu modelu 2D).

Prosty przybliżony sposób poradzenia sobie z uwzględnieniem oddziaływań podał Pierre Weiss. Nazywa to się dziś przybliżeniem pola molekularnego. Otóż orientacja sąsiadów wpływa na energię danego spinu poprzez wartości \pm J. Jeśli spin środkowy zwrócony jest ku górze, to energia oddziaływań z sąsiadami jest równa

E_{+}=-Js_{+}+Js_{-}=-J(s_{+}-s_{-}),

gdzie s_{\pm} to liczba sąsiadów z odpowiednią orientacją. Podobnie

E_{-}=Js_{+}-Js_{-}=J(s_{+}-s_{-}).

Zauważmy, że obie nasz spin środkowy ma takie energie, jakby był w zewnętrznym polu magnetycznym o wartości \mu B=J(s_{+}-s_{-}). Jak dotąd wszystko jest ściśle, ale też i nic nie obliczyliśmy. Krok decydujący i przybliżony polega teraz na uznaniu, że możemy po prawej stronie ostatnich wyrażeń wstawić wartości średnie. Wtedy nasz spin znajduje się niejako w uśrednionym polu zewnętrznym – im bardziej spolaryzowani sąsiedzi, tym większa presja energetyczna na ustawienie się tak jak i oni. Zatem oddziaływania mogą wywierać taki sam skutek jak zewnętrzne pole magnetyczne. Uśrednione wartości liczby sąsiadów każdej orientacji są równe sp_{\pm}, gdzie s jest całkowitą liczbą sąsiadów (dla łańcucha 1D s=2, dla sieci kwadratowej 2D s=4). Możemy teraz wykorzystać wynik dla nieoddziałujących spinów i otrzymać równanie, które zawiera M po obu stronach. Rozwiązując to równanie, dostaje się magnetyzację jako funkcję temperatury w tym przybliżeniu. Wygląda ona następująco (nie ma tu zewnętrznego pola magnetycznego, to, co obserwujemy jest wyłącznie skutkiem oddziaływania spinów):

Temperatura, przy której magnetyzacja spada do zera, to tzw. temperatura Curie (chodzi o doktorat Pierre’a Curie jeszcze przed ślubem z naszą rodaczką Marią Skłodowską). Oczywiście magnetyzacje dodatnie i ujemne są tak samo możliwe. Układ ochładzany poniżej T_{c} ma tutaj dwie możliwości: zależnie od tego, co przeważy, wartości będą dodatnie bądź ujemne. Temperatura Curie równa jest

kT_c=Js.

Opisane zachowanie jest całkiem rozsądne z eksperymentalnego punktu widzenia. Jednak ścisłe rozpatrzenie modelu Isinga dla przypadku łańcucha 1D przynosi niezbyt przyjemny wniosek: układ nie ma w ogóle fazy ferromagnetycznej. A więc w tym przypadku przybliżenie pola molekularnego zawodzi kompletnie. Wynik ten był treścią doktoratu Ernsta Isinga w roku 1924. Podał on też argumenty na rzecz braku uporządkowania dalekiego zasięgu (ferromagnetyzmu) także w przypadku 2D.

Następnym wydarzeniem w dziejach tego modelu był argument Rudolfa Peierlsa opublikowany w roku 1936. Peierls, wychowanek Sommerfelda i Heisenberga, asystent Pauliego w ETH, nie miał po roku 1933 czego szukać w swej ojczyźnie, stając się jeszcze jednym z wielkich uczonych wypchniętych z Niemiec nazistowskich na emigrację. Z czasem pracował on w programie Manhattan i brytyjskim Tube Alloys, otrzymał brytyjski tytuł szlachecki. Niemcy już nigdy nie odzyskały swoich uczonych i swojej pozycji naukowej sprzed wojny. Argument Peierlsa, choć nie do końca prawidłowy w jego sformułowaniu, dowodził, że w dostatecznie niskich temperaturach 2D model Isinga ma fazę ferromagnetyczną.

OPiszemy krótko argument Peierlsa w wersji Wipfa (Statistical Approach to Quantum Field Theory, 2013). Wybierzmy na początek wszystkie spiny do góry, jest to stan o najniższej energii. Stany o orientacji ujemnej będą tworzyły wyspy rozmaitej wielkości, które można zamknąć konturem. Zbiór takich zamkniętych konturów określa jednoznacznie konfigurację spinów. Kontury ważne są dlatego, że po ich obu stronach mamy spiny skierowane antyrównolegle, czyli utworzenie takiego kontury, ściany domenowej, wymaga energii 2Jn, gdzie n to długość konturu.

 

.

Można następnie pokazać, że prawdopodobieństwo utworzenia konturu o długości n jest nie większe niż \exp{\left(-\frac{2Jn}{kT} \right)}. Wynika to z rozkładu Gibbsa, po drodze robi się następującą sztuczkę: zmieniamy znaki wszystkich spinów wewnątrz konturu: sam kontur wówczas znika, natomiast pozostałe energie się nie zmieniają.

Następny krok to wybranie jakiegoś spinu nie leżącego na krawędzi. Chcemy oszacować prawdopodobieństwo, że nasz spin będzie ujemny. Musi on leżeć wewnątrz jakiejś ściany domenowej o pewnej długości n. Możliwe wartości n są parzyste, począwszy od n=4 (samotny spin ujemny). Oszacujmy liczbę konturów A(n) zawierających nasz spin i mających długość n.

 

 

W tym celu prowadzimy od naszego spinu półprostą w prawo (szary kolor na rysunku). Musi ona przecinać jakiś pionowy kontur w jednej z odległości: \frac{1}{2},\frac{3}{2},\ldots, \frac{n-3}{2}. Ostatnia z odległości odpowiada konturowi prostokątnemu o wysokości 1 i długości \frac{n-2}{2}. Mamy więc tutaj (n-1) możliwości. Startując z tego przecięcia i wykonując pętlę, mamy do zrobienia (n-1) kroków, a w każdym nie więcej niż trzy możliwości. Zatem

A(n)\le \frac{n-2}{2} \cdot 3^{n-1}.

Prawdopodobieństwo, że nasz wybrany spin jest ujemny jest więc mniejsze niż

\displaystyle \sum_{n=4}^{\infty}\frac{n-2}{2}3^{n-1} \exp{\left(-\frac{2Jn}{kT}\right)}\le \dfrac{y^2}{3(1-y)^2},

gdzie y=9\exp{(-\frac{2J}{kT})}. Łatwo sprawdzić, że prawa strona nierówności maleje z temperaturą, a więc dla dostatecznie niskiej temperatury prawdopodobieństwo może stać się mniejsze niż \frac{1}{2}. Dotyczy to wszystkich spinów oprócz brzegu. A więc w dostatecznie niskiej temperaturze większość spinów będzie zwrócona tak jak na brzegu, czyli do góry.

W przypadku 2D wystąpuje więc faza ferromagnetyczna wbrew wnioskom Isinga. Onsager potrafił obliczyć funkcję Z=\sum_{\sigma} \exp{-\frac{E_\sigma}{kT}} po wszystkich konfiguracjach \sigma całej sieci. W roku 1948 obliczył też magnetyzację jako funkcję temperatury w tym modelu i napisał wynik na tablicy na dwóch różnych konferencjach. Ma ona następujący kształt.

Mimo upływu lat nie można uzyskać ścisłego rozwiązania 2D szybko, wszystkie metody są dość techniczne. Nie udało się też otrzymać rozwiązania w obecności pola magnetycznego. Także przypadek 3D pozostaje nierozwiązany, i to nie dlatego że nikt nie próbował. Kenneth Wilson, laureat Nobla za zjawiska krytyczne (a więc takie jak w modelu Isinga), wspominał w swoim wykładzie noblowskim, że kiedy jako świeżo upieczony naukowiec zastanawiał się nad przedmiotem badań dla siebie, poszedł zapytać Murraya Gell-Manna i Richarda Feynmana, nad czym aktualnie pracują. Gell-Mann pokazał mu model Isinga i powiedział, że gdyby udało mu się uzyskać rozwiązanie dla przypadku 3D, byłoby miło. Feynman, jak to Feynman – odrzekł, że nic nie robi.