Monte Carlo Simülasyonu ve Atomistik Spin Modelleri

Günümüzde fiziğin çeşitli dallarında yapılan araştırmalar teorik, deneysel ve
bilgisayara dayalı (computational) olmak üzere üç kısıma ayrılabilir. Monte
Carlo simülasyon tekniği 1953’de geliştirilmiş olup [1] temelde manyetik sistemleri
karakterize eden spin Hamiltonyenlerinin modellenmesi olmak üzere çok geniş bir
uygulama yelpazesine sahiptir. Bunların arasında sürekli spin sistemleri,
akışkanlar ve polimerlerin simülasyonlarını saymak mümkündür. Bu çalışmada manyetik
sistemlerin modellenmesi üzerinde duracağız.

Önemli Konfigürasyonların Örneklenmesi (Importance Sampling)

İstatistiksel mekanikte genellikle üzerinde çalıştığımız sistemde iç enerji ve
manyetizasyon gibi termodinamik ve manyetik özelliklerin ortalama değerlerini
hesaplamaya sıklıkla ihtiyaç duyarız.

(1) $\begin{equation*} \langle A \rangle=\frac{\sum_{{s}}A\exp{(-\beta H)}}{\sum_{{s}}\exp{(-\beta H)}}, {s}: \mathrm{sistemin\ \ bulunabilecegi\ \ tum\ \ olasi\ \ durumlar} \end{equation*}$

$N$ tane örgü noktası içeren ve her örgü noktasında iki durumlu bir değişkenin
yer aldığı bir örgüde bu ortalama ifadesi $2^{N}$ tane farklı şekillenimin
katkısını barındırır. Tipik olarak $N\geq40$ için bu ortalamayı analitik olarak
hesaplamak neredeyse imkansız bir hal alır. Bu zorluğu aşmanın bir yolu,
başlangıçta rastgele bir

spin şekillenimi seti'' oluşturarak aranan kesin çözüme
en yakın çözümü belli fiziksel kurallar yardımıyla elde etmektir. Ancak bu
şekillenim setinde bazı şekillenimler bölüşüm fonksiyonuna yeterince büyük katkı
verirken bazı şekillenimlerden ise neredeyse hiç katkı gelmez. Ayrıca söz konusu
şekillenim setinde herhangi bir şekillenimin sahip olduğu enerji gereği (yani
Boltzman ağırlığı ), bölüşüm fonksiyonuna katkı verecek olan
şekillenim sayısı da çok azdır. Bu nedenle bölüşüm fonksiyonu hesaplanırken tüm
şekillenimleri hesaba katmak yerine  ağırlığına göre yeterince
büyük katkı veren şekillenimleri hesaba katmak işi kolaylaştırır.

Importance
sampling” adını verdiğimiz bu yöntemde ilk aşama bir \textbf{Markov zinciri}
oluşturmaktır. Sistemin $t$ anındaki şekillenimi bir şekilde $t-1$ anındaki
şekillenimden elde ediliyorsa, bu şekilde üretilen $n$ tane şekillenimin
oluşturduğu sete Markov zinciri adı verilir. Markov zincirlerinde ergodiklik
hipotezi gereğince, belli bir şekillenimden diğerine sonlu adımda ulaşılmalıdır.

Markov Zincirleri ve Master Denklemi

Bir sistem için tanımlanmış bir süreçte belli bir alt durumunun, sürecin tahmin
edilebilir eylemlerinden ve rastgele elemanlardan yola çıkılarak
kestirilebildiği, ancak davranışı deterministik olmayan süreçlere stokastik
süreçler denir. $t_{1}$ , $t_{2}$ , $t_{3},…$ kesikli zaman aralıklarında
$S_{1}$ , $S_{2}$ , $S_{3},…$ ile verilen sonlu sayıdaki olası durumlar seti ile
tanımlanmış stokastik bir süreç ele alalım ve sistemin $t$ zamanında bulunduğu
durumu $X_{t}$ ile gösterelim. $X_{t_{n}}=S_{i_{n}}$ koşullu olasılığına göre

