Modelo de Ising. Introducción al Método Monte Carlo.

Introducción Teórica:

En el contexto de la Mecánica Estadística de equilibrio se deriva que un sistema en equilibrio termodinámico con (T,m,V) fijos, la probabilidad de encontrar al sistema teniendo N partículas y en un estado de las mismas X = (x(1),..,x(N)) es
PN(X) = Z-1e-bEN(X)-bmN     Z = ¥
å
 N = 0 
 ó
õ 		dX e-bEN(X)-bmN
(1)
donde T es la temperatura, b = 1/(kB T), kB es la constante de Boltzmann, m es el potencial químico, V es el volúmen, Z es la función de partición y EN(X) es la energía correspondiente al estado X de las N partículas. Así, conocida EN(X), formalmente podemos calcular el valor esperado de cualquier variable (que coincidirá con el valor observado en un experimento):
áAñ = Z-1 ¥
å
 N = 0 
 ó
õ 		dX e-bEN(X)-bmN AN(X)
(2)
Asi pues, la Mecánica Estadística conecta las ecuaciones del movimiento microscópicas de un sistema de partículas con sus propiedades termodinámicas macroscópicas. Sin embargo, en la práctica, sólo en unos pocos casos de interés, podemos realizar de forma analítica alguna de esas integrales N dimensionales. De esta forma debemos de acudir a métodos numéricos para conocer el comportamiento del sistema.

La primera estrategia que se nos ocurre es la de aplicar la ley de los grandes números. Ésta nos dice que independientemente de la distribución de probabilidad PN(X), el valor esperado de cualquier observable puede ser hallado mediante sumas parciales de M realizaciones concretas de los mismos, esto es

áA ñ = 1
M
 M
å
 i = 1 
 A(Xi)+O(M-1/2)
(3)
El gran problema de esta estrategia es que los Xi deben de seguir las distribución PN(X) y han de ser estadísticamente independientes. Para generar esas configuraciones podemos utilizar la técnica del rechazo. Ésta, se basa en el hecho obvio de que si f(x) es una cierta densidad de probabilidad (de una sola variable por sencillez), entonces si una y cumple que f(y) = 2f(x) entonces y es dos veces más probable que x. Asi el esquema para generar números aleatorios que sigan una cierta distribución f(x) hemos de seguir el siguiente algoritmo:

Si intentamos utilizar este sencillo algoritmo para generar configuraciones que sigan la distribución PN(X), nos pasaríamos todo el tiempo rechazando configuraciones. Esto es debido a que las configuraciones aleatorias son, normalmente, mucho más numerosas que las demás, por lo tanto, la probabilidad de escoger al azar una configuración aleatoria es prácticamente uno. Además, las configuraciones aleatorias son las que tienen más energía, puesto que la energía E(X) es proporcional al número de partículas, la probabilidad PN(X) en esos casos es cero (numéricamente hablando) y ningún número aleatorio uniforme generado por nosotros va a ser menor o igual a cero.

Para solucionar este problema Metropolis, Rosenbluth, Rosenbluth, Teller y Teller desarrollaron en 1953 un algoritmo que abandona la idea original de generar configuraciones estadísticamente independientes. En este caso, las configuraciones son construidas a través de cadenas de Markov. Una cadena de Markov de eventos sucesivos X1, X2, ..., XN tiene como probabilidad:

PN(X1, X2, ..., XN) = P1(X1)T(X1® X2)T(X2® X3).... T(XN-1® XN)
(4)
donde T(X® Y) es la probabilidad de transición de estando en X ir a Y. Obviamente

å
 Y 
 T(X® Y) = 1
(5)
Si g(X,t) es la probabilidad de que en el paso t de la cadena de Markov nos encontremos en el estado X podemos escribir fácilmente su ecuación de evolución
g(X,t+1) =
å
 X¢ 
 g(X¢,t)T(X¢® X)
(6)
    =
å
 X¢ ¹ X 
 g(X¢,t)T(X¢® X)+g(X,t)T(X® X)
(7)
    = g(X,t)+
å
 X¢ ¹ X 
 g(X¢,t)T(X¢® X)-
å
 X¢ ¹ X 
 g(X,t)T(X® X¢)
(8)
que no es más que la llamada ecuación Maestra. Una solución estacionaria de este proceso de Markov la encontramos pidiendo que
g(X¢)T(X¢® X) = g(X)T(X® X¢)
(9)
donde g(X) es la probabilidad estacionaria de que el sistema se encuentre en el estado X. Esta condición suficiente se denomina condición de balance detallado. Luego ya vemos la estrategia a seguir para generar estados X cuya distribución estacionaria sea PN(X) supuestamente conocida.

Esta dinámica de estados que nos hemos inventado recorre todos los estados accesibles por el sistema y, para tiempos largos, tiende a muestrear los estados con probabilidad PN(X). Obviamente, si elegimos al azar X¢ tendríamos el mismo problema que nos encontrábamos con el algoritmo de rechazo. El truco aquí es que podemos elegir un X¢ que no sea más que una pequeña modificación de X. De esta forma si en algún momento de la cadena de Markov las X son las más probables según PN(X) los siguientes estados generados seguirán estando concentrados alrededor de ese máximo de probabilidad.

Problemas:

Modelo y método numérico:

Lenz en 1927 propuso a su estudiante Ising modelar de forma sencilla unos materiales que sufrían un cambio brusco en sus propiedades magnéticas con la temperatura. A temperaturas altas eran paramagnéticos (magnetización espontánea cero) y a temperaturas suficientemente bajas eran ferromagnéticos (magnetización espontánea no nula). El modelo (hoy llamado de Ising) se define en una red cuadrada bi-dimensional N2 donde en cada nudo de la red existe una variable llamada de espin s(i,j)     i,j = 1,2,..,N que puede tener dos valores 1 y -1. La energía de interacción del sistema es:

E(S) = - 1
2
 N
å
 i = 1 
 N
å
 j = 1 
 s(i,j)[s(i,j+1)+s(i,j-1)+s(i+1,j)+s(i-1,j)]
(13)
donde suponemos que el sistema tiene condiciones de contorno periódicas i.e. s(0,j) = s(N,j), s(N+1,j) = s(1,j), s(i,0) = s(i,N) y s(i,N+1) = s(i,1). Dada una temperatura T, cuando el sistema está en equilibrio termodinámico la probabilidad de encontrar una configuración S viene dada por
P(S) = 1
Z
 e-bE(S)
(14)
Así, notar que cuando T = 0, P(S) = [1/2]Õi,jd(s(i,j)-1)+[1/2]Õi,jd(s(i,j)+1). Esto es, todos los espines tienen valor 1 o todos tienen valor -1. Por otra parte, si T = ¥ entonces P(S) = cte o sea, todas las configuraciones son igual de probables, pero como las configuraciones aleatorias son mucho más numerosas que las demás en el límite termodinámico (N®¥) sólo vemos éstas. Para generar configuraciones típicas con probabilidad de equilibrio P(S) utilizaremos el algoritmo de Metropolis.

Notar que este algoritmo cumple la premisa anteriormente contemplada, la diferencia entre una configuración S y la siguiente después de un paso de Markov S¢ es de un solo espín, s(m,n), todos los demás son iguales. De esta forma, la unidad de tiempo básica es el paso Monte Carlo que equivale a realizar N2 intentos de cambio de un espin. Asi, en un paso Monte Carlo, en promedio todos los espines han intentado cambiar una vez.