Lattice QCD per principianti - 5/N

Digressione numerica - 2/2

Ah, già, eravamo rimasti alle catene di Markov. Nome che come potrete facilmente comprendere mi è molto simpatico. Dunque, abbiamo visto che la distribuzione d'equilibrio del nostro sistema statistico è la distribuzione di Boltzmann: $$ P_{\text{eq}}({\cal C}) = \frac{\text{e}^{-\beta H[{\cal C}]}}{Z} $$ Il nostro obiettivo è quello di generare configurazioni ${\cal C}$ che siano distribuite in accordo a $P_{\text{eq}}$.

Adesso immaginiamo di inizializzare il nostro sistema con una configurazione qualunque, scelta arbitrariamente (purché legale). Chiaramente questa configurazione non seguirà la distribuzione $P_{\text{eq}}$ ma una certa distribuzione arbitraria, quasi sicuramente ignota. Ma questo non ci importa. Adesso definiamo una regola stocastica che ci permetta di passare da una configurazione ${\cal C}$ a un'altra configurazione ${\cal C}'$. La probabilità di transizione $W({\cal C} \to {\cal C}')$ deve soddisfare due requisiti: il primo $$ W({\cal C} \to {\cal C}') \ge 0 $$ è abbastanza ovvio, e il secondo $$ \sum_{ {\cal C}'} W({\cal C} \to {\cal C}') = 1 $$ ci assicura solo che $W$ sia una probabilità di transizione ben normalizzata, che definisce un processo markoviano. Questo significa solo che la probabilità di estrarre la configurazione ${\cal C}'$ dipende solo dalla configurazione precedente, ${\cal C}$, e non da tutta la storia che abbiamo seguito per arrivare a ${\cal C}$.