(2) $\begin{equation*} P(X_{t_{n}}=S_{i_{n}}|X_{t_{n-1}}=S_{i_{n-1}},X_{t_{n-2}}=S_{i_{n-2}},…,X_{t_{ 1}}=S_{i_{1}}) \end{equation*}$

yazabiliriz. Denklem (2)’den sistemin bir önceki zaman adımındaki
$X_{t_{n-1}}$ durumunun $S_{i_{n-1}}$ şekillenimi ile verildiğini anlıyoruz. Bu
tür süreçlere Markov süreçleri diyoruz. Denklem (2) ile verilen koşullu
olasılık, sistemin tüm durumlarından bağımsız ise ve son durumun olasılığı
sadece bir önceki durumun olasılığına bağlı ise,

(3) $\begin{equation*} P=P(X_{t_{n}}=S_{i_{n}}|X_{t_{n-1}}=S_{i_{n-1}}) \end{equation*}$

yazabiliriz. Bu durumda denklem (3)’e karşılık gelen ${X_{t}}$
durumlarının dizisine Markov zinciri denir. Böylece, denklem (3) ile
verilen koşullu olasılık

(4) $\begin{equation*} W_{ij}=W(S_{i}\rightarrow S_{j})=P(X_{t_{n}}=S_{j}|X_{t_{n-1}}=S_{i}) \end{equation*}$

ile verilen, $i$ durumundan $j$ durumuna geçişin olasılığı olarak
tanımlanabilir. Denklem (4) de tüm geçiş olasılıkları gibi

(5) $\begin{equation*} W_{ij}\geq0,\quad \sum_{j}W_{ij}=1 \end{equation*}$

şartlarını sağlamalıdır. Böylece, sistemin $t_{n}$ zaman adımında $S_{j}$
durumunda bulunma olasılığı $P(X_{t_{n}}=S_{j})$ ‘yi

(6) $\begin{equation*} P(X_{t_{n}}=S_{j})=P(X_{t_{n}}=S_{j}|X_{t_{n-1}}=S_{i})P(X_{t_{n-1}}=S_{i})=W_{ ij}P(X_{t_{n-1}}=S_{i}) \end{equation*}$

şeklinde tanımlayabiliriz.

$P(X_{t_{n}}=S_{j})=P(S_{j},t)$ için denklem (6)’nın zamana göre
değişimi master denklemi olarak adlandırılır.

(7) $\begin{equation*} \frac{dP(S_{j},t)}{dt}=-\sum_{i}W_{ji}P(S_{j},t)+\sum_{i}W_{ij}P(S_{i},t) \end{equation*}$

Denklem (7)’deki olasılık ifadesi korunumludur. Ayrıca, $j$ durumundan
$i$ durumuna geçerken $j$ ‘de bulunma olasılığının azalması veya tam tersi
şekilde $i$ durumundan $j$ durumuna geçerken $j$ durumunun olasılığının artması
nedeniyle denklem (7)’de tanımlanan master denklemi aynı zamanda bir
süreklilik denklemidir. Basitçe, bu denklem $S_{j}$ durumundan,
$S_{j}\rightarrow S_{i_{1}}$ , $S_{j}\rightarrow S_{i_{2}}$ , $S_{j}\rightarrow S_{i_{3}}$ gibi olası geçişleri, $j$ durumundan uzaklaşmanın toplam olasılığı
$\sum_{i}W_{ij}P(S_{j},t)$ şeklinde ifade etmektedir.

Metropolis Algoritması

Klasik Metropolis yöntemine göre \cite{metropolis} bir $t$ zamanındaki şekillenimi $t-1$ zamanındaki şekillenimden türetmek
için, ilk ve son durumlar arasındaki enerjilerin farkına bağlı olarak yazılan
bir geçiş olasılığı kullanılır. Üretilen durumlar zaman sırası ile dizilmiş,
ancak deterministik olmayan bir dizi oluşturur. Böyle bir sistemin zamana bağlı
davranışı

(8) $\begin{equation*} \frac{\partial P_{n}(t)}{\partial t}=-\sum_{n\neq m}\left[ P_{n}(t) W_{n\rightarrow m}-P_{m}(t)W_{m\rightarrow n}\right] \end{equation*}$

