Il metodo Monte Carlo è un insieme di algoritmi che permette di simulare un gran numero di sistemi di fisica statistica. In questo articolo ci concentreremo su un algoritmo specifico, quello di Metropolis, e vedremo prevalentemente il caso particolare del modello di Ising.
Cosa calcoliamo
Consideriamo un sistema di spin, ad esempio il modello di Ising. La funzione di partizione è data da:
$$Z = \sum_{\{s\}} e^{-\beta H(\{s\})}$$
dove $H$ è l’Hamiltoniana del modello e $\{s\}$ è una configurazione di spin, per cui stiamo sommando su tutte le configurazioni di spin. La media di un’osservabile $O$ è data da:
$$\langle O \rangle = \frac{1}{Z}\sum_{\{s\}} O(\{s\}) e^{-\beta H(\{s\})}\tag{1}$$
Alcune scelte comuni sono $O=H$ oppure $O = \sum_i s_i$, cioè l’energia e la magnetizzazione del sistema. In altre parole stiamo calcolando la media di $O$ rispetto alla distribuzione di probabilità data da
$$p(\{s\}) = \frac{1}{Z}e^{-\beta H(\{s\})}$$
Questa è una vera e propria distribuzione di probabilità, perché soddisfa
$$0 \leq p(\{s\}) \leq 1\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \sum_{\{s\}} p(\{s\}) = 1$$
Data $p$, la media di un’operatore è data da $\langle O \rangle = \sum_{\{s\}} O(\{s\}) p(\{s\})$. Chiaramente in linea di principio potremmo calcolare $p(\{s\})$ analiticamente e poi effettuare il calcolo della media di qualsiasi operatore. Tuttavia un calcolo del genere è possibile in pratica solo su reticoli piccoli: ad esempio nel modello di Ising su un reticolo con $N$ siti abbiamo due possibili valori per sito, e quindi in totale $2^N$ configurazioni. L’esponenziale cresce molto rapidamente e quindi è in pratica impossibile effettuare il calcolo analiticamente se non per reticoli molto piccoli. Questo è un grosso ostacolo in pratica, perché come sappiamo molte proprietà dei modelli, come ad esempio le transizioni di fase, avvengono solo nel limite termodinamico $N \to \infty$ e quindi ci serve avere reticoli grandi. I metodi Monte Carlo servono ad aggirare questo problema.
Algoritmi Monte Carlo
I vari algoritmi Monte Carlo si basano tutti sulla stessa idea di base per aggirare la crescita esponenziale delle configurazioni. Invece di calcolare la probabilità esplicitamente, l’algoritmo genera una serie di configurazioni di spin
$$\{s^{(1)}\}, \{s^{(2)}\}, \ldots, \{s^{(T)}\}\tag{2}$$
dove $\{s^{(i)}\}$ è la $i-$esima configurazione di spin generata dall’algoritmo, e in questo caso abbiamo generato $T$ configurazioni totali. Nel modello di Ising, una configurazione di spin è un assegnamento di uno tra $-1, +1$ ad ogni spin su ogni sito del reticolo.
Le configurazioni $(2)$ non sono prese a caso, ma sono scelte secondo la probabilità $p(\{s\})$. Ovvero le configurazioni con $p$ maggiore appaiono con maggiore probabilità, mentre quello con $p$ piccolo appaiono raramente. Se il numero di configurazioni $T$ è abbastanza grande, il campione di configurazioni $(2)$ assomiglierà molto alla distribuzione $p(\{s\})$ che cerchiamo. A questo punto la media di un’osservabile $O$ può essere calcolata tramite
$$\langle O \rangle \approx \frac{1}{T}\sum_{i=1}^T O(\{s^{(i)}\})\tag{3}$$
Poiché le configurazioni sono state campionate secondo la probabilità $p(\{s\})$ possiamo cioè effettuare direttamente la media: le configurazioni più probabili appariranno più spesso e quindi contribuiranno di più, cosicché se $T$ è grande e il campionamento è fatto bene la $(3)$ è una buona approssimazione della $(1)$.
La domanda ora diventa: come si fa il campionamento delle configurazioni? In particolare ci serve un metodo per produrre configurazioni con probabilità $p$, senza però dove calcolare $p$. Altrimenti, se dovessimo calcolare $p$ esplicitamente siamo daccapo, perché per calcolare $p$ dobbiamo calcolare $Z$, e per calcolare $Z$ dobbiamo sommare tutte le configurazioni.
Il punto fondamentale di un algoritmo Monte Carlo è che è possibile produrre una serie di configurazioni secondo una probabilità $p$ senza calcolare $p$. Questo fatto è abbastanza sorprendente, ma è possibile e vedremo come nel prossimo articolo.