Ora possiamo pensare di applicare $W$ non alle configurazioni, ma a $P$ stessa; possiamo farlo perché possiamo pensare di applicare $W$ a un intero ensemble di configurazioni distribuite secondo $P$. Potremo quindi scrivere qualcosa del genere: $$ W[P({\cal C}) \to P({\cal C}')] = \sum_{\cal C} W({\cal C} \to {\cal C}') P({\cal C}) $$ Immaginiamo che $W$ soddisfi queste due condizioni:

• $P_{\text{eq}}$ è un punto fisso della trasformazione, cioè $W$ preserva $P_{\text{eq}}$; in formule $$ \sum_{\cal C} W({\cal C} \to {\cal C}') P_{\text{eq}}({\cal C}) =  P_{\text{eq}}({\cal C}') $$
• $W$ è ergodica: questo significa che a partire da una determinata configurazione legale del sistema possiamo raggiungere qualsiasi altra configurazione legale del sistema in un numero finito di passi.

Ebbene, si può dimostrare che queste due semplici condizioni sono sufficienti a garantire che se partiamo da una $P$ qualunque, dopo un certo numero di applicazioni di $W$ finiremo per atterrare su $P_{\text{eq}}$, e una volta atterrati lì non ci muoveremo più. In altre parole dopo un certo periodo di così detta "termalizzazione", cominceremo a estrarre configurazioni distribuite secondo Boltzmann.

Lasciando da parte la condizione di ergodicità, che è comunque abbastanza semplice da implementare, si può dimostrare che la condizione di bilancio dettagliato, ossia $$  \frac{W({\cal C} \to {\cal C}')}{W({\cal C}' \to {\cal C})}  = \frac{P_{\text{eq}}({\cal C}')}{P_{\text{eq}}({\cal C})} = \text{e}^{-\beta \Delta H} \quad\quad \Delta H = H[{\cal C}'] - H[{\cal C}]$$ è sufficiente per far sì che $W$ soddisfi la prima condizione. Notate, incidentalmente, che nel rapporto la funzione di partizione si cancella: non ne avremo mai bisogno per implementare $W$.

L'algoritmo di Metropolis

Il modo più semplice e comunemente usato per implementare una $W$ che soddisfi il bilancio dettagliato è quello di usare l'algoritmo di Metropolis: $$ W({\cal C} \to {\cal C}') = \text{min}\{ 1, \exp(-\beta\Delta H)  \}$$ Vediamo perché questa scelta soddisfa il bilancio dettagliato:

1. Se passando da ${\cal C}$ a ${\cal C}'$ l'energia diminuisce o rimane uguale, allora $\exp(-\beta\Delta H) \ge 1$ e la nuova configurazione sarà sempre accettata: dunque $$  \frac{W({\cal C} \to {\cal C}')}{W({\cal C}' \to {\cal C})}  = \frac{1}{\text{e}^{\beta\Delta H}} = \text{e}^{-\beta\Delta H}$$
2. Se invece l'energia aumenta, allora avremo direttamente $$  \frac{W({\cal C} \to {\cal C}')}{W({\cal C}' \to {\cal C})}  = \text{e}^{-\beta\Delta H}$$ Fatto.

Vediamo ora come funziona in pratica, nel caso di Ising, l'algoritmo di Metropolis.

1. Scegliete una configurazione di spin arbitraria. In genere si parte o da una configurazione "fredda" (tutti gli spin allineati) o da una configurazione "calda" (il valore degli spin è scelto a caso con distribuzione piatta);
2. Selezionate un particolare spin, $s_x$, a caso tra tutti gli spin disponibili, e flippatelo: $s_x \to -s_x$; questa sarà la nuova configurazione di prova;
3. Calcolate $\Delta H$;
4. Estraete un numero random $r$ con distribuzione piatta, $0 \le r \lt 1$;
5. Se $r < \exp(-\beta\Delta H)$ accettate la nuova configurazione, altrimenti tenetevi la vecchia; notate che in questo secondo caso dovete considerare la vecchia configurazione come nuova a tutti gli effetti pratici, e quindi dovrà essere tenuta in conto per le medie statistiche;
6. Ripetete i passi da (2) a (5) un numero sufficiente di volte perché tutti gli spin del sistema siano visitati, in media, una volta; in pratica se avete $N$ spin ripeteteli $N$ volte; questo costituisce uno sweep  Monte Carlo;
7. Aspettate $N_t$ sweep di termalizzazione: servono per far arrivare la configurazione alla distribuzione d'equilibrio;
8. Fate altri $M$ sweep, e a ogni sweep calcolate le vostre osservabili, in modo da poter costruire le medie Monte Carlo, gli errori statistici eccetera. 

Se fate i bravi, una delle prossime volte fornirò un codice python per simulare il modello di Ising bidimensionale con l'algoritmo di Metropolis. Per il momento vi dovrete accontentare di un po' di disegnetti, che servono a spiegare il punto (7) dell'algoritmo: come possiamo sapere quando il sistema ha termalizzato? Nella prima figura mostro l'andamento dell'energia, in funzione del tempo Monte Carlo, cioè del numero di sweep, per un sistema con $32\times 32$ siti alla temperatura critica: la curva verde si riferisce a una partenza "fredda", quella fucsia a una partenza "calda". Come vedete dopo 100 sweep le curve "si sovrappongono", e possiamo pensare che il sistema abbia termalizzato, nel senso che ha "dimenticato" la situazione di non-equilibrio di partenza.

Naturalmente la situazione cambia drasticamente se aumentiamo le dimensioni del sistema: questa è la stessa figura di prima ma per un sistema $256\times 256$.

Come avrete notato la convergenza tra le due curve è molto più lenta. Questo è dovuto alle proprietà specifiche del sistema alla temperatura critica, il cui il modello di Ising in $d=2$ mostra una transizione di fase del secondo ordine e la lunghezza di correlazione va a infinito. Ma avremo modo di riparlarne.

Avvertenza per i più scaltri: l'algoritmo di Metropolis non è esattamente il più furbo per simulare Ising, ora come ora sono noti algoritmi molto più efficaci, ma questi sono dettagli tecnici, di cui forse parlerò, forse no.