ile verilen master denklemi kullanılarak belirlenebilir. Denklem (8)’de
$P_{n}(t)$ , sistemin $t$ zamanında $n$ durumunda bulunma olasılığı,
$W_{n\rightarrow m}$ ise $n\rightarrow m$ geçişi olasılığıdır. Denge durumunda
$\partial P_{n}(t)/\partial t=0$ olur ve böylece denklem (8)’in
sağındaki iki terim birbirine eşit olur. Bu durumda elde edilen eşitlik
ayrıntılı denge koşulu (detailed balance condition) olarak adlandırılır ve

(9) $\begin{equation*} P_{n}(t)W_{n\rightarrow m}=P_{m}(t)W_{m\rightarrow n} \end{equation*}$

ile verilir. Klasik bir sistemde $n.$ durumun ortaya çıkma olasılığı ise

(10) $\begin{equation*} P_{n}(t)=e^{-E_{n}/k_{B}T}/Z \end{equation*}$

şeklindedir. Paydadaki bölüşüm fonksiyonundan ötürü bu olasılığın değerini kesin
olarak belirlemek genellikle olanaksızdır. Bu zorluğu aşmak için bir Markov
zinciri oluşturulur. Başka bir deyişle, her bir yeni durum, doğrudan bir önceki
durumdan türetilir. $n.$ durum $m.$ durumdan türetilirse, bağıl olasılığı
belirlemek için denklem (10)’a göre eski ve yeni durumlar arasındaki
enerji farkını bilmek yeterli olur:

(11) $\begin{equation*} \triangle E=E_{n}-E_{m} \end{equation*}$

Denklem (9) ile verilen detaylı denge koşulunu sağlayan her geçiş
olasılığı kabul edilebilir. Bu olasılığın ilk tanımı Metropolis tarafından

(12) $\begin{eqnarray*} W_{n\rightarrow m}&=&\tau_{0}^{-1}\exp(-\triangle E/k_{B}T), \quad \triangle E>0\ &=&\tau_{0}^{-1}, \qquad\qquad\qquad\qquad \triangle E<0 \end{eqnarray*}$

şeklinde yapılmıştır. Burada $\tau_{0}$ niceliği bir spinin değerinin değişmesi
(spin-flip) için geçen süredir ve genellikle 1 seçilir.

Metropolis algoritmasının uygulama adımları aşağıdaki gibidir [2,3]

Başlangıç şekillenimi oluşturulur
$i.$ örgü noktası seçilir
$i.$ örgü noktasındaki değişkenin işaret değiştirmesi durumunda oluşan enerji farkı $\triangle E$ hesaplanır
rastgele bir $r$ sayısı üretilir, $0<r<1$
Eğer $r<\exp(-\triangle E/k_{B}T)$ ise yeni şekillenim kabul edilir
Bir sonraki örgü noktasına geçilir ve (3)’den itibaren prosedür tekrar edilir

Markov zincirinin her bir adımında termodinamik değişkenlerin ortalaması hesaplanır.
İlk birkaç iterasyonda sistem henüz dengeye ulaşmamış olabilir. Bu nedenle ortalama hesaplarken bu adımların ürettiği şekillenimler hesaba katılmaz.
Sistemi dengeye getirmek için gerekli zaman şu şekilde bulunur: Simülasyon, aynı parametreler için birkaç farklı şekillenimden başlatılarak her bir başlangıç şekillenimi için sistemin manyetizasyonunun ortak saturasyon değerine kadar olan zaman adımları ihmal edilir.
Genellikle $N=L\times L$ kadar adımda sistemi dengeye getirmek yeterli olur [3]

Denklem (9)’un tek çözümü Metropolis fonksiyonu değildir. Alternatif
olarak Glauber geçiş fonksiyonu [4]

(13) $\begin{equation*} W_{n\rightarrow m}=\tau_{0}^{-1}\left[1+\sigma_{i}\tanh(E_{i}/k_{B}T)\right] \end{equation*}$

