Repositorio Institucional de Documentos

Trabajo Fin de M´ aster Desarrollo de un algoritmo de simulaci´on geom´etrico para modelos de gauge acoplados a un campo de Higgs Autor: Eduardo Royo Amondarain Directores: Vicente Azcoiti P´ erez Eduardo Follana Ad ´ ın M´aster en F´ısica y Tecnolog´ıas F´ısicas Zaragoza, Junio 2014 Facultad de Ciencias Resumen En este trabajo se plantea como objetivo la implementaci´on de un algoritmo de simulaci´on geom´etrico para el modelo Z2gauge-Higgs. Tras una exposici´on de los resultados relevantes en el campo, se introduce el modelo de estudio y se prueba su dualidad con el modelo de Ising con campo magn´etico en d= 2. Posteriormente se describe en detalle el funcionamiento de los algoritmos implementados y se presentan los resultados obtenidos por estos, incluyendo una comparativa de la eficiencia computacional de cada uno. Adicionalmente, se estudia el modelo de Ising en su sector antiferromagn´etico, corroborando la existencia de una l´ınea de transici´on que persiste para valores del campo no nulos, a diferencia del caso ferromagn´etico. ´ Indice 1. Introducci´on 1 2. Fundamentos te´oricos 3 2.1. El modelo Z2gauge-Higgs.................... 3 2.2. Caso d= 2: Dualidad con el modelo de Ising . . . . . . . . . . 6 2.3. Observables de inter´es: e, m, CV, χ ............... 8 3. M´etodos de Monte Carlo: Algoritmos de Metropolis 9 3.1. Aspectos generales . . . . . . . . . . . . . . . . . . . . . . . . 10 3.2. Algoritmos implementados . . . . . . . . . . . . . . . . . . . . 12 3.2.1. Algoritmo geom´etrico . . . . . . . . . . . . . . . . . . 13 3.2.2. Algoritmo PS . . . . . . . . . . . . . . . . . . . . . . . 15 3.2.3. Algoritmo h´ıbrido . . . . . . . . . . . . . . . . . . . . 18 4. Resultados num´ericos 21 4.1. Caso h=0 ............................ 21 4.2. Caso h6=0 ............................ 24 5. Conclusiones y perspectivas 29 1. Introducci´on Hace ya algunos a˜nos, los trabajos de Prokof’ev y Svistunov en modelos cu´anticos [1] y posteriormente en sistemas de espines cl´asicos (Ising, XY, Potts) en dos y tres dimensiones [2], introdujeron un nuevo m´etodo que permite abordar num´ericamente el estudio de modelos f´ısicos tanto en materia condensada como en altas energ´ıas, mediante diferentes generalizaciones que se han desarrollado en la ´ultima d´ecada [3–5]. 1 Dicho m´etodo consiste en general en reemplazar la funci´on de partici´on original, dada por una suma sobre configuraciones de los campos involucrados en el modelo particular, por una reformulaci´on de esta, usualmente un desarrollo (ya sea de acoplo fuerte, altas temperaturas o en potencias del par´ametro de hopping en lattice QCD). En un segundo paso, se extiende el nuevo espacio de configuraciones introduciendo defectos en la red, y haciendo que se propagen mediante un algoritmo de gusano, llamado as´ı por la forma de las cofiguraciones que genera. Esto permite extraer informaci´on directa de las funciones de correlaci´on a dos puntos. M´as all´a de que pueda mejorar la eficiencia en el c´alculo de algunos observables, este enfoque permite superar, al menos en algunos casos, dos problemas fundamentales en el ´area de la simulaci´on de sistemas f´ısicos. Por un lado ya en [2] y posteriormente en an´alisis m´as detallados de Deng y colaboradores en [6] y Wolff en [7], el algoritmo no presenta apenas critical slowing down (CSD), es decir, el exponente din´amico zcon el cual escala el tiempo de autocorrelaci´on integrado1es pr´oximo a cero. Como consecuencia, se hace posible el c´alculo de diferentes propiedades en los puntos cr´ıticos de los modelos, por ejemplo exponentes cr´ıticos, algo imposible para los algoritmos que presentan CSD. No menos importante es el llamado problema del signo, que afecta a los casos en los que la acci´on del sistema es compleja, lo que ocurre por ejemplo en QCD a densidad finita (µ > 0) o en sistemas con t´erminos topol´ogicos (θterm en QCD). Se trata de un problema computacional de la m´axima dificultad, NP-hard, lo que implica que una soluci´on general para este (en tiempo polin´omico) probar´ıa2P=NP, tal y como se explica en [9]. Aunque en principio es un asunto t´ecnico, se trata en la pr´actica de un obst´aculo insalvable en muchos casos. Si bien este algoritmo no soluciona el problema del signo de forma general, la reformulaci´on de la funci´on de partici´on llevada a cabo en trabajos posteriores, ha permitido evitarlo completamente en algunos modelos, recientemente en [4,5] con Z3gauge-Higgs, [10] con U(1) gauge-Higgs y [11] con el modelo de Ising con campo magn´etico imaginario. El objetivo de este trabajo es el desarrollo de un algoritmo de simulaci´on geom´etrico (en el sentido que se ver´a m´as adelante), en principio en un sistema suficientemente sencillo, que pueda servir para testear el algoritmo, y que a su vez tenga inter´es en altas energ´ıas. El modelo Z2gauge-Higgs es el elegido por cumplir lo anterior, ya que en d= 2 es dual al modelo de 1τint da una medida de la autocorrelaci´on de la serie de datos obtenida en una simulaci´on de Monte Carlo. Si el tama˜no total de la muestra es N, ´unicamente Neff =N/2τint datos son independientes, de cara a calcular los errores de los estimadores. Como τint ∝Lz, los valores habituales de zen los algoritmos locales (z≃2) hacen inviable simular tama˜nos del ret´ıculo grandes. Para un an´alisis m´as profundo de τint consultar la secci´on 3.1 ´o [8]. 2P=NP es uno de los problemas del milenio propuestos por el Clay Mathematics Institute, premiado con un mill´on de d´olares. 2 Ising con campo magn´etico en una red cuadrada. Dado que el algoritmo del gusano original [2] se realiza en su forma m´as sencilla precisamente en este ´ultimo modelo, pero sin campo magn´etico, parece un buen campo de pruebas para intentar generalizar el algoritmo. Adem´as, gracias a esta dualidad, podemos abordar tambi´en el modelo de Ising antiferromagn´etico, un sistema que pese a haber sido estudiado con anterioridad [12–14], incluso localizando mediante m´etodos de Monte Carlo [15,16] una l´ınea de transici´on3, no ha sido estudiado con detalle en las regiones cr´ıticas mediante m´etodos num´ericos. En resumen, buscamos desarrollar nuevos algoritmos en un modelo test, Z2gauge-Higgs, que nos permitan superar las dificultades habituales al aplicar m´etodos de Monte Carlo a modelos tanto de altas energ´ıas como de materia condensada. Por otro lado, dada la dualidad de este modelo con el modelo de Ising con campo magn´etico, aprovecharemos para estudiar el sector antiferromagn´etico de este ´ultimo. 2. Fundamentos te´oricos En esta secci´on describiremos el modelo que usaremos a lo largo de todo el trabajo, as´ı como la dualidad que presenta con el modelo de Ising con campo magn´etico en una red cuadrada. Presentamos tambi´en los observables de inter´es que analizaremos en las siguientes secciones mediante m´etodos de Monte Carlo. 2.1. El modelo Z2gauge-Higgs Partiremos del modelo Z2gauge-Higgs (tambi´en Ising gauge-Higgs) introducido por Wegner en [17] y Balian y colaboradores en [18]. En este modelo tenemos una red c´ubica, en principio de dimensi´on arbitraria dy tama˜no L, con condiciones de contorno peri´odicas. El n´umero de puntos N(sites), que indexaremos con n= 1, . . . , N, es N=Ld. Las variables fundamentales ser´an los links σnµ ∈Z2que los unen (µindica la direcci´on espacio-temporal) que se acoplar´an a los campos de materia zn∈Z2, que residen en los propios sites. Ser´a de utilidad definir tambi´en las llamadas variables plaqueta σnµν := σnµσn+µ,νσn+ν,µσn,ν,(1) como el producto de cuatro links contiguos, formando un cuadrado. Con estos elementos se construye la siguiente acci´on: S=−βGX n,µ,ν µ<ν σnµν −βIX n,µ znσnµzn+µ.(2) 3En esto difiere del modelo ferromagn´etico, que no tiene transici´on de fase para h6= 0. 3 Podemos simplificar la expresi´on trabajando con variables invariantes gauge (o lo que es lo mismo, escogiendo la gauge unitaria). Formalmente esto se traduce en introducir un cambio de variable wnµ := znσnµzn+µ.(3) Por construcci´on wnµν =σnµν (las plaquetas son invariantes gauge), y la acci´on queda igual que en (2), salvo por la ausencia de las variables zn. De este modo, tomamos directamente las σnµ como variables invariantes gauge y tenemos que Sse reduce a S=−βGX n,µ,ν µ<ν σnµν −βIX n,µ σnµ.(4) Para escribir la funci´on de partici´on Zasociada a esta acci´on, nos ayudamos del hecho de que las variables link (y tambi´en las variables plaqueta, por ser producto de las anteriores) est´an en Z2(es decir, son iguales a ±1), verific´andose entonces la siguiente identidad elemental, eβσ = cosh β+ sinh βσ = cosh β(1 + tanh βσ),(5) que nos permite escribir Z=Pconf exp(−S) de la siguiente forma Z= (cosh βG)d(d−1)N/2(cosh βI)dN × X {σnµ}Y n,µ,ν µ<ν (1 + tanh βGσnµν)Y n,µ (1 + tanh βIσnµ).(6) Para simplificar el trabajo conviene realizar un cambio de variable x= tanh βG, y= tanh βI,(7) y prescindir del factor constante que precede a la suma sobre configuraciones. De este modo, la funci´on de partici´on con la que trataremos viene dada por Z=X {σnµ}Y n,µ,ν µ<ν (1 + xσnµν)Y n,µ (1 + yσnµ).(8) El siguiente paso es clave: la funci´on de partici´on es una suma sobre todas las configuraciones posibles de links {σnµ}. Antes de sumar, tenemos un polinomio en xeyde la forma (1 + xσP1)···(1 + xσPNP)(1 + yσl1)···(1 + yσlNl),(9) donde Nl(NP) es el n´umero total de links (plaquetas). En este punto se produce cierto abuso de notaci´on, ya que tanto para links como para plaquetas 4 se emplea la letra σ, de ah´ı el sub´ındice loPpara diferenciar ambos casos. Recordar adem´as c´omo una plaqueta es el producto de cuatro links contiguos. Expandiendo el producto, tendr´ıamos un total de 2Nl+NPsumandos. Podemos escribir algunos de ellos para ejemplificar: 1. y2σl1σl2(Escogiendo los dos primeros links, y el resto 1’s) 2. xy σPiσlj(Escogiendo la plaqueta iy el link j) 3. xy4(Escogiendo una plaqueta y los cuatro links asociados) Sin embargo, al realizar la suma sobre configuraciones, dado que σ∈Z2, los sumandos que contengan cualquiera de las variables link un n´umero impar de veces se ver´an anulados, ya que X σi=±1 σi(···)=0.(10) De los tres ejemplos anteriores, ´unicamente el tercero (t´ermino proporcional a la identidad del grupo, que en este caso es Z2), contribuye a la funci´on de partici´on. En general, un sumando cualquiera no se anular´a si no contiene ninguna variable link suelta, o dicho de otro modo, si todas las variables link de dicho sumando aparecen un n´umero par de veces. En caso contrario, la restricci´on (10) evita que el t´ermino contribuya. Notar que cada link puede aparecer una vez por el t´ermino que va con y, y hasta 2(d−1) veces m´as, una por cada plaqueta que lo contiene. De lo anterior se deduce que cada sumando en (9) viene caracterizado por una lista de plaquetas y links, que denominaremos activos. De estos sumandos, ´unicamente contribuyen a la funci´on de partici´on (8) aquellos en los que cada link aparece un n´umero par de veces, y todos estos sumandos son simplemente multiplicados por un factor 2Nal hacer la suma sobre configuraciones, ya que son proporcionales a la identidad. Lo anterior nos permite reinterpretar cada sumando como una nueva configuraci´on gde la funci´on de partici´on (8). Cada una de estas configuraciones viene caracterizada en principio por una lista de plaquetas activas y links activos. Si ahora definimos la frontera de las plaquetas activas como el conjunto de links compartidos por un n´umero impar de plaquetas, la restricci´on de que cada link debe aparecer un n´umero par de veces se traduce en que una configuraci´on v´alida viene dada por un conjunto de plaquetas activas, en cuya frontera residen links activos. Por tanto cada configuraci´on gviene caracterizada ´univocamente por una lista de plaquetas activas, que denotaremos por {σnµν}, y su peso depende del n´umero de estas y del n´umero de links que contiene su frontera. Es posible dar una interpretaci´on m´as gr´afica a lo anterior, visualizando cada configuraci´on como (hiper)superficies, cerradas o abiertas, formadas por plaquetas activas. En las fronteras de las superficies abiertas es donde 5 encontrar´ıamos los links activos. En el sentido anterior introducimos dos magnitudes, ´area y per´ımetro, como A≡N´umero de plaquetas activas, P≡N´umero de links activos. Con todo lo expuesto anteriormente, podemos reformular la funci´on de partici´on original como una suma de grafos gadmitidos, cada uno con un peso w(g) = xA(g)yP(g), es decir Z=X g∈G w(g) = X g∈G xA(g)yP(g)(11) Siendo Gel conjunto de grafos permititdos, y donde hemos prescindido del factor com´un constante 2N. Destacar que existen ejemplos recientes en la literatura en los que se emplea este tipo de reformulaci´on de la funci´on de partici´on, ver por ejemplo [3] para el modelo U(1) gauge o [4] para Z3 gauge-Higgs, as´ı como [11] para el modelo de Ising con campo imaginario. Podemos preguntarnos ahora cu´antos grafos permitidos contribuyen a la funci´on de partici´on (11). Dado que una configuraci´on gviene caracterizada por una conjunto de plaquetas activas, habr´a tantas c´omo en {σnµν}, donde cada σnµν puede estar activa o no. Por tanto, tendremos 2NPgrafos g posibles. Es decir, hemos pasado de 2Nl= 2dN t´erminos en la funci´on de partici´on original a 2NP= 2d(d−1)N/2contribuciones en la funci´on reformulada. En cualquier caso, incluso para tama˜nos de la red modestos, contamos con un n´umero de grafos demasiado grande como para sumar todas sus contribuciones de forma num´erica; a modo de ejemplo, un ret´ıculo bidimensional 10 ×10 tendr´ıa 2100 sumandos. Por este motivo, recurriremos al empleo de m´etodos de Monte Carlo, explicados en la secci´on 3. 2.2. Caso d= 2: Dualidad con el modelo de Ising El modelo Z2gauge-Higgs es dual al modelo de Ising con campo magn´etico para d= 2 [18]. Es relativamente sencillo ver esto partiendo de los resultados anteriores. Hemos visto c´omo las configuraciones de grafos en (11) se correspond´ıan una a una con las configuraciones de plaquetas {σnµν}. Definimos entonces los espines sien el ret´ıculo dual (que es el ret´ıculo cuyos sites estan situados en el centro de las plaquetas del ret´ıculo original). As´ı pues, para cada i∈[1, N] establecemos la siguiente correspondencia, σi12 =activa ↔si=−1, σi12 =inactiva ↔si= +1.(12) Ahora podemos encontrar significado a las variables AyPempleadas en 11, y expresarlas en t´erminos de estos espines duales. El ´area Aes el n´umero de 6 plaquetas activas, es inmediato entonces que podemos expresar Acomo A=1 2 N X i=1 (1 −si).(13) Por otra parte, la frontera de cada configuraci´on viene definida en d= 2 por el conjunto de Plinks que pertenecen a una y s´olo una de las Aplaquetas activas. Para que un link que una los puntos iyjs´olo pertenezca a una plaqueta, debe cumplirse que sisj=−1, por tanto Ppuede escribirse P=1 2X <ij> (1 −sisj),(14) donde < ij > denota a las parejas de puntos de la red vecinos. Si ahora introducimos ambas expresiones en la funci´on de partici´on (11) llegamos a Z= 22NxN 2yNX {si} x−1 2Pisiy−1 2P<ij> sisj.(15) Considerando x, y ≥0 (es decir, βG, βI≥0), tenemos que la funci´on de partici´on viene dada por Z= 22NxN 2yNX {si} exp  −1 2log xX i si−1 2log yX <ij> sisj .(16) Salvo un factor constante, est´a es la funci´on de partici´on del modelo de Ising ferromagn´etico en una red cuadrada, ZIsing =X {si} exp FX <ij> sisj+hX i si .(17) El acoplo Fy campo magn´etico hse relacionan con nuestras variables seg´un las siguientes expresiones F=−1 2log y, h =−1 2log x. (18) De cara a estudiar el diagrama de fases de este modelo, dado que ZIsing es invariante bajo el cambio h↔ −h, la funci´on de partici´on ser´a sim´etrica bajo x↔1/x, y por tanto bastar´a estudiar la regi´on 0 < x ≤1. Adem´as para βI∈[0,∞) tenemos y∈[0,1) →F > 0, es decir: El modelo Z2gauge-Higgs es dual al modelo de Ising ferromagn´etico para d= 2. No obstante, una vez desarrollados los algoritmos para estudiar el primero, nada nos impide emplear y > 1 (F < 0) para analizar el modelo antiferromagn´etico, que tiene su propio inter´es como ya se ha comentado anteriormente. 7 Activar una plaqueta # vecinos 0 1 2 3 4 ∆A1 ∆P4 2 0 −2−4 α xy4xy2x x/y2x/y4 Tabla 1: Probabilidades del paso de Metropolis elemental al activar una plaqueta con un n´umero de vecinos activos dado. Desactivar una plaqueta # vecinos 0 1 2 3 4 ∆A-1 ∆P−4−2 0 2 4 α1/xy41/xy21/x y2/x y4/x Tabla 2: Probabilidades del paso de Metropolis elemental al desactivar una plaqueta con un n´umero de vecinos activos dado. por α=p(gnew) p(g)=x∆Ay∆P,(35) definiendo ∆A≡A(gnew)−A(g), ∆P≡P(gnew)−P(g). Sabiendo esto, introducimos el paso elemental de Metropolis como el intento de activar/desactivar una plaqueta (o en t´erminos de la dualidad, flipar un esp´ın). Al hacerlo, siempre incrementaremos o disminuiremos el ´area en una unidad. El incremento en el per´ımetro depende del n´umero de plaquetas vecinas activas en el momento de realizar el cambio. Para el caso d= 2, en el que el n´umero de vecinos activos est´a entre 0 y 4, escribimos explicitamente en las tablas 1 y 2 las probabilidades αcorrespondientes a activar o desactivar una plaqueta. El paso elemental que hemos definido intenta cambiar el estado de una ´unica plaqueta. Para conseguir un paso de Metropolis erg´odico4, es decir, para poder alcanzar cualquier configuraci´on partiendo de una dada, podemos 4En el caso x=y= 1 el sistema hace un flip colectivo. Para subsanar este caso, por otro lado trivial, puede cambiarse el paso secuencial por una elecci´on de las plaquetas aleatoria. No obstante en el caso general es preferible la elecci´on secuencial por ser m´as eficiente. 14 realizar sucesivamente Npasos elementales, secuencialmente. Esto es, en cada pasada o sweep intentamos activar/desactivar cada una de las plaquetas del sistema, una a una y por orden, pudiendo llegar as´ı a cualquier punto del espacio de configuraciones con probabilidad no nula. 3.2.2. Algoritmo PS El siguiente algoritmo que describiremos es una implementaci´on propia del algoritmo original de Prokof’ev y Svistunov [2], en adelante algoritmo PS. El objetivo es tener una referencia con la que poder comparar y en la que de antemano sabemos que no se da critical slowing down [7]. Adem´as, entender qu´e posee de diferente frente a los algoritmos locales habituales, puede ayudarnos a la hora de desarrollar generalizaciones que conserven algunos de sus beneficios. La idea fundamental de este algoritmo es cambiar el espacio de configuraciones original, en nuestro caso el conjunto de grafos G, por un espacio extendido Gext, de forma que G ⊂ Gext. En la pr´actica, esto permite al algoritmo pasar de unas configuraciones a otras dentro de G(tambi´en llamadas configuraciones de vac´ıo en la literatura) a trav´es de nuevos caminos que se abren gracias a las configuraciones a˜nadidas en Gext. Estas ´ultimas configuraciones se construyen mediante la inclusi´on de defectos en la red, entendidos como violaciones de las restricciones que se aplican para determinar qu´e grafos son v´alidos y cu´ales no. Por ejemplo, en el modelo Z2gauge-Higgs, las configuraciones de vac´ıo, que hemos llamado anteriormente grafos v´alidos, son superficies formadas por variables plaqueta (activas) y variables link en su frontera, definiendo frontera como los links que pertenecen a un n´umero impar de plaquetas (con d= 2 esto se reduce a los links que pertenecen a una ´unica plaqueta). Introducir defectos significar´ıa por ejemplo definir Gext como todas las configuraciones que cumplen la restricci´on anterior, salvo en dos links l1yl2. Si l1=l2, el defecto se corrige, y recuperamos las configuraciones de vac´ıo. En el modelo de Ising bidimensional sin campo magn´etico, el caso m´as sencillo expuesto en [2], el espacio de configuraciones base es en principio {si}, siendo la funci´on de partici´on ZIsing =X {si} exp  FX <ij> sisj .(36) Repitiendo el razonamiento de la dualidad expuesta en 2.2 en sentido opuesto, se encuentra que el modelo puede ser expresado en t´erminos de variables link σnµ. Para ello se emplea la identidad (5) de forma que ZIsing = 2N(cosh F)NlX {si}Y <ij> (1 + sisjtanh F),(37) 15 con NyNlel n´umero de sitios y links de la red, respectivamente. Al realizar la suma sobre todas las configuraciones, ´unicamente sobreviven aquellas en las que los espines de cada sitio de la red aparecen un n´umero par de veces. Haciendo σnµ =−snsn+µ, la condici´on se traduce en que las configuraciones v´alidas son caminos de links cerrados, llegando finalmente a ZIsing ∝(cosh F)NlX g∈G w(g) con w(g) = (tanh F)P,(38) con P=Pnµ(1 + σnµ)/2 el per´ımetro, o n´umero de links activos (iguales a 1, tal y c´omo los hemos definido). Ahora podr´ıamos realizar una rutina de Monte Carlo para muestrear el sistema, o tambi´en hacer equivaler cada configuraci´on de caminos cerrados con una configuraci´on de plaquetas activas, pero nos quedamos con la formulaci´on en t´erminos de links por resultar m´as adecuada para extender el espacio de configuraciones. Para ello, consideramos la funci´on de correlaci´on de dos espines en el modelo de Ising, G(i1−i2) = hsi1si2i, por definici´on igual a g(i1−i2)/Z con g(i1−i2) = X {si} si1si2exp  FX <ij> sisj .(39) Podemos proceder igual que con (36), llegando a g(i1−i2) = X {si} si1si2Y <ij> (1 + sisjtanh F).(40) Si nos preguntamos ahora qu´e tipos de configuraciones g∈ Gext contribuyen a la funci´on g(i1−i2) tras la suma sobre {si}, la respuesta es casi id´entica a la del caso anterior: todas en las que en cada sitio de la red confluyan un n´umero par de links (G), m´as aquellas en las que incidan un n´umero par de links sobre todos los sitios de la red salvo i1ei2. Por tanto en Gext, adem´as de las configuraciones de caminos de links cerrados, se permite un ´unico camino adicional, cuyos extremos sean i1ei2. Cabe destacar que g(0) y Zposeen id´enticas configuraciones. Esto hace posible que llevando a cabo un proceso de Monte Carlo para g(i1−i2), se acumule la estad´ıstica necesaria para calcular G(i1−i2) = g/Z. Los estimadores empleados en el conjunto Gext ser´an g(i) = δi,i1−i2, Z =δi1,i2,(41) siendo Gel ratio entre ambos. Una vez estimada la funci´on de correlaci´on a dos puntos, podemos obtener dos observables de los que ya se ha hablado 16 con anterioridad, energ´ıa y susceptibilidad magn´etica, esta ´ultima definida como en (25): e=−1 4X |i|=1 g(i) Z, χ0=X i g(i) Z.(42) La energ´ıa tambi´en puede obtenerse a trav´es del per´ımetro en las configuraciones de vac´ıo. Derivando la funci´on de partici´on en (38), obtenemos un segundo estimador para la energ´ıa −e≡1 2N dlog Z dF = tanh F+1−(tanh F)2 tanh F P 2N,(43) ´unicamente v´alido para i1=i2. Tambi´en es posible obtener el calor espec´ıfico en las configuraciones de vac´ıo a partir del valor esperado de PyP2. Derivando la expresi´on de la energ´ıa anterior respecto a F, obtenemos, tras un poco de ´algebra: CV= (1 −a2) + 1 2Na2−1 a2hPi+1 2N1−a2 a2hP2i−hPi2, (44) donde hemos introducido a= tanh Fpara simplificar la notaci´on. El algoritmo PS o algoritmo del gusano, es una forma de avanzar por el espacio de configuraciones Gext. Dicho de otro modo, define el paso de Metropolis a realizar para muestrear Gext. El paso m´as elemental del algoritmo se define de forma que si i16=i2, intenta desplazar el defecto i1a un sitio de la red adyacente, elegido al azar, cambiando el estado del link que une el sitio original con el nuevo. De (38) se deduce que el cambio se aceptar´a con probabilidad m´ın {1,tanh F},(45) si se intenta activar un link y m´ın {1,(tanh F)−1},(46) en el caso de que se intente desactivar un link. De aqu´ı viene la denominaci´on de algoritmo del gusano, ya que la evoluci´on de Monte Carlo consiste en caminos abiertos que van abri´endose paso por el ret´ıculo. Adem´as, si i1=i2, se intenta cambiar la posici´on de ambos defectos con cierta probabilidad p(arbitraria a la hora de dar el resultado correcto) a otro sitio de la red elegido al azar. Dependiendo de c´omo combinemos estos pasos elementales obtendremos diferentes versiones del algoritmo con diferente coste computacional. Hasta tres variantes diferentes han sido implementadas: 17 PS original. El algoritmo tal cual es presentado en [2]: si i1=i2intenta cambiar la posici´on de los defectos con probabilidad p(p= 0,5, aunque tras varias pruebas se ha comprobado que esto tiene poca influencia en el algoritmo), y con probabilidad 1 −phacer un desplazamiento de i1. Si i16=i2, siempre se intenta desplazar i1. Wolff. El algoritmo descrito por Wolff en [7], que consta de dos pasos. En el primero, siempre intenta desplazar i1. El segundo ´unicamente se lleva a cabo si i1=i2, intentando cambiar la posici´on de los defectos con probabilidad p(con probabilidad 1 −pno hace nada, a diferencia del PS original). Bic´efalo. Esta es una propuesta propia, realizada con el objetivo de lograr un mejor rendimiento del algoritmo (en t´erminos de su coste computacional). A diferencia de los anteriores, ´unicamente se emplea el paso elemental del desplazamiento de un defecto. En cada paso de Metropolis, primero se intenta desplazar i1, y acto seguido i2. De esta forma, el camino abierto de la configuraci´on o gusano, avanza por ambos extremos (de ah´ı la denominaci´on de bic´efalo). En las tres variantes el modo de tomar medidas es id´entico: se empieza la simulaci´on con una configuraci´on de vac´ıo, y se deja avanzar hasta que se vuelve a encontrar otra, es decir, hasta que i1vuelve a ser igual a i2. En un sweep como este, ´unicamente la ´ultima configuraci´on contribuye a la funci´on de partici´on, y la estimaci´on de los observables (42) se vuelve inmediata: χ es el n´umero de pasos de Metropolis dados en el sweep y eel n´umero de veces que la distancia entre i1ei2ha sido la unidad (|i|= 1) en el mismo. Dado que de las tres variantes implementadas, la que menos coste computacional tiene es la que hemos denominado bic´efala (ver figura 2), es esta la que usaremos para llevar a cabo los c´alculos de la secci´on 4, tanto para el c´alculo de los observables con h= 0, como para su comparaci´on con los otros algoritmos. 3.2.3. Algoritmo h´ıbrido Conociendo el funcionamiento de los dos algoritmos anteriores, cabe preguntarse si una mezcla de ambos podr´ıa abordar modelos m´as generales que el algoritmo PS, como nuestro Z2gauge-Higgs, y a la vez conservar las buenas propiedades de este. Un ejemplo de tal intento es el de Mercado y colaboradores en [5], donde abordan los modelos Z3yU(1) gauge-Higgs con d= 4. Siguiendo su filosof´ıa, vamos a implementar una versi´on del bautizado como surface worm algorithm (SWA), adaptado al modelo Z2 gauge-Higgs en dos dimensiones. Nuestro objetivo ser´a estudiar la presencia o no de critical slowing down, algo que queda sin determinar en la referencia mencionada. 18 10 100 16 32 64 128 Neff /tCPU (s−1) L - Tama˜no del ret´ıculo Figura de m´erito para las diferentes variantes Wolff PS Bic´efalo Figura 2: Comparaci´on del coste computacional de las diferentes variantes de algoritmos PS en el punto cr´ıtico del modelo de Ising. La figura de m´erito presentada es el n´umero de medidas efectivas estad´ısticamente independientes obtenidas en cada segundo de c´omputo, es decir Neff /tCPU =N/(2∗τint ∗tCP U ). Los tiempos de autocorrelaci´on son los correspondientes al estimador (43). Si bien la dependencia con Les similiar para los 3 casos, la variante bic´efala es m´as eficaz en todos los tama˜nos ensayados. 19 Activar una plaqueta Desactivar una plaqueta # links activados 2 0 -2 2 0 -2 ∆A1 -1 ∆P2 0 −2 2 0 −2 α xy2x x/y2y2/x 1/x 1/xy2 Tabla 3: Probabilidades del paso de Metropolis elemental al activar una plaqueta. El espacio de configuraciones Ges en este caso el mismo que en el algoritmo geom´etrico (3.2.1). Adem´as se a˜naden las configuraciones que cumplen las restricciones propias de Gen todos los links salvo en uno, LV, para formar Gext. Los observables que se miden en este caso son tambi´en los mismos que para el algoritmo geom´etrico. Dicho de otro modo, s´olo se mide en las configuraciones de vac´ıo, empleando las configuraciones de Gext ´unicamente para conseguir una evoluci´on de Monte Carlo diferente. El paso de Metropolis de este algoritmo se estructura como sigue: 1. Partimos de una configuraci´on v´alida (de vac´ıo). Se elige un link al azar, LV, y se intenta insertar un defecto (cambiar su estado de activo a inactivo, o viceversa), con probabilidad m´ın {1, α}con α=(y, para activar un link , 1/y, para desactivar un link . Si se inserta el defecto se sigue con el siguiente punto, en caso contrario el paso de Metropolis acaba aqu´ı. 2. En esta situaci´on se tiene una configuraci´on con un defecto en LV. Se escoge una direcci´on y sentido al azar. Si es paralela al defecto LV, se intenta subsanar con las mismas probabilidades que en el punto uno. Si se tiene ´exito, hemos llegado a una configuraci´on de vac´ıo y el gusano acaba, en caso contrario se repite el punto 2. Si la direcci´on no es paralela a LV, se intenta propagar el segmento mediante la adici´on de un segmento. Este segmento esta formado por una plaqueta, adyacente a LV, y dos links escogidos al azar entre los tres restantes. El link que queda sin escoger ser´ıa (si se acepta el cambio) la nueva ubicaci´on del defecto. Las probabilidades dependen del n´umero de links y plaquetas a activar/desactivar, y vienen detalladas en la tabla 3. Tanto si se acepta el segmento como si no, se vuelve a repetir el punto 2. Entre los cambios que puede llevar a cabo este paso, est´a el de activar / desactivar una plaqueta (introducir un defecto, a˜nadir un segmento y sanar 20 el defecto), por lo que la ergodicidad del algoritmo est´a asegurada. Adem´as, en cada sweep (despu´es del cual se toman las medidas de los observables, definidos en 2.3) llevamos a cabo tantos de los pasos anteriores como sites haya en la red. 4. Resultados num´ericos En esta secci´on se muestran los resultados obtenidos por los algoritmos presentados con anterioridad. En primer lugar veremos los correspondientes al caso particular x= 1, es decir, al caso sin campo magn´etico. Comparando los valores de los estimadores de cada algoritmo, podremos verificar que los hemos implementado correctamente. Hecho esto compararemos el coste computacional de cada algoritmo en el punto cr´ıtico. Para acabar, mostraremos los resultados obtenidos para el espacio de par´ametros completo, esto es: para h6= 0, incluyendo el caso antiferromagn´etico F < 0. 4.1. Caso h= 0 En este caso podemos comparar los datos de los tres algoritmos entre s´ı, analizando la dependencia de diferentes observables con el acoplo Fen un entorno del punto cr´ıtico ferromagn´etico Fc=1 2log √2−1≈0,44068679350977 . . . (47) Los casos de la energ´ıa, el calor espec´ıfico CVy la susceptibilidad magn´etica χ0, se presentan en las figuras 3, 4 y 5, respectivamente. El valor absoluto de la magnetizaci´on se obtiene ´unicamente del algoritmo geom´etrico y del h´ıbrido, y se presenta en la figura 6. Vemos en todas ellas c´omo existe un buen acuerdo entre los diferentes algoritmos, dados los errores estad´ısticos asociados (las barras de error corresponden en todos los casos a 1σ, y se han calculado siguiendo los procedimientos indicados en 3.1). Adem´as tambi´en se han comparado los valores de la energ´ıa y la susceptibilidad en el punto cr´ıtico con los obtenidos por Wolff en [7], obteniendo un acuerdo igual de satisfactorio. Considerando estas pruebas como un buen test de que los algoritmos est´an correctamente implementados, podemos comparar la eficiencia de estos en el punto cr´ıtico. Para ello introducimos la siguiente figura de m´erito: Err2 O×tCPU .(48) Dado que los errores de cada observable decrecen con 1/√N, la cantidad definida tiende a una constante que nos dice, para un tiempo de c´omputo 21 -0.95 -0.9 -0.85 -0.8 -0.75 -0.7 -0.65 -0.6 -0.55 -0.5 -0.45 0.34 0.36 0.38 0.4 0.42 0.44 0.46 0.48 0.5 0.52 0.54 hei(F) Acoplo F Curvas hei(F) para h= 0 H´ıbrido Geom´etrico Ps-Bic´efalo (42) Ps-Bic´efalo (43) Figura 3: Energ´ıa een funci´on del acoplo F, ret´ıculo 32 ×32. 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 0.34 0.36 0.38 0.4 0.42 0.44 0.46 0.48 0.5 0.52 0.54 CV(F) Acoplo F Curvas CV(F) para h= 0 H´ıbrido Geom´etrico Ps-Bic´efalo Figura 4: Capacidad calor´ıfica CVen funci´on del acoplo F, ret´ıculo 32 ×32. 22 0 100 200 300 400 500 600 700 800 900 1000 0.34 0.36 0.38 0.4 0.42 0.44 0.46 0.48 0.5 0.52 0.54 χ0(F) Acoplo F Curvas χ0(F) para h= 0 H´ıbrido Geom´etrico Ps-Bic´efalo Figura 5: Susceptibilidad χ0en funci´on del acoplo F, ret´ıculo 32 ×32. 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.34 0.36 0.38 0.4 0.42 0.44 0.46 0.48 0.5 0.52 0.54 h|m|i(F) Acoplo F Curvas h|m|i(F) para h= 0 H´ıbrido Geom´etrico Figura 6: Magnetizaci´on h|m|i en funci´on del acoplo F, ret´ıculo 32 ×32. 23 Referencias [1] N. Prokof’ev, B. Svistunov and S. Tupitsyn, Worm Algorithm in Quantum Monte Carlo Simulations,Phys. Lett. A 238 (256) (1998). [2] N. Prokof’ev and B. Svistunov, Worm Algorithms for Classical Statistical Models,Phys. Rev. Lett. 87 (16) 160601 (2001). [3] V. Azcoiti, E. Follana, A. Vaquero, and G. Di Carlo, Geometric algorithm for abelian-gauge models,JHEP 0908 008 (2009). [4] C. Gattringer and A. Schmidt, Gauge and matter fields as surfaces and loops: An exploratory lattice study of the Z(3) gauge-Higgs model,Phys. Rev. D 86 094506 (2012). [5] Y. D. Mercado, C. Gattringer, and A. Schmidt, Surface worm algorithm for abelian Gauge–Higgs systems on the lattice Comput. Phys. Commun. 184 1535 (2013). [6] Y. Deng, T. M. Garoni and A. D. Sokal, Dynamic critical behavior of the worm algorithm for the Ising model,Phys. Rev. Lett. 99 110601 (2007). [7] U. Wolff, Simulating the all-order strong coupling expansion I: Ising model demo,Nucl. Phys. B 810 491 (2009). [8] U. Wolff, Monte Carlo errors with less errors,Comput. Phys. Commun. 156 143 (2004). [9] Matthias Troyer and U. J. Wiese, Computational Complexity and Fundamental Limitations to Fermionic Quantum Monte Carlo Simulations, Phys. Rev. Lett. 94 170201 (2005). [10] Y. D. Mercado, C. Gattringer, and A. Schmidt, Dual Lattice Simulation of the Abelian Gauge-Higgs Model at Finite Density: An Exploratory Proof of Concept Study Phys. Rev. Lett. 111 141601 (2013). [11] V. Azcoiti, G. Cortese, E. Follana, M. Giordano, A geometric Monte Carlo algorithm for the antiferromagnetic Ising model with “topological” term at θ=π,Nuclear Physics B 883 656 (2014) [12] D. C. Rapaport and C. Domb, The smoothness postulate and the Ising antiferromagnet,J. Phys. C 4 2684 (1971). [13] E. Muller-Hartmann and J. Zittartz, Interface Free-Energy and Transition-Temperature of Square-Lattice Ising Antiferromagnet at Finite Magnetic-Field,Z. Phys. B 27 261 (1977). [14] X.-Z. Wang and J. S. Kim, The Critical Line of an Ising Antiferromagnet on Square and Honeycomb Lattices,Phys. Rev. Lett. 78 413 (1997). 30 [15] D. C. Rapaport, Monte Carlo Study of the Phase Boundary of the Ising Antiferromagnet,Phys. Lett 65 A 147 (1978). [16] K. Binder and D. P. Landau, Phase diagrams and critical behavior in Ising square lattices with nearestand next-nearest-neighbor interactions,Phys. Rev. B21 1941 (1980). [17] F. J. Wegner, Duality in Generalized Ising Models and Phase Transitions without Local Order Parameters, J. Math. Phys. 12 2259 (1971). [18] R. Balian, J. R. Drouffe and C. Itzykson, Gauge Fields on a lattice. II. Gauge-invariant Ising model,Phys. Rev. D 11 2098 (1975). [19] N. Metropolis et al, Equation of State by Fast Computing Machines,J Chem. Phys. 21 6 1087 (1953). [20] R. H. Swendsen and J.-S. Wang, Nonuniversal critical dynamics in Monte Carlo simulations,Phys. Rev. Lett. 58 86 (1987). [21] N. Madras and A. D. Sokal, The pivot algorithm: A highly efficient Monte Carlo method for the self-avoiding walk,J. Stat. Phys. 50 109 (1988). [22] B. Efron, The Jackknife, the Bootstrap and Other Resampling Plans (Society for Industrial and Applied Mathematics [SIAM], Philadelphia, 1982). [23] R.G. Miller, The Jackknife–A review,Biometrika 61 1 (1974). [24] V. Azcoiti, G. Di Carlo, and A. F. Grillo, New proposal for including dynamical fermions in lattice gauge theories: The compact-QED case, Phys. Rev. Lett 65 2239 (1990) [25] V. Azcoiti, V. Laliena, X. Q. Luo et al, Microcanonical fermionic average method for Monte Carlo simulations of lattice gauge theories with dynamical fermions,Phys. Rev. D 48 402 (1993) 31