13.2: El modelo de Ising

13.2.1 Declaración del problema

Para comprender mejor la formulación general anterior del método MCMC, apliquémoslo al modelo 2D Ising, un modelo sencillo e instructivo que se utiliza comúnmente para enseñar conceptos de mecánica estadística. El sistema se describe mediante un conjunto de N “giros”, dispuestos en una celosía cuadrada 2D, donde el valor de cada giro\(S_{n}\) es\(+1\) (spin up) o\(-1\) (spin down). Esto describe un hipotético material magnético bidimensional, donde la magnetización de cada átomo está restringida a apuntar hacia arriba o hacia abajo.

Cada estado puede ser descrito por una cuadrícula de\(+1/=1\) valores. Por ejemplo, para una\(4\times 4\) cuadrícula, un estado típico se puede representar como

\[\begin{align}\{S\} = \left\{\begin{aligned}&+1\: +1\: +1\: -1 \\ &+1\: -1\: -1\: +1 \\ &-1\: +1\: +1\: -1 \\ &-1\: +1 \: -1 \: -1\end{aligned}\right\},\end{align}\]

y el número total de estados posibles es\(2^{16} = 65536\).

La energía de cada estado viene dada por

\[E(\{S\}) = - J \sum_{\langle ij\rangle} S_i S_j,\]

donde\(\langle ij\rangle\) denota pares de espines, en sitios adyacentes etiquetados\(i\) y\(j\), que son adyacentes entre sí en la cuadrícula (sin contar doble). Asumiremos condiciones de contorno periódicas en los bordes de la celosía. Así, por ejemplo,

\[\begin{align}\left\{\begin{aligned}&+1\: +1\: +1\: +1 \\ &+1\: +1\: +1\: +1 \\ &+1 \: +1\: +1\: +1 \\ &+1\: +1\: +1 \: +1\end{aligned}\right\} \quad E = -32J.\end{align}\]

\[\begin{align}\left\{\begin{aligned}&+1\: +1\: -1\: +1 \\ &+1\: +1\: +1\: +1 \\ &+1\: +1\: +1\: +1 \\ &+1\: +1\: +1\: +1\end{aligned}\right\} \quad E = -24J.\end{align}\]

Para cada estado, podemos calcular varias cantidades de interés, como el giro medio

\[S_{\mathrm{avg}}(\{S\}) = \frac{1}{N}\, \sum_i S_i.\]

Aquí,\(\mathrm{avg}\) denota el promedio sobre la celosía, para una configuración de giro dada. Entonces nos interesa el promedio termodinámico\(\langle S_{\mathrm{avg}}\rangle\), que se obtiene promediando\(S_{\mathrm{avg}}\) sobre un conjunto termodinámico de configuraciones de espín:

\[\left\langle S_{\mathrm{avg}}\right\rangle \;\; = \sum_{\mathrm{possible}\;\mathrm{states} \; \{S\}} S_{\mathrm{avg}}(\{S\}) \; \pi(\{S\}),\]

donde\(\pi(\{S\})\) denota la probabilidad de una configuración de giro:

\[\pi(\{S\}) = \frac{1}{Z}\, \exp\left(- \frac{E(\{S\})}{kT}\right).\]

13.2.2 Simulación Metropolis Montecarlo

Para aplicar el método MCMC, diseñamos un proceso de Markov utilizando el algoritmo Metrópolis discutido anteriormente. En el contexto del modelo de Ising, los pasos son los siguientes:

1. En el paso\(k\), elige aleatoriamente uno de los giros\(i\),, y considera un movimiento candidato que consiste en voltear ese giro:\(S_i \rightarrow -S_i\).
2. Calcular el cambio de energía que resultaría de voltear giro\(i\), relativo a\(kT\), es decir, la cantidad:
   \[Q \equiv \frac{\Delta E}{kT} = - \left[\frac{J}{kT} \sum_{j\;\mathrm{next}\;\mathrm{to}\;i} S_j\right] \Delta S_i,\]
   donde\(\Delta S_i\) está el cambio en\(S_{i}\) debido al giro, que es\(-2\) si\(S_{i}=1\) actualmente, y\(+2\) si \(S_{i}=-1\)actualmente. (La razón por la que calculamos\(Q \equiv \Delta E/kT\)\(\Delta E\), más que, es para mantener las cantidades en nuestro programa sin dimensiones, y para evitar tratar con números de punto flotante muy grandes o muy pequeños. Tenga en cuenta también que podemos hacer este cálculo sin sumar sobre toda la celosía; solo necesitamos encontrar los valores de los giros adyacentes al giro que estamos considerando voltear.)
   • Si\(Q \le 0\), acepta el spin-flip.
   • Si\(Q>0\), acepta el spin-flip con probabilidad\(\exp(-Q)\). De lo contrario, rechace el flip.
3. Esto nos dice el estado al paso\(k+1\) de la cadena de Markov (si el giro se volteó, o se quedó como estaba). Usa esto para actualizar nuestra media móvil de\(S_{\mathrm{avg}}\) (o cualquier otra media que nos interese).
4. Repetir.

El método MCMC consiste en aplicar repetidamente el proceso de Markov anterior, partiendo de algún estado inicial. Podemos elegir ya sea un estado inicial “perfectamente ordenado”, donde\(S_i = +1\) para todos los giros, o un estado “perfectamente desordenado”, donde cada uno\(S_{i}\) se asigna ya sea\(+1\) o\(-1\) aleatoriamente.

En algunos sistemas, la elección del estado inicial es relativamente poco importante; puedes elegir lo que quieras, y dejarlo a la cadena de Markov para que alcance la distribución estacionaria. Para el modelo de Ising, sin embargo, existe una razón práctica para preferir un estado inicial “perfectamente ordenado”, por la siguiente razón. Dependiendo del valor de\(J/kT\), el modelo de Ising se asienta en una fase “ferromagnética” donde los giros están mayormente alineados, o una fase “paramagnética” donde los giros son mayormente aleatorios. Si el modelo está en fase paramagnética y comienzas con un estado inicial ordenado (ferromagnético), es fácil que la celosía de giro se “derrita” en estados desordenados volteando giros individuales, como se muestra en la Figura\(\PageIndex{1}\):

Figura\(\PageIndex{1}\): Progreso de una simulación Montecarlo de un modelo de Ising con un estado inicial ordenado, para una\(30\times 30\) celosía con\(J/kT = 0.25\). La red de espín ordenada “se funde” en una configuración desordenada, que es el equilibrio termodinámico para este valor de\(J/kT = 0.25\).

En la fase ferromagnética, sin embargo, si comienzas con un estado inicial desordenado, la celosía de giro se “congelará” alineando los giros adyacentes. Cuando esto sucede, se formarán dominios grandes con giros opuestos, como se muestra en la Figura\(\PageIndex{2}\). Estos dominios separados no se pueden alinear fácilmente volteando giros individuales, y como resultado la cadena de Markov queda “atrapada” en esta parte del espacio estatal durante mucho tiempo, al no poder acceder al conjunto de estados más energéticamente favorable donde la mayoría de los espines forman un solo dominio alineado. (La simulación eventualmente se despegará, pero solo si esperas mucho tiempo). La presencia de dominios sesgará el cálculo de\(S_{\mathrm{avg}}\), porque los giros en diferentes dominios se cancelarán. De ahí que en esta situación es mejor iniciar la simulación MCMC en un estado ordenado.

Figura\(\PageIndex{2}\): Progreso de una simulación Monte Carlo de un modelo de Ising con una configuración inicial desordenada, para una\(30\times 30\) celosía con\(J/kT=1\). A medida que la red de espín desordenada “se congela”, forma dominios de larga duración que pueden interferir con los cálculos de\(S_{\mathrm{avg}}\).