şeklindedir. Yüksek sıcaklıklarda (paramanyetik fazda) Metropolis algoritması
her bir adımda seçilen spinin çevrilmesine izin verir ve sistem iki durum
arasında salınım yapar. Bu nedenle Metropolis algoritması yeterince yüksek
sıcaklıklarda ergodiklik özelliğini yitirirken Glauber dinamiği ise tüm sıcaklık
spektrumu için ergodik kalır.

Ising modeli: Manyetik ve Termodinamik Özelliklerin Hesaplanması

Hamiltonyen:

(14) $\begin{equation*} \mathcal{H}=-J\sum_{<ij>}s_{i}s_{j} \end{equation*}$

Manyetizasyon:

(15) $\begin{equation*} M=\frac{1}{N}\left\langle\sum_{i}s_{i}\right\rangle \end{equation*}$

İç enerji:

(16) $\begin{equation*} E=\left\langle \mathcal{H}\right\rangle \end{equation*}$

Duygunluk:

(17) $\begin{equation*} \chi=\frac{1}{k_{B}T}\left(\left\langle M^{2}\right\rangle-\left\langle M\right\rangle^{2}\right) \end{equation*}$

Isı kapasitesi:

(18) $\begin{equation*} C_{h}=\frac{1}{k_{B}T^{2}}\left(\left\langle E^{2}\right\rangle-\left\langle E\right\rangle^{2}\right) \end{equation*}$

Binder kümülantı:

(19) $\begin{equation*} V_{L}(L,T)=1-\frac{\left\langle M^{4}\right\rangle}{3\left\langle M^{2}\right\rangle^{2}} \end{equation*}$

Şekil 1 (a): Ferromanyetik $(J>0)$ düzende farklı sıcaklıklar için şekillenimler. $T=1.0$ , $T=2.0$ , $T=T_{c}$ ve $T=3.0$ .

Şekil 1 (b): Antiferromanyetik $(J<0)$ düzende farklı sıcaklıklar için
şekillenimler. $T=1.0$ , $T=2.0$ , $T=T_{c}$ ve $T=3.0$ . Kırmızı
kareler $\uparrow$ , mavi kareler $\downarrow$ spin yönelimlerini göstermektedir.

Sayısal Sonuçlar

Denklem (14) Hamiltonyeni ile tanımlanmış olan iki boyutlu Ising modelin
manyetik ve termodinamik özelliklerinin örgü parametresi $L$ ve sıcaklık ile
değişimi (15)-(19) eşitlikleri kullanılarak hesaplanabilir.

Kare örgü $(z=4)$ için sistemin çeşitli sıcaklıklardaki ferromanyetik $(J>0)$ ve antiferromanyetik
$(J<0)$ şekillenimleri Şekil-(??)’de
görülmektedir. Ferromanyetizma ve antiferromanyetizma, Pauli ilkesi
ve Coulomb etkileşmelerine dayanan olgulardır. $s_{1}$ ve $s_{2}$ gibi iki spin
ele alırsak, bu iki spin arasındaki değiş-tokuş etkileşmesi $\pm Js_{1}s_{2}$
olacaktır. Bu durumda $J>0$ için değiş-tokuş enerjisi $-Js_{1}s_{2}$ negatif
değer alır. Bu nedenle spinler birbirleri ile paralel olarak yönelmeyi tercih
ederler ve kritik sıcaklığın altında sistemin net manyetizasyonu sıfırdan farklı
olur. $J<0$ için ise $-Js_{1}s_{2}$ enerjisi pozitif olacağı için spinler,
minimum enerji koşulunu sağlamak adına birbirleri ile anti-paralel olarak
yönelirler ve durum antiferromanyetizma olarak isimlendirilir. Antiferromanyetik
sistemlerin taban durumda $(T=0)$ net manyetizasyonu sıfırdır.

Ferromanyetizma ve antiferromanyetizama için Şekil-1’de gösterilen örgü
şekillenimleri sabit sıcaklıktaki Markov zincirinin bir halkası
olarak düşünülebilir. Bu şekillenimlerden üretilen Markov zincirlerinin sabit
sıcaklıktaki evinimleri yani zaman serileri Şekil-2’de
görülmektedir [5].
Düşük sıcaklıklarda manyetizasyon saturasyon değeri civarındadır. Artan
sıcaklıkla birlikte sistemin manyetizasyonu saturasyon değerinden uzaklaşır ve
sistemin sıcaklığı kritik sıcaklık civarına ulaştığında $+1$ ve $-1$ değerleri
arasında büyük dalgalanmalar ortaya çıkar. Şekil-1’den de görüldüğü gibi kritik
sıcaklık civarında sonlu ve değişken ebatlarda irili ufaklı kümelenmeler oluşur.
Sistem sıcaklığı kritik sıcaklığın üstüne çıktığında ise artık net manyetizasyon
sıfır değeri etrafında salınım yapar. Bu nedenle ortalama manyetizasyon değeri
sıfır olur ve sistem paramanyetik faza geçmiş olur.

Şekil 2: Ferromanyetik sistem için farklı sıcaklıklarda Markov zincirleri

Bir Markov zinciri boyunca sabit sıcaklıkta elde edilen termodinamik ve manyetik
özelliklerin ortalamaları ise o sıcaklıktaki ortalama değerleri
verir. Şekil-3 ‘de kare örgü için manyetizasyon, iç enerji, ısı sığası ve manyetik duygunluğun sıcaklıkla
değişiminin farklı örgü parametreleri $(L)$ için davranışı görülmektedir. $L$
değeri büyüdükçe hesaba katılan spin sayısı artar ve sistemin davranışı
termodinamik limite yaklaşır. Benzer şekilde, $L$ arttıkça ısı sığası ve manyetik duygunluk kritik sıcaklıkta
ıraksama davranışı sergiler.

Şekil 3: Kare örgü için manyetizasyon, iç enerji, ısı sığası ve manyetik duygunluk eğrilerinin sıcaklıkla değişimi. Sol paneldeki kesikli çizgi kesin çözümü göstermektedir.

Kritik Üsteller

Manyetizasyon ve korelasyon uzunluğu kritik üstelleri:

(20) $\begin{equation*} M\sim|T-T_{c}|^{\beta} \end{equation*}$

(21) $\begin{equation*} \xi\sim|T-T_{c}|^{-\nu} \end{equation*}$

Bu iki denklem birleştirilirse

(22) $\begin{equation*} M\sim\xi^{-\beta/\nu} \end{equation*}$

elde edilir. Termodinamik limitte $T=T_{c}$ ‘de $\xi\rightarrow\infty$ olur.
Ancak üzerinde çalıştığımız sistem sonlu ( $L\times L$ kare örgü) olduğu için
kritik sıcaklık noktasında korelasyon uzunluğu $\xi\rightarrow L$ olur.
Dolayısıyla lineer boyutları $L_{1}$ ve $L_{2}$ olan kare örgüler için

(23) $\begin{eqnarray*}\nonumber M_{1}&=&C_{1}L_{1}^{-\beta/\nu}\ \nonumber M_{2}&=&C_{2}L_{2}^{-\beta/\nu} \end{eqnarray*}$

yazılabilir. Bu iki eşitlik oranlanırsa

(24) $\begin{equation*} \log\left(\frac{M_{1}}{M_{2}}\right)=\log\left(\frac{C_{1}}{C_{2}}\right)-\frac{ \beta}{\nu}\log\left(\frac{L_{1}}{L_{2}}\right) \end{equation*}$

olduğu kolayca görülebilir. $T=T_{c}$ ‘de $C_{1}=C_{2}$ olduğu için

(25) $\begin{equation*} f(M,L)=\frac{\log\left(M_{1}/M_{2}\right)}{\log\left(L_{1}/L_{2}\right)}=-\frac{ \beta}{\nu} \end{equation*}$

Denklem (25)’in sol tarafı farklı $L_{1}$ ve $L_{2}$ çiftleri için
çizilirse $T=T_{c}$ ‘de bu eğriler birbirini keser. Kesim noktasının $x$ eksenini
kestiği yer kritik sıcaklığa, $y$ eksenini kestiği yer ise kritik üstele
karşılık gelir. Kritik üstel hesaplamanın bir diğer yolu denklem (22)’de
$\xi\rightarrow L$ seçerek
$T=T_{c}$ ‘de manyetizasyonun $L$ örgü parametresi ile değişimine bakmaktır.
Kritik sıcaklığı yüksek hassasiyetle hesaplamanın bir yolu ise denklem
(19) ile verilen Binder kümülantının $(V_{L})$ sıcaklıkla değişimini
farklı $L$ değerleri için incelemektir. Şekil-4’de (sol panel)
kümülantın sıcaklıkla değişimi görülmektedir. Farklı $L$ değerleri için çizilen
eğrilerin kesim noktası sistemin kritik sıcaklığına karşılık gelir. Sistemin
kritik sıcaklığı belirlendikten sonra

(26) $\begin{eqnarray*} \nonumber M(T)&\sim&L^{-\beta/\nu}\ \nonumber\chi(T)&\sim&L^{\gamma/\nu}\ C(T)&\sim&L^{\alpha/\nu} \end{eqnarray*}$

ile verilen eşitliklerden yararlanılarak net manyetizasyon, manyetik duygunluk
ve ısı kapasitesi için $\beta$ , $\gamma$ ve $\alpha$ üstelleri de elde
edilebilir. Kritik üstellerin analitik değerleri, $d$ boyut olmak üzere

(27) $\begin{equation*} \alpha+2\beta+\gamma=2 \end{equation*}$

(28) $\begin{equation*} d\nu=2-\alpha \end{equation*}$

Şekil 4: Sol panel: denklem (19)’den hesaplanan Binder kümülantının
farklı $L$ değerleri için sıcaklıkla değişimi. Sağ panel: denklem (25)
ile verilen $f(M,L)$ fonksiyonunun sıcaklıkla değişimi. Eğrilerin birbirini
kestiği noktanın $x$ bileşeni sistemin kritik sıcaklığını, $y$ bileşeni ise
$\beta$ üstelini verir.
eşitlikleri ile birbirine bağlıdır. Denklem (26)’dan manyetizasyon ve
manyetik duygunluk için elde edilen sonuçlar Şekil-5 ve
Şekil-6’da gösterilmiştir. Diğer yandan, Şekil-7’den de
görüldüğü gibi, ısı sığasının $L$ ile değişimi logaritmik olmadığı için
$log(C)-log(L)$ düzlemindeki davranışı doğrusal değil, eğrisel bir yapıya
sahiptir. Bu durum $\alpha=0$ olmasından kaynaklanır [6].
Kritik üsteller için $\nu=1$ kabul edilerek elde edilen sonuçlar
Tablo-1’de verilmiştir.

Şekil 5: Net manyetizasyonun $T=T_{c}$ ‘de örgü parametresi $L$ ile değişimi. Sol
panel doğrusal, sağ panel ise logaritmik ölçekte çizilmiştir.

Şekil 6: Manyetik duygunluğun $T=T_{c}$ ‘de örgü parametresi $L$ ile değişimi.
Sol panel doğrusal, sağ panel ise logaritmik ölçekte çizilmiştir.

Şekil 7: Isı kapasitesinin $T=T_{c}$ ‘de örgü parametresi $L$ ile değişimi. Sol
panel doğrusal, sağ panel ise logaritmik ölçekte çizilmiştir.

Tablo 1: 2 boyutlu Ising modelde termal ve manyetik özelliklerin kritik
üstelleri. MC-1: denklem (26)’dan MC-2 ise denklem (25)’den elde
edilen sonuçları göstermektedir.

Kaynaklar

[1] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, E. Teller, J. Chem. Phys. 21 (1953) 1087.
[2] D. P. Landau, K. Binder, A Guide to Monte Carlo Simulations in Statistical Physics, SE (2005) p.68
[3] M. E. J. Newman, G. T. Barkema, Monte Carlo Methods in Statistical Physics, Clarendon Press (Oxford) (2001)
[4]R. J. Glauber, J. Math. Phys. 4 (1963) 294.
[5] N. J. Giordano, Computational Physics, (1997) p.217
[6]H. Gould, J. Tobochnik, An Introduction to Computer Simulation Methods, SE (1996) p.572.