Repositorio Institucional de Documentos

2021 116 Eduardo Royo Amondarain Non-perturbative physics in lattice gauge theories Director/es Azcoiti Pérez, Vicente Follana Adín, Eduardo © Universidad de Zaragoza Servicio de Publicaciones ISSN 2254-7606 Eduardo Royo Amondarain NON-PERTURBATIVE PHYSICS IN LATTICE GAUGE THEORIES Director/es Azcoiti Pérez, Vicente Follana Adín, Eduardo Tesis Doctoral Autor 2021 UNIVERSIDAD DE ZARAGOZA Escuela de Doctorado Programa de Doctorado en Física Repositorio de la Universidad de Zaragoza – Zaguan http://zaguan.unizar.es Non-perturbative physics in lattice gauge theories Doctoral dissertation of Eduardo Royo Amondarain Supervised by Vicente Azcoiti P´ erez Eduardo Follana Ad´ın UNIVERSIDAD DE ZARAGOZA Departamento de F´ ısica Te´ orica & Centro de Astropart´ ıculas y F´ ısica de Altas Energ´ ıas 2020 Repositorio de la Universidad de Zaragoza – Zaguan http://zaguan.unizar.es Aknowledgements Agradecimientos Llevar a buen puerto el trabajo que aqu´ı se presenta, y que ha sido origen de no pocos desvelos, habr´ıa sido imposible de no contar con el apoyo de muchas personas. Dada la dificultad de nombrar a todas ellas en tan pocas l´ıneas, querr´ıa empezar por dar las gracias de forma general a quien, desde su propio ´ambito y en mayor o menor medida, haya contribuido en esta empresa. Sin perjuicio de lo anterior, y rogando se disculpe cualquier omisi´on que pueda producirse a continuaci´on, en lo que sigue quisiera destacar algunos nombres que considero especialmente relevantes. En lo acad´emico, resulta ineludible reconocer la labor del Profesorado con el que he tenido la oportunidad de coincidir en las diferentes etapas del Sistema Educativo. Buenos ejemplos ser´ıan mi profesora de Lengua en 1ode la E.S.O. o mi tutor durante ese mismo curso, de Dibujo, sin olvidar la profesora que, en 2ode la E.S.O., me involucr´o en la Olimpiada Matem´atica. Del mismo modo, ya durante la etapa universitaria, pude disfrutar de un gran elenco de docentes que supieron transmitir su pasi´on por esta disciplina, con menci´on especial para el Departamento de F´ısica Te´orica, cuyas asignaturas siempre dejaban con ganas de seguir indagando en la materia. El agradecimiento al Departamento es doble, de hecho, pues su acogida y acompa˜namiento a lo largo de estos a˜nos han sido dignas de elogio. Ha sido para m´ı un aut´entico honor habitar, siquiera temporalmente, los pasillos de Te´orica, en los que siempre me han hecho sentir como en casa. A caballo entre lo acad´emico y lo personal, tengo que dar las gracias a mis directores, Vicente y Eduardo, cuya gu´ıa por los no-siempre-claros caminos del doctorado ha sido encomiable. De ellos he podido aprender mucho, pero por si esto fuera poco, su contribuci´on trasciende ampliamente lo profesional. Adem´as de ser buenos colegas, en el sentido m´as formal del t´ermino, tambi´en han demostrado ser buenos amigos. En este mismo sentido, no podr´e agradecer lo suficiente la excelente acogida de Giuseppe durante mis estancias en L’Aquila. Muy dif´ıcilmente se puede dar en la vida con una persona de mayor calidad humana. No quiero olvidarme tampoco de otros compa˜neros ocasionales de viaje, como Alex o Matteo, ni por i ii supuesto de mi compa˜nero de fatigas y despacho, Filiberto. Ya en el terreno personal, gracias a todas las amistades que han estado ah´ı en tantos momentos. A mis compa˜neros del Stadium Casablanca, con los que durante muchos a˜nos he podido disfrutar del balonmano. A los del Instituto y a los de la carrera, con menci´on especial al eje reformista-populista. A las Colonas, reino animal incluido, que sol´ıamos ser de Cat´an—y que cada vez somos m´as. Y, por supuesto, gracias a mi familia. Gracias a Menchu, Musta, Nora y Avi, que siempre me han hecho sentir como uno m´as. Gracias a Arrate, Jes´us, Alejandro y Gato, porque decir que siempre me hab´eis apoyado incondicionalmente ser´ıa quedarse corto. Hab´eis sido el mejor apoyo que un hijo o un hermano puede tener. Gracias tambi´en, Robin, por sacarme a pasear todos los d´ıas, pues no hay mejor manera de darle una vuelta a un problema rebelde. Finalmente y por encima de todo, gracias a Cecilia y a Yasmina, sobre cuyo trabajo de cuidados descansa el grueso de esta tesis. Ambas hab´eis sido condici´on sine qua non para que este trabajo haya podido salir adelante, tanto como lo sois para m´ı. Va por vosotras. iii Este trabajo ha sido financiado por el Ministerio de Econom´ıa y Competitividad a trav´es de la ayuda predoctoral para la formaci´on de doctores BES-2013-063567. De igual modo, cabe mencionar los proyectos del Ministerio de Econom´ıa y Competitividad/Fondo Europeo de Desarrollo Regional FPA2012-35453 y FPA201565745-P, y de la Diputaci´on General de Arag´on/Fondo Social Europeo 2015-E24/2. 4CONTENTS underlying interactions of quarks and gluons in the QCD framework, constitutes a very active field of research, that includes a large variety of approaches [24]. The interest that this matter generates is divided between the behavior of the coupling in the infrared region—where nowadays there is no consensus about αS(q2→0) tending to zero, freezing or even diverging—and at large momenta, where perturbative QCD can be applied and both experimental and theoretical methods try to provide the most accurate approximation. In this context, lattice-based strategies have been capable of delivering results both in the infrared region and in the high energy regime, where in fact they provide the most precise determination of αS(MZ) [2]. Our work can be framed precisely into these approaches that come from lattice QCD, and it relies upon a ghost-gluon vertex computation, as in [25–27]. This thesis is orgnized as follows. In Chapter 1 we review the fundamentals of some of the topics presented in this introduction, including a brief historical review of the lattice approach. A motivation for the inclusion of a θterm in the QCD action is developed in Chapter 2, together with a brief review of the existing approaches to this problem. Chapter 3, based on the work presented in [28], is dedicated to the Ising model within an imaginary magnetic field. Hereafter we address the study of the massive Schwinger model with a θterm. The 1-flavor results, relying on the work presented in [29], are reported in Chapter 4, whereas the 2-flavor case is treated with a pseudofermions approach in Chapter 5. Finally, Chapter 6 covers the computation of αs(q2) via the ghost-gluon vertex. Lastly, our conclusions are summarized in the homonymous chapter and technical details concerning the computation of the cumulant expansion performed in Chapter 3 are given in Appendix A. Introducci´on Varias d´ecadas han pasado desde que la cromodin´amica cu´antica (QCD, por sus siglas en ingl´es) se estableci´o como la teor´ıa que describe las interacciones fuertes. Es ampliamente aceptada como una de las teor´ıas m´as exitosas de la f´ısica moderna y ha sido puesta a prueba de forma exhaustiva, tanto desde el punto de vista te´orico como del experimental. A altas energ´ıas, QCD es asint´oticamente libre, lo que significa que sus constituyentes fundamentales, quarks ygluones, interaccionan con una intensidad que decrece conforme la energ´ıa alcanza escalas m´as altas. En esta situaci´on, resulta factible aplicar la teor´ıa de perturbaciones para resolver las interacciones a corta distancia. Al estudiar secciones eficaces en procesos de altas energ´ıas, es necesario tener en cuenta las interacciones producidas tanto a corta como a larga distancia. Sin embargo, con la ayuda de los teoremas de factorizaci´on [1], es posible combinar QCD perturbativa con cierto input no perturbativo, proveniente de fuentes emp´ıricas o te´oricas, y obtener as´ı predicciones que pueden ser confrontadas experimentalmente. De hecho, este programa ha sido llevado a cabo por una amplia comunidad de cient´ıficos y cient´ıficas, trabajando en cientos de Universidades, laboratorios y habitualmente organizados en grandes colaboraciones internacionales, incluyendo varios aceleradores de part´ıculas. La importante cantidad de evidencias recogidas en este proceso ha permitido a QCD convertirse en un componente muy fiable del actual Modelo Est´andar de f´ısica de part´ıculas. Adem´as, QCD perturbativa es hoy en d´ıa un campo de investigaci´on muy activo, toda vez que centros como el Gran Colisionador de Hadrones proveen de nuevos datos experimentales cada a˜no [2], y debe invertirse un gran esfuerzo te´orico en el c´alculo de los primeros t´erminos de las expansiones correspondientes. Por otro lado, para escalas de energ´ıa no tan elevadas, la interacci´on fuerte no puede ser reducida a una serie convergente de diagramas de Feynman. De hecho, una de sus propiedades caracter´ısticas es el llamado confinamiento de color. Esto significa que los quarks se encuentran siempre (excluyendo r´egimenes de alta densidad) en estados ligados, llamados hadrones, que son neutrales respecto del color. En esta situaci´on puramente no perturbativa, hay pocas t´ecnicas que puedan analizar la teor´ıa con ´exito. Probablemente la que mejor establecida est´a es QCD 5 6CONTENTS en el ret´ıculo, denominada comunmente lattice QCD. Desde el trabajo fundacional de Wilson en 1974 [3], el ´exito del m´etodo ha ido creciendo con el tiempo. Si bien durante los primeros a˜nos realizar los c´alculos necesarios para extraer resultados significativos de QCD parec´ıa muy lejano, el progresivo refinamiento de los algoritmos junto con el crecimiento exponencial de la capacidad computacional mundial dio la vuelta a la situaci´on. Muchos hitos han sido ya alcanzados: simulaciones precisas incluyendo los efectos de loops de quarks virtuales [4], la determinaci´on del espectro de hadrones ligeros con errores sistem´aticos totalmente controlados [5] o, m´as recientemente, la computaci´on de los splittings de isosp´ın (esto es, de las diferencias de masa entre neutr´on y prot´on, u otros canales hadr´onicos) con gran acuerdo con los datos experimentales, incluso excediendo su precisi´on en algunos casos [6]. Todos ellos son buenos ejemplos del ´exito de este enfoque. Por los motivos citados, QCD es considerada como la teor´ıa que describe correctamente la interacci´on fuerte, tanto para altas como para bajas energ´ıas, y lattice QCD es reconocido por la comunidad como un m´etodo ab initio fiable que tiene una interacci´on ´util con lo experimental, parafraseando a Wilson [7]. Es tentador pensar que, con la evoluci´on actual de la potencia de c´alculo, ser´ıa simplemente cuesti´on de tiempo que el enfoque de la lattice enfrentase y resolviese cada uno de los problemas no perturbativos que todav´ıa esperan una soluci´on. Por supuesto, las cosas no son tan sencillas, e incluso cuando para un subconjunto de problemas bastar´ıa con contar con equipos m´as potentes, existen temas fundamentales que hoy en d´ıa constituyen preguntas abiertas. Al menos dos problemas comparten este estatus: el comportamiento de la materia a densidad bari´onica finita—incluyendo su diagrama de fases de temperatura-densidad—y los estudios que involucran efectos topol´ogicos en QCD. Aunque, como veremos, los intentos han sido numerosos, los avances en ambos campos han sido escasos. La principal dificultad detr´as de este modesto progreso en ambas ´areas es la misma: la acci´on de la teor´ıa es compleja, y no existe reformulaci´on conocida que pueda evitar la aparici´on de un problema de signo severo (SSP, por sus siglas en ingl´es). El SSP puede ser definido como el problema de evaluar num´ericamente la integral de una funci´on muy oscilatoria, que adem´as depende de un gran n´umero de variables. Pertenece a la clase de problemas NP-completos [8], lo que significa que encontrar un algoritmo que superase el SSP en tiempo polin´omico equivaldr´ıa a probar que P=NP, uno de los siete Problemas del milenio propuestos por el Clay Mathematics Institute. En otras palabras, si un algoritmo as´ı existiese, resolver´ıa todos los problemas de clase NP en tiempo polin´omico. Pero, hasta el d´ıa de hoy, el dilema P vs NP contin´ua sin respuesta, y consecuentemente no existe soluci´on general para el SSP que permita la aplicaci´on de las t´ecnicas de Montecarlo usuales. As´ı pues, los esfuerzos dedicados a intentar superar esta dificultad se concentran en la elaboraci´on de diferentes m´etodos, adaptados a las CONTENTS 7 particularidades del problema en cuesti´on. Este es el caso de QCD con componente imaginaria en la acci´on, que aparece cuando est´an presentes los t´erminos correspondientes al potencial qu´ımico (µ > 0) o a efectos topol´ogicos (θ6= 0). Con el objetivo de avanzar en la comprensi´on de QCD con acci´on compleja, varias propuestas se han desarrollado a lo largo de las ´ultimas d´ecadas. La din´amica de Langevin compleja [9–11], los m´etodos desarrollados por Azcoiti et al [12,13] y, m´as recientemente, los dedales de Lefschetz [14–16] y el m´etodo de la densidad de estados [17–19], son ejemplos relevantes de algunos de estos intentos. En general, estas estrategias han permitido estudiar un n´umero considerable de toy models y, en algunos casos, han obtenido resultados interesantes para QCD con potencial qu´ımico finito. Sin embargo, en el caso de θQCD no se ha podido obtener casi ning´un progreso, debido bien a las diferentes limitaciones fundamentales de las que los m´etodos citados adolecen, bien por las dificultades pr´acticas que sus respectivas implementaciones implican. En este contexto, la parte principal de esta tesis se ha dedicado al estudio de modelos que sufren de un SSP, como el modelo de Ising bidimiensional con campo magn´etico puramente imaginario, o el modelo de Schwinger masivo con un ´unico flavor y un t´ermino θ. En el primer caso, estudiamos el conocido modelo por medio de t´ecnicas anal´ıticas, explorando una regi´on del espacio de par´ametros (temperatura real y campo magn´etico imaginario) algo desatendida en la literatura, posiblemente debido a la dificultad de aplicar t´ecnicas tanto anal´ıticas como num´ericas. De acuerdo con los pocos trabajos que exploran este sistema [20,21], se espera una estructura de fases con cierta riqueza. Nuesto trabajo pretende probar la validez de uno de los m´etodos de Azcoiti et al [13,21] en este escenario y, dado que puede arrojar algo de luz en una regi´on cr´ıtica que sufre de un SSP (que ha obstaculizado la consecuci´on de una soluci´on num´erica) esperamos que pueda servir como referencia para otros m´etodos que aspiren a superar el problema del signo en cualquier teor´ıa gauge en el ret´ıculo. Con el objetivo de enfrentar sistemas similares a QCD con un t´ermino θ, y de este modo desarrollar los m´etodos que lidian con el SSP, hemos estudiado el modelo de Schwinger masivo con un flavor y t´ermino θ, que se corresponde con QED en dimensi´on 1 + 1. Como veremos, comparte un gran n´umero de propiedades con QCD, incluyendo el confinamiento y, en cierto modo, la libertad asint´otica, por lo que de hecho es ampliamente usado como su toy model. Adem´as, definir la carga topol´ogica en este modelo es casi trivial, en contraste con cualquiera de las definiciones usuales para este observable en QCD, que resultan mucho m´as intrincadas. Este hecho nos ha permitido explorar el modelo con un coste computacional factible, obteniendo resultados compatibles con la predicci´on anal´ıtica de Coleman [22] y, lo que es m´as importante, poninedo a prueba el m´etodo desarrollado en [13] en una teor´ıa gauge con fermiones, lo que constituye un paso 8CONTENTS importante en el camino a su aplicaci´on en QCD. Como derivado de la l´ınea de trabajo anterior, y empujados por la necesidad de una mayor optimizaci´on de los algoritmos anteriores, tambi´en hemos analizado la versi´on de 2 flavors del modelo de Schwinger. En este caso, se ha evitado el c´alculo completo del determinante fermi´onico siguiendo un enfoque basado en la t´ecnica de los pseudofermiones [23]. M´as all´a del estudio de sistemas afectados por un problema de signo, otro tema, incluido dentro de lattice QCD, se ha tratado en esta tesis: el running coupling αS. La dependencia de αS(q2) con el momento transferido q, que codifica las interacciones subyacentes de quarks y gluones en el marco de QCD, constituye un campo de investigaci´on muy activo, que incluye una gran variedad de metodolog´ıas [24]. El inter´es que esta materia genera se divide entre el comportamiento de este acoplamiento en la regi´on infrarroja, donde hoy en d´ıa no existe consenso sobre si αS(q2→0) tiende a cero, se congela o incluso diverge, y a momentos altos, donde QCD perturbativa puede ser aplicada y m´etodos tanto te´oricos como experimentales intentan proveer la aproximaci´on m´as precisa. En este contexto, las estrategias basadas en el ret´ıculo han sido capaces de ofrecer resultados para ambos casos, consiguiendo de hecho la determinaci´on m´as precisa para αS(MZ) [2]. Nuestro trabajo puede ubicarse precisamente entre los enfoques que vienen de la lattice, y se apoya en un c´alculo del v´ertice ghost-gluon, como en [25–27]. Esta tesis se organiza como sigue. En el Cap´ıtulo 1 repasamos lo fundamental de algunos de los temas presentados en esta introducci´on, incluyendo un breve resumen hist´orico sobre el origen de lattice QCD. Una motivaci´on para la inclusi´on del t´ermino θen la acci´on de QCD se desarrolla en el Cap´ıtulo 2, junto con un breve repaso de de los enfoques existentes para este problema. El Cap´ıtulo 3, basado en el trabajo presentado en [28], se dedica al modelo de Ising en un campo magn´etico puramente imaginario. A continuaci´on enfrentamos el estudio del modelo de Schwinger masivo con t´ermino θ. Los resultados para el caso de un ´unico flavor, presentados en [29], se ofrecen en el Cap´ıtulo 4, mientras que el caso de 2 flavors se trata con un m´etodo basado en pseudofermiones en el Cap´ıtulo 5. Finalmente, el Cap´ıtulo 6 cubre la computaci´on de αs(q2) mediante el v´ertice ghost-gluon. Para acabar, nuestras conclusiones se resumen en el cap´ıtulo hom´onomio, y algunos detalles t´ecnicos, que conciernen al c´omputo de la expansi´on de cumulantes realizada en el Cap´ıtulo 3, se discuten en el Ap´endice A. Chapter 1 The lattice approach In this chapter we cover the essential points of the lattice approach to QCD, including a brief historical review of its birth and evolution over the past few decades. The main aspects of the formalism are explained, discussing the strenghs and limitations of Monte Carlo methods when studying lattice gauge theories. Finally, some considerations about the type of errors associated with this methodology are discussed, recalling how we can control them and, eventually, in which way we can provide a precise estimation of a given observable. 1.1 The dawn of color The appearance of quarks and gluons as the fundamental constituents of baryons and mesons is relatively recent. In order to give the necessary context, it is desirable to go back to the middle of the last century. The discovery of the pion through cosmic ray experiments in 1947 [30], which was soon followed by those of the first strange particles, the kaon and the Λ, marked the beginning of a tendency that continued during the 50s, by means of different experiments involving particle colliders. This experimental fact, i.e., the discovery of a large number of new particles and resonances that were somehow related, claimed for an explanation in terms of a reduced set of degrees of freedom. It was in 1961 when Murray Gell-Mann (who had previously introduced the strangeness as a quantum number, conserved by both electromagnetic and strong interactions) provided a successful explanation of the so-called particle zoo, in his famous The Eightfold Way [31]. By extending the SU(2) isospin symmetry, Gell-Mann proposed a SU(3) flavor symmetry, broken by mass differences, that was capable to organize all the observed hadronic states. Moreover, this symmetry predicted the existence of the Ω−baryon, a particle that was observed three years later, with a mass that matched accurately the value anticipated by the model [32]. Precisely Gell-Mann in 1964, and independently 9 10 CHAPTER 1. THE LATTICE APPROACH Zweig [33], proposed that both mesons and baryons were composed of more fundamental constituents (called quarks by the mentor, and aces by his pupil) which hold fractional electric charge. At this point, some questions still remained to complete the quark puzzle. Particularly relevant was the fact that the Ω−particle, composed of three squarks in its ground state, should have a symmetric wave function. This was in contradiction with Pauli exclusion principle, which required an antisymmetric wave function under the exchange of two quarks, which were assumed to be fermions of spin 1/2. Another concern involved the decay amplitudes predicted by the quark model, which differed from the values measured in electron-positron colliders. Both issues were addressed by Gell-Mann, Fritzsch and Bardeen, who during 1971 and 1972 introduced a new exact SU(3) symmetry for the quarks, named color [34–36], which was exactly conserved. Color was soon interpreted as a gauge symmetry and a field theory was constructed by Gell-Mann, Fritzsch and Leutwyler [37] and independently by Gross and Wilczek [38] in 1973, who emphasized (as Politzer in [39]) one of the characteristic properties of the new theoretical artifact: asymptotic freedom. Later, Fritzsch and Gell-Mann would finally give the theory its modern name: Quantum Chromodynamics. In order to complete the picture of QCD, another phenomenon had to be explained: why experiments only measure color singlets, i.e., baryons and mesons are free of any color charge, and individual quarks are not present in nature as free particles. To explain this confinement of quarks and gluons into hadrons, Wilson demonstrated in 1974 how lattice gauge theories confine charged states in the strong coupling limit [3]. This work set the basis for the qualitative understanding of color confinement—and in this way closes the early period in which QCD was constructed—but its influence over the theory would transcend by far this objective, since the regularization developed by Wilson opened a whole new field within high energy physics. With respect to confinement itself, it should be noted that, although there exist broad numerical evidence of its validity, today a rigorous mathematical proof is still lacking. 1.2 A discretized spacetime The formulation of Quantum Chromodynamics as the dynamical theory of the strong interaction, together with the understanding of asymptotic freedom and confinement, is without doubt one of the key milestones reached during the past century. Even so, and due to the strongly coupled nature of the theory, the predictive power of QCD was severely limited during its first steps. Specially at low energies, perturbative methods—developed with great success for processes involving Quantum Electrodynamics—were of no use, making observables such as the 1.2. A DISCRETIZED SPACETIME 11 hadron spectrum unreachable. In this context, Wilson introduced—in which is considered to be the foundational work of the lattice approach [3]—a novel non-perturbative regularization of QCD. The essential ideas of the procedure, applicable to any gauge theory, have endured up to the present day, and are relatively simple. Starting from Feynman’s path integral formalism, spacetime is discretized in a four-dimensional (Euclidean) hypercube. All objects are now defined in the sites of the hypercube, or the links joining them. This procedure is done in such a way that gauge symmetry is exactly preserved. Although Lorentz (or Euclidean) invariance is lost for any finite lattice spacing a, Wilson argued that this obstacle could be overcome by means of a renormalization-group approach, possible if there exists a critical point at some value of the gauge-coupling of the theory. At this point, and even if lattice QCD had elucidated the confinement phenomenon, it was unclear how the approach could be exploited in order to calculate relevant physical observables. The turning point came by the end of the decade, when a procedure widely used at the time in Statistical Mechanics was introduced into the realm of Field Theory. In 1978, Wilson [40] proposed to apply Monte Carlo methods, namely a Metropolis algorithm, to lattice gauge theories, in order to elucidate numerically the confinement phenomenon. With a detailed prescription of how observables should be measured in such a framework, he also pointed out that the first calculations were already ongoing for the SU(2) gauge theory. The first numerical results, that demonstrated the potential of the new technique, came by Creutz, Jacobs and Rebbi in 1979 [41], who performed a simulation of the gauge Z(2) model in four dimensions, observing a first-order transition that supported the confinement hypothesis by countering a previous conjecture due to Migdal [42]. In the same year, Wilson [43] presented another renormalizationgroup approach to the SU(2) gauge theory, importing block-spin techniques from Statistical Mechanics. Already in 1980, Creutz [44] provided further evidence of both confinement and asymptotic freedom in the SU(2) gauge theory, ensuring in this way the possibility of taking the continuum limit by holding constant a physical observable like the string tension. These novel ideas crystallized in several works over the years to follow, including the pioneering computations of the rho meson mass for a discrete approximation of the SU(2) theory, by Weingarten [45], or the calculation of several hadron masses due to Hamber and Parisi [46], performed already with SU(3) as the gauge group. In both cases, the limits imposed by computational resources were alleviated by the use of the so-called quenched approximation, which neglects entirely the effect of virtual quark loops—with the associated addition of uncontrolled errors. The inclusion of full dynamical fermions, i.e., taking into account both gluon and quark dynamics, was achieved first in 1983 by Azcoiti and Nakamura [47], who, leaning 12 CHAPTER 1. THE LATTICE APPROACH on the pseudofermions method proposed by Fucito et al [23], computed the mass splitting of the rho and omega mesons for the icosahedral approximation of SU(2) as the gauge group. A similar computation, also relying in the pseudofermionic approach, was made by Hamber in 1985 [48], already using SU(3) as the gauge symmetry group. These early studies, even when they were made with very modest computational power by today standards, and with a number of sources of uncontrolled systematic errors—especially significant for those made within the quenched approximation—were successful in proving the feasibility of the lattice approach. During the subsequent decades, computing resources kept growing steadily. In parallel, new algorithmic and theoretical advances were developed, which ultimately made possible breakthroughs such as the computation of the light hadron spectrum from first principles [5] or, more recently, the determination of the isospin mass splittings [6]. At present time, there exist many research groups actively working in this area; state-of-the-art results can be reviewed annually in the proceedings of the International Symposium on Lattice Field Theory. For more details about the historical development of the lattice approach, including a more technical discussion on the issue, we refer the interested reader to the extensive review of Fodor and Hoelbling [49]. 1.3 The formalism At this point, and before discussing some of the caveats of the lattice approach, it seems desirable to sketch at least the basic premises of the method. In any case, for a more exhaustive introduction to the topic, we recommend the interested reader the lecture notes of Davies [50] or the more recent book by Gattringer and Lang [51], which in fact has served as a reference for some of the topics presented in this section. 1.3.1 A finite path integral The starting point of the approach is the path integral formalism [52,53], which, for a given quantum field theory, allows to write its partition function as a functional integral Z=ZDΦe−S[Φ],(1.1) where Sis the euclidean action of the theory—even being considered at imaginary time, physical information can be recovered as long as the theory fulfills a set of axioms [54,55]—and the integration is meant to be performed over all possible configurations of the fields of the theory, denoted generically by Φ. If we particularize 1.3. THE FORMALISM 13 for QCD, the fields to be considered are quarks ψ, antiquarks ¯ ψand gluons Aµ, and (1.1) can be formulated schematically as ZQCD =ZDψD¯ ψDAµe−SQCD[ψ, ¯ ψ,Aµ],(1.2) with vacuum expectation values of operators O(ψ, ¯ ψ, Aµ) being given by hOi =1 ZQCD ZDψD¯ ψDAµO(ψ, ¯ ψ, Aµ)e−SQCD[ψ, ¯ ψ,Aµ].(1.3) The precise meaning of the integration symbol present in (1.2) and (1.3) is a subtle mathematical issue; at this point the theory would be ill-defined and a regularization is mandatory in order to extract any physical information and avoid divergences. To this end, there are several possibilities that lean on perturbative expansions in the coupling, such as Pauli-Villars (or dimensional) regularization [56], which have been widely used both in QED and QCD. The prescription of Wilson that was introduced earlier [3], is however the only fully non-perturbative regularization known, allowing to study quantum field theories from first principles, especially QCD beyond short-distance interactions, where confinement—together with the highly non-trivial structure of the QCD vacuum—poses an insurmountable obstacle to perturbative-based approaches. Thereby, going ahead with QCD, we discretize the continuum four-dimensional euclidean space-time in an hypercubic grid with constant lattice spacing a; this parameter will be the regulator of the theory, in which limit a→0 its original version is recovered. Moreover, we restrict the full spacetime to a simple 4-d box of finite extent, limiting in this way the spatial volume and the imaginary-time evolution of the system. By doing so, the number of variables to consider becomes finite and the original problem can begin to be pictured as computationally treatable. If we return to (1.2), the fields ψ(x) and ¯ ψ(x) now take values only for x=a(n1, n2, n3, n4), where niare integers satisfying 0 ≤nia < Li,Libeing the spatial (temporal) extent of the 4-d box in each dimension. With respect to the gluon field Aµ(x), it is convenient to postpone briefly its definition on the lattice. First, we should note that the action of the theory, SQCD in (1.2), is the spacetime integral of the following euclidean1lagrangian density, LQCD =¯ ψ(iγµDµ+m)ψ+1 2g2Tr (FµνFµν),(1.4) 1The ordinary real-time version of LQCD, i.e., with Minkowski spacetime metric, is just ¯ ψ(iγµDµ−m)ψ−Tr (Fµν Fµν)/2g2.Note that we are using a rather compact notation, where only spacetime indexes are explicit, while those corresponding to color, spin and flavor are kept implicit. A more detailed notation could label quark (and antiquark) field components by ψf α,c, since ψis a 3-color vector, a 4-Dirac spinor and a nf-flavor vector. In the same way, mis a diagonal matrix in flavor space, containing each one of the quark masses. 20 CHAPTER 1. THE LATTICE APPROACH prescriptions that minimize further the order of the errors, the original proposal of Wilson provides a good balance between the precision achieved and the complexity of its definition. In any case, improved gauge actions are built by adding higher order terms to (1.17), involving loops larger than the plaquette that are suppressed by powers of a. This type of gluonic actions are state-of-the-art in modern lattice computations. 1.3.3 Interlude: chiral symmetry in the continuum Just before addressing the construction of the fermionic matrix M, we found necessary to make at least a shallow review of the role of chiral symmetry in QCD, as its formulation on the lattice has deep and direct implications in how dynamical fermions can enter into the lattice action. In its continuum formulation, the action of QCD for Nfflavors is given by the spacetime integral of the lagrangian (1.4). Focusing on the fermionic part and keeping a compact notation, where ψand ¯ ψfields are understood to be Nf-vectors in flavor space, we have SQCD-F =Zd4x¯ ψ(iγµDµ+m)ψ. (1.18) For the particular case of zero fermion mass, i.e., if all the components of the diagonal matrix mvanish, there exist two sets of transformations that leave the above action invariant. The first family is composed by the so-called vector transformations, given by ψ→eiα1ψ¯ ψ→¯ ψe−iα1,(1.19) ψ→eiαTiψ¯ ψ→¯ ψe−iαTi,(1.20) where the Timatrices are the generators of the flavor group, SU(Nf), and so, the index iruns from 1 to N2 f−1. When considering altogether (1.19) and (1.20), the group extends to U(Nf). The other family of transformations can be constructed by including a γ5matrix in the previous ones, and thus they are given by ψ→eiαγ51ψ¯ ψ→¯ ψe−iαγ51,(1.21) ψ→eiαγ5Tiψ¯ ψ→¯ ψe−iαγ5Ti.(1.22) This last set receives the name of axial transformations. When considered together with vector transformations, they receive the name of chiral transformations; in the same way, the massless limit—in which they preserve SQCD—is called chiral limit. To differentiate the two sectors that compose chiral symmetry, it is customary to label them with a Vor a Asubscript. In this way, the full chiral group is given by SU(Nf)V×SU(Nf)A×U(1)V×U(1)A.(1.23) 1.3. THE FORMALISM 21 It is important to note that axial transformations are symmetries of the massless action, but this symmetry is explicitly broken for finite mass quarks. On the contrary, vector transformations (1.20) preserve the action also in the degenerate mass case—i.e., Nfspecies of equal non-zero mass—conforming the well-known isospin symmetry. Furthermore, (1.19) is always a symmetry of the action, which implies baryon number conservation. Equally relevant is the so-called axial anomaly. Although U(1)Ais a symmetry of the massless action, it is explicitly broken in the fully quantized theory. This reduces the full expression (1.23) to the following symmetry group: SU(Nf)V×SU(Nf)A×U(1)V.(1.24) This group is usually expressed in slightly different terms. To this aim, it suffices to recall that chiral symmetry splits both the quark fields and the action into two separate pieces, commonly named leftand right-handed. Defining the projectors PLand PRas PL=1−γ5 2PR=1+γ5 2,(1.25) it is possible to partition the full flavor space, since introducing the projected quark and antiquark fields ψL,R ≡PL,Rψ¯ ψL,R ≡¯ ψPR,L,(1.26) allows to express the quark (antiquark) field as the sum of ψLand ψR, separating in this way the action into several components: SQCD-F =Zd4x¯ ψL(iγµDµ)ψL+¯ ψR(iγµDµ)ψR+¯ ψLmψR+¯ ψRmψL. (1.27) From the above expression it trivially follows that, in the chiral limit, only the first two components survive. Remarkably, they are completely decoupled in this limit, interacting only through the mass terms in the above formula. For this reasons, the chiral symmetry group (1.24) is often reformulated in terms of its leftand right-handed degrees of freedom, i.e., as SU(Nf)L×SU(Nf)R×U(1)V.(1.28) In any case, when considering quarks of finite but degenerate mass, the axial sector symmetry is broken explicitly—just the SU(Nf) part, since the one corresponding to U(1)Awas already broken by the anomaly—and the above group gets reduced to its vector sector U(Nf)V, or following the structure of (1.24) SU(Nf)V×U(1)V.(1.29) 22 CHAPTER 1. THE LATTICE APPROACH Eventually, if one considers Nfflavors with different masses, the above symmetry gets reduced to the tensor product of Nfcopies of U(1)V. Before closing this interlude, we want to stress a couple of issues concerning chiral symmetry. The first one is that chiral symmetry—or more precisely, its axial sector—is spontaneously broken in QCD at zero temperature. Let us consider the action with just two flavors: up and down. Since these quarks have a very small mass compared to the QCD scale, its explicit symmetry breaking, corresponding to the shift between (1.28) and (1.29) for Nf= 2, should give rise to experimental effects. In particular, some particles—as, e.g., parity partners— should have almost-degenerate masses. However, this expected symmetry does not match the observations; mass differences due to the explicit breaking of the QCD action are significantly smaller than the experimental measures. The origin of these discrepancies is, precisely, that the global symmetry (1.28) is spontaneously broken, since the vacuum of the theory is not invariant under the corresponding transformations. In other words, although the action is almost preserving chiral symmetry, the ground state of the system is not. The second point to discuss is just a corollary of the previous argument: since there is a continuous symmetry being spontaneously broken, Goldstone’s theorem predict the appearance of several massless bosons. This is in fact the case, since if we consider just quarks uand d, we should expect three light bosons— since the symmetry is slightly explicitly broken—which can be identified with the three pions π0,π±. Including also the strange quark s, would account for eight bosons, which can also be identified with the four kaons K0,¯ K0, K±and the η meson, which sums up to the three pions. Note that if we overlook the anomaly, in the last case 9 light mesons should appear, and the η0particle should be also considered within this picture. It is precisely the anomaly what allows to explain the η0−π0mass difference, as was stressed in the previous subsection. With these considerations, we can now proceed to talk of fermions in the lattice, with a last word of caution that summarizes the previous argument. The lightest particles of the QCD spectrum gain their low masses thanks to the spontaneous breaking of chiral symmetry. If we were to tweak the theory and remove the chiral invariance from the action, the explicit breaking of this symmetry would prevent these particles from being pseudo Goldstone bosons, and consequently their masses would be expected to increase. 1.3.4 Dynamical fermions: doubling and chiral symmetry In order to implement full QCD in the lattice, we need to choose an action Slatt which, as in (1.12), can be decomposed into two components. First, a purely gluonic term Sg, for which a viable candidate has already been reviewed in subsection 1.3.2. In second place, a term involving quark and antiquark fields is needed. Now 1.3. THE FORMALISM 23 that the basics of chiral symmetry in continuum QCD have been discussed, we can face the task of constructing this fermionic term, completing in this way the lattice action. Thereby, in this subsection we will analyze the difficulties in proposing a proper lattice version of the fermionic matrix M, which defines how quarks and antiquarks interact. While the construction of the Wilson gauge action (1.17) was not very troublesome, the situation changes substantially when trying to proceed with the fermionic part of the action in a similar way. The simplest approach to discretize the continuum term consists in replacing the first order derivatives in LQCD by a standard symmetric difference. In this way, the continuum fermionic term ¯ ψ(iγµDµ+m)ψbecomes 1 2a¯ ψ(n)iγµUµ(n)ψ(n+ ˆµ)−U† µ(n−ˆµ)ψ(n−ˆµ)+mψ(n),(1.30) which is indeed gauge-invariant, as bilocal products of quark-antiquark fields are connected by the corresponding link variables. Unfortunately, this procedure introduces the so-called doublers: for each of these naive fermions included in the lattice action, 16 copies appear in the continuum limit, 15 of them being unphysical. In order to fix this situation, an alternative is to use Wilson fermions, which add an extra term to the naive formulation. Its aim is to assign a divergent mass, of order O(a−1), to each of the unwanted doublers, decoupling them from the theory in the continuum limit. This mechanism is sufficient to fix the continuum limit of the action, however, the axial sector of chiral symmetry gets explicitly broken in the lattice by the extra term—not only the anomalous part U(1)A, which in fact should be broken, but the whole group U(Nf)A. In continuum QCD, as it was discussed in the previous subsection, chiral symmetry plays a central role, since its non-anomalous sector is spontaneously broken, producing approximate Goldstone bosons. With an explicitly broken chiral symmetry, this mechanism is not possible anymore: particles such as pions, kaons or the η, which conform the so-called meson octet, acquire masses well above its original values. This behavior hinders actual simulations from reaching the physical point. In fact, during decades of lattice numerical works, the achieved mass of the pion—a very relevant quantity, being the lightest hadron in the spectrum—has been very far from its physical value, which makes necessary to perform extrapolations to the physical point. The previous discussion rises a natural question: is it possible to find a fermion formulation that, whithout introducing unphysical multiplicities, preserves chiral symmetry? For a long time, it was believed that, unfortunately, it was not. In fact, far from being a technical complication, the doublers issue has its root in the chiral anomaly of QCD. In the naive lattice formulation of (1.30) the anomaly disappears, canceled exactly by the extra doublers. Adding the Wilson term, which decouples the extra particles, can be interpreted as introducing a lattice version of 24 CHAPTER 1. THE LATTICE APPROACH the anomaly by hand—with one undesirable effect: it also breaks the non-anomalus part of the symmetry, spoiling the spontaneous breaking mechanism. In fact, as was proven in 1981 by Karsten and Smit [67], this is a general result: either the axial anomaly is canceled by the presence of extra fermions, or it is introduced in the lattice by a term that, necessarily, explicitly breaks chiral symmetry, including its non-anomalous sector. In the same way, a contemporary result due to Nielsen and Ninomiya—the so-called no-go theorem—states that it is not possible to construct a lattice formulation of QCD that has at the same time absence of doublers, chiral symmetry, and locality [68, 69]. Furthermore, there also exists a scheme-independent version of this result, with slightly different conditions [70], that makes the impossibility of regularizing a theory with chiral fermions a rather profound question, not exclusive of the lattice approach. As a consequence, a choice must be made between preserving chiral symmetry in its continuum form or avoiding unwanted degrees of freedom; every fermion formulation on the lattice suffers from one pathology or the other. But, even though the above results are correct, a workaround fortunately exists. Almost 40 years ago, Ginsparg and Wilson realized that a remnant of chiral symmetry could be identified within the lattice formulation, in such a way that the axial anomaly is still properly preserved and, at the same time, no unphysical degrees of freedom are introduced [71]. While the continuum form of chiral symmetry—i.e., the preservation of the transformation sets (1.19) to (1.22)—requires the massles Dirac operator to anticommute with γ5, the remnant symmetry proposed in [71] verifies a broader restriction, namely Dγ5+γ5D=aDγ5D, (1.31) where D≡iγµDµstands for the massless Dirac operator. This condition is in fact different from its continuum counterpart for every finite lattice spacing a, although they converge in the a→0 limit. In this sense, (1.31) can be interpreted as an extension of the continuum definition, with a singular advantage: it is able to evade the Nielsen-Ninomiya no-go theorem, since only the continuum form of chiral symmetry—with a vanishing right-hand side in (1.31)—is affected by it [72]. Moreover, a modified chiral rotation can be defined in the lattice, in such a way that the gauge fields transform essentially as its continuum versions, preserving relevant results as, i.e., the index theorem [73]. However, the solution given in [71] was not constructive, in the sense that a particular form of Dcould not be found at that moment; it would took almost two decades to find a practical implementation of these ideas. Within the above scenario, other alternatives would be explored over the years to come. In fact, a number of these strategies were developed and nowadays are part of the current lore of the field. They can be classified in groups according to 1.3. THE FORMALISM 25 how they deal with the doubling problem. Some approaches choose to give up chiral symmetry; this is the case of Wilson fermions. In this category are also included clover fermions, a Symanzik improved version of Wilson fermions that removes O(a) discretization errors [74], and twisted mass fermions, which consider pairs of mass degenerate Wilson fermions together with an isospin mass splitting term [75]. An alternative approach consists in preserving chiral symmetry, while assuming some of the doublers degeneracy. Here the prime example are Kogut-Susskind or staggered fermions [76]. This variant is constructed from the naive formulation, but distributes the usual Dirac 4-spinor components over neighboring sites, in such a way that the doublers degeneracy gets reduced from 16 to 4 species. Moreover, for some observables it is possible to remove the residual degrees of freedom by taking the fourth root of the fermion determinant—a procedure commonly referred to as rooting, originally proposed by Marinari, Parisi and Rebbi, in the context of the massive Schwinger model [77]. Although its validity was at first controversial, both numerical evidence [78] and theoretical arguments [79] support that this technique leads to the correct theory as long as continuum and chiral limits are taken precisely in this order. As with Wilson fermions, there also exist Symanzik improved versions of the staggered action; two of them, widely used by the lattice community, are the Asqtad action [80]—for a-squared tadpole improved—and the HISQ action [81]—for highly improved staggered quarks. In Chapter 4 we make use of staggered quarks, while on Chapter 6 the analyzed configurations were generated with the HISQ action. For the sake of completeness—even if it is less related with the work developed in this thesis—there is still another class of fermions that deserves at least a mention: the particular solutions of the Ginsparg-Wilson equation (1.31). They preserve chiral symmetry without the doublers pathology, so in this sense they are the best possible fermions, and they should be preferred when compared to other alternatives. However, all their implementations suffer from the same illness—they are by far the most expensive fermions in computational terms. As a consequence, its use is reserved to situations where chiral symmetry is needed with considerable precision. There exist two different solutions that are broadly used: they are called overlap [82–84] and domain-wall [85,86] fermions. In the first case, the fermionic matrix is constructed by operating with a Wilson-like matrix—the result being a costly non-sparse matrix. For domain-wall fermions, an extra fifth dimension of infinite extent is needed to preserve chiral symmetry. Since in practice this is implemented as an additional dimension in a finite lattice, a mild violation of the symmetry still remains, although it can be controlled. 26 CHAPTER 1. THE LATTICE APPROACH 1.3.5 Monte Carlo: ensembles and observables Summing up all the previous considerations, we have now a fully regularized version of QCD. Thereby, the original expressions for the partition function (1.2) and the vacuum expectation value (1.3) of a given observable acquire now a well-defined mathematical meaning. If we now want to consider, e.g., a purely gluonic observable O(U), we should compute the following integral: hOi =RDUO(U)e−Slatt(U) RfflathcalDUe−Slatt(U).(1.32) The above integration measure DUis composed by the product of as many individual Haar measures as links are in the given lattice, which amount to V≡QiLi times the dimension of the lattice. So, for almost any size that we can think of, (1.32) contains a highly multidimensional integral, with a vast configuration space that, even for small lattices, can not be fully explored by any computational means. There exist however a class of numerical approaches, commonly referred to as Monte Carlo methods, that are capable of properly sampling these kind of spaces—to apply them to the lattice regularization was, in fact, one of the key proposals of Wilson [40]. In order to compute the expectation value of (1.32), a first Monte Carlo approach could consist in generating a random sample of link configurations, i.e., a collection or ensemble of gauge fields Ui, each one containing the information of each individual link Uµ(n) in the lattice. Then, hOi could be computed as a weighted average over the previous ensemble, with weight e−Slatt(U). In a similar fashion, statistical errors could be computed by standard techniques. However, this approach is too naive for the problem at hand, since the exponential factor in (1.32) highly suppresses a great majority of the configuration space elements; in other words, only a small subset of the whole space contributes significantly, and a naive random sampling of the space would miss the area of interest. In order to fix this issue, importance sampling is to be applied: again, an ensemble of configurations Uiis generated, but they should be selected with probability e−Slatt(Ui). Then, provided a well-distributed ensemble of Nelements is available, an estimator for the expectation value of a given observable can be straightforwardly computed as ¯ O=1 N N X iO(Ui).(1.33) The original problem is now shifted to the generation of the Uicollection. This can be achieved by starting from a given initial configuration, say U1, and following a Markov chain process, in which the selection of the next configuration is regulated 1.3. THE FORMALISM 27 by a probability that depends only on the immediate previous state. The idea is that, even if one starts far from the meaningful configurations, the Markov process should drive the system to an equilibrium state, in which the distribution e−Slatt(Ui) is reproduced. To guarantee that, the transition probabilities from each Uito any other Ujneed to be adequately defined. It is sufficient (but not necessary) to require the detailed balance condition to be fulfilled, so if we label the transition probability of the chain as p(Ui|Uj), p(Ui|Uj)e−Slatt(Uj)=p(Uj|Ui)e−Slatt(Ui)(1.34) would be required for every iand j. In this case, if the system is also ergodic— meaning every configuration is accessible from any other in a finite number of steps—it is guaranteed that, with independence of the initial configuration chosen, the Markov process will reach equilibrium and, consequently, a well-distributed ensemble will be generated. Then, to elaborate a particular algorithm that follows the described Markov process consists then in specifying how the next element of the chain is selected, and which transition probabilities connect the configuration space. With respect to the latter point, most of the algorithms used by the lattice community verify the detailed balance condition. Maybe the most elementary of them is Metropolis algorithm, which accepts a change between Uiand Ujwith probability max {1, e−∆Slatt }. This is exactly the approach followed in Chapter 3. More involved alternatives include the heatbath algorithm [44], which often includes some overrelaxation steps [87–89]. Furthermore, for the more general case in which fermions are included into the action, one of the most widely used algorithms is the Hybrid Monte Carlo [90,91], which at each update combines a microcanonical evolution with a final metropolis step to accept or reject the proposed change. In any case, the stochastic nature of Monte Carlo methods introduce some uncertainty into the computed observables, in the form of statistical errors associated to the corresponding estimators. Nevertheless, these can be dealt with without much difficulty, just by taking into account two fundamentals. First, that the generated ensemble has to be thermalized, meaning that sufficient iterations need to be spent to reach the equilibrium distribution of the Markov chain. This can be achieved by monitoring a set of observables that allow to determine when the Markov evolution is stationary. Secondly, and even more important, is the inherent correlation of the configurations generated by the update process—generally a local algorithm that needs a high number of iterations to produce a new configuration far from the original. In this case, the use of standard error analysis tools, such as jacknife binning, or even the direct computation of autocorrelation times for each observable, allows to estimate in a reliable way the statistical errors of any computed observable. 28 CHAPTER 1. THE LATTICE APPROACH 1.4 Reaching the continuum and sources of error Up to this point, we have covered how continuum QCD can be regularized into a spacetime lattice where, taking advantage of Monte Carlo methods, it is possible to estimate vacuum expectation values of given observables. The computed observables will depend in general on the dimensions of the lattice and on the bare values of the gauge coupling βand the Nfmasses of the fermion species involved. Although knowing these values can be enough for some applications, in most calculations the objective is to obtain a physically measurable quantity that can be eventually contrasted with experiment. This is not the case of the present thesis, since the work being presented is not directly concerned about the last part of a general lattice calculation; however, in order to give a general perspective of what a complete computation would require, and for the sake of completeness, hereafter we proceed with a brief overview of the issues involving the last steps of the lattice approach. Obtaining a physical prediction from bare lattice results requires a somewhat involved limit procedure; not only the lattice spacing ashould vanish—this would be the continuum limit—but also the bare couplings should go to their physical values. In other words, the physical point needs to be reached; if not, the simulated theory would be describing an alternative scenario, one with, e.g., different hadron masses. In fact, the couplings of the theory are not directly observable (since confinement prevents quark fields to manifest outside a color-singlet, the mass of the quark is not an observable and depends on the renormalization scheme). In place, other observables, such as hadron masses ratios—which of course depend on the aforementioned parameters—are to be computed in the lattice, so they can be compared with the experimental values. Thus, in order to compute a given quantity, it is necessary to keep track of several additional observables, one for each of the free parameters—or flavor species—that are included into the calculation. Morover, it should be noted that the lattice spacing ais not a free parameter of the lattice theory, but another quantity that needs to be measured within the lattice. The above scenario can be summarized as follows: we need to take the continuum limit, a→0, while driving a set of physically measurable quantities to their experimentally determined values. Moreover, since doing the computations on a lattice of infinite extent is out of reach, the procedure is to be performed on a 4D box of finite physical size.5In order to take the a→0 limit, it is required to repeat the computations at different lattice sizes, keeping the physical size of the box constant. In this way, lattices with different values for the lattice spacing a—and consequently for the number of nodes N—are computed, while the volume of the 5In fact, repeating the calculations in several box sizes allows to extrapolate the results to the thermodynamic limit. 1.4. REACHING THE CONTINUUM AND SOURCES OF ERROR 29 box a4Nis kept constant. It should be recalled that while Nis a free parameter, ais an observable that needs to be measured (depending not only on the gauge coupling but also on the inclusion of fermionic species). So, as a first consequence, keeping the physical volume constant is not a completely trivial issue. The same problematic that appears when trying to fix the volume of the 4D box applies to every physical observable that is intended to be kept constant (or driven to a suitable value, obtained from experimental sources). To achieve this, essentially two options are possible: either the free parameters are tuned in an iterative process—one lattice at a time up to the physical point—or the computations are repeated for several values of these parameters, and then interpolated to the desired point. A less ideal variation of the latter is to reach the physical point by means of an extrapolation. This was in fact the rule during several decades of lattice numerical computations, since the different types of fermion discretizations (reviewed in Section 1.3.4) have many difficulties in reaching a sufficiently low mass for the pion.6 As a final remark, it is important to stress that any full computation on the lattice should include a careful analysis of all sources of errors. First, statistical errors are introduced by the Monte Carlo evaluation process. These are however relatively easy to deal with, since they can be estimated by standard methods and, even more importantly, can be systematically reduced by increasing the computation time. A more challenging obstacle is in estimating the different issues leading to systematic errors, specially when a high precision is desired—and previously unnoticed systematic errors can become big enough to be taken into account. A non-exhaustive but common record of the different sources of systematics considered would include the errors associated with the tuning or interpolation to the physical point and the continuum limit process. Usually less severe are those related with the thermodynamic limit—finite volume effects—or with the exclusion from the action of electromagnetic effects, which, given that the quarks are electrically charged particles, should be accounted for. In general, any approximation assumed should account for its corresponding error—as, e.g., the case of some lattice fermion types which consider uand dquark bare masses as degenerate, and thus are required to account for the corresponding isospin breaking effects. 6Although Ginsparg-Wilson fermions preserve chiral symmetry, they suffer from a large computational overhead that spoils their advantage with respect to other fermion types, when reaching a physical pion mass is considered. 36 CHAPTER 2. TOPOLOGY IN QCD pieces, Lθ-QCD =LQCD +iθ g2 64π2µνρσFa µνFa ρσ,(2.8) where LQCD stands for the standard expression (1.4), which is both CP conserving and real. In contrast, the second piece in the right hand side of (2.8) breaks CP— as was stressed in the previous section—and, more importantly for the upcoming discussion, is a purely imaginary number. The topological nature of the θterm has direct physical effects: the vacuum of the theory cannot be constructed as a quantum fluctuation around a classical definite state [97]. Therefore, in order to properly analyze in which way QCD depends on the θparameter, a non-perturbative approach is required. Under this circumstances, it would seem adequate to deal with the inclusion of the extra term Lθwithin the lattice regularization formalism. In principle, it would suffice to find a proper discretization for the topological charge Q, and include its effects in the corresponding Monte Carlo algorithm. Inconveniently, this is far from being enough if non-vanishing values of the vacuum angle are to be explored, since the θterm amounts to a complex phase in the action that gives rise to a severe sign problem. Usual importance sampling methods are of no use in this situation, and workarounds need to be found. At present, much of the achieved progress in this topic involves computations of topological quantities at θ= 0, where standard MC methods are still applicable. Although the definition of the topological charge on the lattice is a subtle question, it is possible to compute quantities such as the topological susceptibility χ, even with some difficulties [102]. However, little progress has been achieved in the last decades in what concerns the study of the θ > 0 case. It is worth noting that the inclusion of the θterm is not the only example of a sign problem being induced by a complex component into the action of QCD. Another major example is given by finite density QCD, which takes place when a non-zero chemical potential term is included into the fermionic matrix—a required addition when considering high baryonic densities. This, in fact, is a relevant scenario which affects several areas. It is needed in the studies of astrophysical objects such as neutron stars. Moreover, in early Universe investigations, dealing with extremely high values of both temperature and density is required. And, last but not least, particle collider facilities reproduce these conditions in heavy ion experiments. Although, thanks to asymptotic freedom, perturbative expansions can offer an insight for some limiting cases—namely high energy, which translates in high Tor high µ—and the µ= 0 case is accessible to the standard importance sampling techniques, almost the full µ−Tphase diagram is beyond the reach of these approaches. In both of the cases above, considering a complex action results in the appearance of a sign problem, which, as it was warned at the beginning of this thesis, 2.3. TRYING TO OVERCOME THE SIGN PROBLEM 37 constitutes one of the Millenium problems—in particular one that is not expected to be solved with a positive outcome, which for our interests would be P=NP. Therefore, the community has directed their efforts towards the development of different alternatives that try to evade, or in some cases ameliorate, the sign problem in the systems of interest, i.e., in QCD or QCD-like models. In the next section, we mention some of the more popular strategies, including a brief review of the methods developed by Azcoiti et al [12, 13], which in fact are applied in Chapter 4, when reconstructing the θdependence of the Schwinger model. 2.3 Trying to overcome the Sign Problem In order to progress in the study of complex action systems, different methods have been developed over the years. In the cases where the sign problem is mild enough, such as in finite density QCD for small values of the chemical potential µ, several strategies can be followed with success. Perhaps the more straightforward approach is the technique known as reweighting. It introduces an auxiliary partition function with non-negative density, which is used to reformulate the expectation value of a given observable—computed in the original ensemble—in terms of other expected values that are to be determined within the auxiliary ensemble. Indeed, this change does not eliminate the problem, since the computational cost of this method escalates exponentially with the volume of the system; however, it is useful when applied to small systems, or in certain regions of the parameter space, such as QCD with a small chemical potential µ, where the effects of the sign problem are far from being severe. In this scenario, alternative approaches include Taylor expansions around µ= 0 and analytic continuations from purely imaginary chemical potential, although its scope is heavily bounded by the presence of a SSP; for an extended discussion on this topic, we recommend the reader [103,104]. Apart from the above methods, that provide some insight for mild sign problem regions—but fail to deliver otherwise—other procedures that have, in principle, greater scope, have been developed over the last decades. In chronological order, the first is Complex Langevin. The original works of Parisi and Klauder proposed to generalize the Langevin equation formalism to the case of a complex-valued distribution [9,10]. Even when rigorous proofs were lacking—e.g., the existence of a stationary solution was a conjecture—the technique seemed to give correct results in some cases [105]. However, in other systems the algorithm failed to converge, or it did but to the wrong limit [106]. Interest in the method was rebounded when Berges and Stamatescu realized that the instabilities of the Langevin evolution can be dealt with if the stepsize is reduced enough [11]. In fact, using an adaptative stepsize eliminates this problem [107]. Notwithstanding its recent successes, the approach still has some caveats that need to be addressed, since in some scenarios 38 CHAPTER 2. TOPOLOGY IN QCD it continues to converge to a wrong result. Current efforts are devoted to identify under which conditions the proper solution can be reached, and what mechanisms can be used to monitor if the Langevin evolution is to be trusted; for a review on this topic, we refer to [108]. With a special focus on systems with a topological term in the action, a different strategy was proposed by Azcoiti et al in 2002 [12]. In summary, the topological charge dependence on θis reconstructed from the probability distribution function at θ= 0, which is to be computed from simulations at purely imaginary values of θ, possible since in this case the action becomes real. These results are to be adjusted to a suitable analytical expression that, once integrated, allows to obtain q(θ) with the help of multiprecision algorithms. This approach proved to work well in a number of systems, including the one dimensional Ising model and the U(1) compact model in two dimensions, giving predictions for CP3and, in a later work, for CP9—a model that, as QCD, exhibits confinement and asymptotic freedom [109]. Nevertheless, the scope of the approach gets hampered by the fact that it cannot reconstruct a non-monotone order parameter [21]. In other words, if a given model breaking CP gets this symmetry restored at θ=π, then q(θ) needs to decrease at some point; in this situation, the method fails to reproduce this behavior and a flattening is observed—which, on the other hand, it is generally avoided in models with (spontaneously) broken symmetry at θ=π. Initially with the aim of crosschecking the above method, an alternative approach, using the same input—i.e., Monte Carlo simulations at imaginary values of the vacuum angle θ—was proposed in 2003 [13]. In this case, an additional assumption needs to be made, namely that in order to reconstruct the full θdependence, no critical points are allowed, except at most at θ=π. This poses a severe obstacle to its applicability in finite density QCD, where a rich phase diagram is expected, but is in principle well adapted to θQCD, where only one phase transition—if any—is surmised, precisely at θ=π[102]. Besides the later assumption, a particular extrapolation is required, that needs certain observables to vary as slowly as possible. For this reason, the method is expected to work well, among others, in asymptotically free gauge theories. This approach, which has already obtained good results in a variety of models [109–111], has been applied in Chapter 4 to reconstruct the θdependence of the massive Schwinger model; for a more extensive review, covering the specifics of the actual implementation, we refer the reader to Section 4.3. Finally, other recent approaches to complex action systems can be mentioned, such as Lefschetz thimbles [14–16] or the density-of-states or LLR method [17–19]. Their details are beyond the scope of this thesis, since their application to QCD with a topological term seems, for now, remote. Chapter 3 Ising model with a θterm In this chapter we present our work on the two-dimensional antiferromagnetic Ising model with a purely imaginary magnetic field, which can be interpreted as a toy model for the usual θphysics, and that was published in [28]. Our motivation, as it was anticipated at the beginning of this thesis, is twofold. First, we pretend to provide a benchmark calculation in a system which suffers from a strong sign problem, so that our results can be used to test Monte Carlo methods developed to tackle such problems. In second place, we want to test the predictions of the method developed by Azcoiti et al. [13] regarding this system [21], since its performance with the expected non-trivial phase diagram could entail additional obstacles for its reliable application in other scenarios. In the forecoming sections, we justify the choice of model and review their fundamentals. Then, we discuss the analytical techniques applied, including the exact computation of the first eight cumulants of the expansion of the effective Hamiltonian in powers of the inverse temperature, which allows to calculate physical observables for a large number of degrees of freedom with the help of standard multi-precision algorithms. Finally, we report accurate results for the free energy density, internal energy, standard and staggered magnetization, and the position and nature of the critical line, which confirm the mean-field qualitative picture of [21], and which should be quantitatively reliable, at least in the hightemperature regime, including the entire critical line. 3.1 Why Ising? As has been argued along the first chapters of this work, numerical simulation of systems with a severe sign problem is one of the major challenges for high-energy theorists—a statement which is also valid for their solid-state colleagues. If we denote the microscopic states of a given physical system by s, and the thermody39 40 CHAPTER 3. ISING MODEL WITH A θTERM namics of such system is described by a partition function of the form Z=X s P(s),(3.1) we say that the system in question presents a sign problem if the “weights” P(s) are not real and positive: This implies that we cannot interpret P(s) as a proper probability distribution, and the standard, efficient Monte Carlo algorithms cannot be applied. Not all sign problems are equally severe. Let us restrict ourselves for simplicity to the case where the P(s) are real but not positive definite1. One can easily devise a reweighting algorithm that uses the absolute value |P(s)|as the weight of each state, and shifts the sign of P(s) into the observables. Now a standard Monte Carlo method is applicable, and in the limit of infinite statistics we should obtain the correct result. With finite statistics, however, a key quantity is the thermodynamic average of the sign of each contribution to the partition function, that is, hsign(P(s))i. If this quantity goes to zero exponentially with the volume, hsigni ∝ e−αV , then we would need an exponential amount (in the volume of the system V) of statistics to get correct results, which is of course impossible in practice. In this case we say that the sign problem is severe. Beyond QCD at finite baryon density or QCD with a topological term in the action, there exist other physically relevant systems which suffer from a SSP. Some of the most popular examples include chains of quantum spins with antiferromagnetic interactions, the two-dimensional O(3) non linear sigma model with a topological term or the Hubbard model. The existence of a SSP is the main reason for the little progress made on the theoretical understanding of these physical systems outside of phenomenological models. In order to check novel Monte Carlo methods designed to tackle such problems, it is highly desirable to have a set of benchmark calculations as extensive as possible. For very few systems an analytic solution is known, for example, the one-dimensional antiferromagnetic Ising model with an imaginary magnetic field, the two-dimensional compact U(1) model with topological term, or the twodimensional Ising model with an imaginary magnetic field h=iπ/2. In a few other cases the sign problem can be avoided by reformulating the physical system with new degrees of freedom, taking advantage of the fact that a good choice of these degrees of freedom provides an equivalent physical system free from the sign problem, which can therefore be simulated by standard methods; see e.g. [112] for a recent discussion on this dualization approach. Unfortunately this idea works only in a few cases which, until now, are not the most interesting physical systems—indeed none of the examples previously mentioned have been solved with this idea. Within the above scenario, our intention in this work is to provide a benchmark calculation for a system for which we do not have an analytic solution available, 1The discussion for complex weights does not add any fundamental difficulty. 3.1. WHY ISING? 41 nor a reformulation that avoids the sign problem. We study the two-dimensional antiferromagnetic Ising model with a purely imaginary magnetic field, which can be thought of as a toy model for the usual θphysics. Indeed the Euclidean partition function for QCD with a nonvanishing θterm can be written in the form ZV(θ) = X n pV(n)eiθn (3.2) where n, the topological charge, is an integer, and pV(n) is, up to a normalization, the probability of the topological sector nat θ= 0. This has the same structure as the partition function of the antiferromagnetic Ising model in an external purely imaginary magnetic field, as we will see in detail later on, and we expect that the SSP in both systems should also be similar. This system was studied in [20] by locating the zeros of the partition function in the complex temperature-magnetic field plane, and they found, for purely imaginary magnetic field, a rich phase structure with two phases characterized by a vanishing (nonvanishing) staggered magnetization, separated by a phase transition line. We study this system by an exact cumulant expansion to eighth order, followed by the analytic computation of the partition function and other physical quantities for a large number of degrees of freedom with the help of a standard multiprecision algorithm. This amounts essentially to the computation of the effective Hamiltonian up to order T−8, and therefore is expected to work well in the high-temperature regime, and we provide strong evidence that this is indeed the case. Our results are consistent with [20], and extend the results of [21], obtained through the application of algorithms developed in [12,13], and through a meanfield analysis. We are able to obtain a more precise quantitative determination of the transition line separating the paramagnetic and antiferromagnetic phases of the model. For some systems with a SSP, we know a priori that the partition function will be positive, for example systems in thermal equilibrium with a (Hermitian) Hamiltonian description. Such is the case in a quantum field theory with a θterm. In the toy model we study here, although we do not have a rigorous proof in this case,2we have evidence that, at least in the region where the approximation we use is valid, the partition function is indeed positive (it is trivially always real). Such evidence is twofold. First, we can prove rigorously that up to the fifth cumulant, the partition function is indeed positive. Unfortunately we have not been able to extend this proof to higher cumulants, but in our multiprecision calculations with up to eight cumulants, we have never seen an instance where 2This would imply a nontrivial restriction on the position of the Lee-Yang zeros for the antiferromagnetic Ising model. To the best of our knowledge, very little is rigorously known about such zeros. 42 CHAPTER 3. ISING MODEL WITH A θTERM the partition function is negative or vanishes. This is highly nontrivial: If instead of a constant imaginary magnetic field we try, for example, to put a staggered imaginary field in our lattice (this is of course equivalent to the ferromagnetic model with a constant imaginary field), we immediately get a fluctuating sign for the partition function. Second, there have been studies locating the Lee-Yang zeros of the antiferromagnetic two-dimensional Ising model up to 142lattices [113], and in 12 ×13 lattices [20]. Up to that size there is no sign of any zeros cutting the imaginary axis at any temperature. Whereas this by no means amounts to a rigorous proof, we believe it provides a strong indication that, at least in the region of interest for our work, this model should have a positive partition function. Hereafter, Section 3.2 is devoted to formulate the model and to recall the main ingredients and results of the mean-field approximation developed in [21]. In Sec. 3.3 we introduce the cumulant expansion, report the analytical results for the first eight cumulants in the two-dimensional model, and write the analytical expressions for the free energy and mean values of interesting physical quantities. The results for the staggered magnetization, susceptibility, and phase diagram of the model are reported in Sec. 3.4, where we also compare our results at h= 0 and iπ/2 with the analytical solutions of [114–116]. In Sec. 3.5 we report our conclusions. The technical details of the analytical computation of the cumulant expansion can be found in Appendix A. 3.2 Two-dimensional Ising model The Ising model [20, 114–119] has been studied for a long time now, and it has known analytical solutions in the one-dimensional case at any external magnetic field h[117], and in two dimensions only for the case without magnetic field h[114] and for h=iθ/2 = iπ/2 [115, 116]. The model with a pure imaginary magnetic field suffers from a SSP in any number of dimensions. In addition to that, the expected phase diagram for d≥2 is non trivial [21], making the reconstruction of the θdependence of the observables even more challenging. All this makes the model a good theoretical laboratory to test new methods designed to deal with the SSP. It is therefore worthwhile to carry out a detailed study of this model at purely imaginary magnetic field, particularly because little progress has been achieved on reconstructing the θdependence of the observables, apart from the analysis of [21] and the recent study in [120]. 3.2. TWO-DIMENSIONAL ISING MODEL 43 The partition function of the model, following the conventions of [21], is: Z=X {si} exp FX <ij> sisj+iθ1 2X i si!.(3.3) The half magnetization M 2≡1 2X i si,(3.4) is an integer taking any value between −N/2 and N/2, where Nis an even number denoting the total number of spins in the lattice. It is in this sense that we identify M/2 with a topological charge and regard the imaginary magnetic field term in the action as a θterm. It is important to mention that, from now on, we will consider only the antiferromagnetic case F < 0, since the model with imaginary field does not define a unitary theory for arbitrary values of the ferromagnetic coupling [115,121]. As we shall see in detail in the next section, by dividing the rectangular lattice into two sublattices, introducing the respective magnetizations M1and M2, making a cumulant expansion and keeping only the first cumulant, we arrive at the following approximation to the partition function (where ddenotes the dimensionality of the lattice): Z1c(F, θ) = X {si} exp iθM1+M2 2+ 4Fd NM1M2.(3.5) We recall now the mean-field analysis carried out in [21]. The resulting partition function, ZMF (F, θ) = X {si} exp iθM1+M2 2−Fd N(M1−M2)2,(3.6) is different from Eq. (3.5). However, it can be seen to give the same qualitative results for the observables and the phase diagram. In this regard, we will consider the first-cumulant expansion Z1cas a mean-field approximation to Z, and the general expansion itself as an improvement of it, at least for small F, where the expansion is expected to converge. Applying standard saddle-point techniques to the mean-field partition function [21], one obtains the F−θphase diagram shown in Fig. 3.1. A second order critical line, dFc=1 2cos2θc 2,(3.7) separates two different phases: a staggered one, with hmsi 6= 0, for F > Fc(θ), and a paramagnetic one, with hmsi= 0, for F≤Fc(θ). 44 CHAPTER 3. ISING MODEL WITH A θTERM 0 0.1 0.2 0.3 0.4 0.5 0π/2π hmsi= 0 hmsi 6= 0 d|F| θ dFc(θ) Figure 3.1: Phase diagram of the mean-field approach of [21] to the antiferromagnetic Ising model in the F−θplane. 3.3 Cumulant expansion and observables Our interest is focused on the antiferromagnetic model, where the staggered magnetization is a good order parameter. From now on we will work with a rectangular two-dimensional lattice, although the method is easily generalizable to any number of dimensions. We divide the lattice into two sublattices Ω1and Ω2in a chessboard fashion. In the two-dimensional lattice this means that if iand jindex, respectively, the row and the column of a given spin, this spin will be in the first (second) sublattice if the sum i+jis even (odd). For simplicity we will require both lengths of the lattice to be even. Denoting by Nthe total number of points in the lattice, we define the magnetization densities m1and m2as mj≡Mj N/2≡Pi∈Ωjsi N/2j= 1,2,(3.8) and the density of staggered magnetization is ms≡m1−m2 2.(3.9) Let us denote by g(m1, m2) the number of microstates with magnetization 3.3. CUMULANT EXPANSION AND OBSERVABLES 45 densities m1and m2in sublattices Ω1and Ω2, respectively, that is, g(m1, m2) = X {si} δ X i∈Ω1 si−M1!δ X i∈Ω2 si−M2!.(3.10) A trivial computation gives: g(m1, m2) = N/2 N1+N/2 N2+,(3.11) with Nj+≡N(1 + mj)/4 for j= 1,2. Now, by restricting ourselves to the set of configurations with given magnetization densities m1and m2, it is straightforward to define the expectation value of a general observable O({si}) within this subset— i.e., at fixed m1, m2—as: hOim1,m2≡1 g(m1, m2)X {si} δ(X i∈Ω1 si−M1)δ(X i∈Ω2 si−M2)O({si}).(3.12) Then, the sum over all possible spin configurations in the original partition function (3.3) can be partially summed up—at least formally—grouping together sectors with equal magnetization densities m1, m2. By doing so, and taking into account the above definitions, the reformulated partition function takes the following form: Z=X m1,m2 g(m1, m2)*exp iθ 2X i si+FX <ij> sisj!+m1,m2 .(3.13) The θterm in Eq. (3.13) is just iθ (m1+m2)N/4, and therefore constant at fixed m1and m2; we can take it out of the expectation value, arriving at Z=X m1,m2 g(m1, m2)e1 4Niθ(m1+m2)*exp FX <ij> sisj!+m1,m2 .(3.14) We cannot evaluate exactly the expectation value in Eq. (3.14), as that would be equivalent to solving exactly the model for arbitrary values of the external field. Instead we perform a cumulant expansion and truncate at a given order. Let us recall the definition: etX≡exp ∞ X n=1 κn tn n!!,(3.15) where the nth cumulant κnis an nth degree polynomial in the first nnoncentral moments of X, given by the following recursion formula: κn=µ0 n− n−1 X m=1 n−1 m−1κmµ0 n−m, µ0 n≡ hXni.(3.16) 52 CHAPTER 3. ISING MODEL WITH A θTERM 0.5 0.6 0.7 0.8 0.9 0 0.1 0.2 0.3 0.4 0.5 0.6 e|ns (F) |F| k=1 k=4 k=8 Analytic Figure 3.6: Nonsingular part of the internal energy at θ=π, N = 2000. −0.25 −0.2 −0.15 −0.1 −0.05 0.1 0.2 0.3 0.4 0.5 0.6 cv(F) |F| k=1 k=4 k=8 Analytic Figure 3.7: Specific heat at θ=π, N = 2000, plotted against the analytical expression. 3.5. CONCLUSIONS 53 0 0.5 1 1.5 2 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 cv(F) |F| k=1, N=2000 k=4, N=2000 k=8, N=2000 k=8, N=6000 Analytic Figure 3.8: Specific heat at θ= 0, plotted against the analytical solution. At θ= 0, Fc= log(1 + √2)/2≈0.4407. 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0.28 0.29 0.3 0.31 0.32 0.33 0.34 0.35 0.36 hm2 si |F| N = 400 N = 800 N = 1600 N = 3200 Figure 3.9: hm2 sicurves at θ= 2, k = 8. Solid lines are just a guide to the eye. 54 CHAPTER 3. ISING MODEL WITH A θTERM 0 0.5 1 1.5 2 2.5 3 3.5 4 0.24 0.26 0.28 0.3 0.32 0.34 0.36 0.38 0.4 dhm2 si/dθ |F| N = 100 N = 200 N = 400 N = 800 N = 1600 N = 3200 Figure 3.10: Scaling of dhm2 si/dθ at θ= 2, k = 8. Solid lines are a guide to the eye. 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 π/4π/2 3π/4π0 |F| θ k = 1 k = 4 k = 8 Matveev and Shrock, 2008 Figure 3.11: The critical line Fc(θ), computed as the maximum of dhm2 si/dθ at N= 2000. The maximal Fpoints obtained in [20] are also shown. 3.5. CONCLUSIONS 55 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.28 0.29 0.3 0.31 0.32 0.33 0.34 0.35 0.36 cv(F) |F| N = 400 N = 800 N = 1600 N = 3200 Figure 3.12: Specific heat cvwith k= 8 and θ= 2. Solid lines are just a guide to the eye. magnetization as an order parameter. The finite-size scaling suggests that the two phases are separated by a continuous phase transition line. The position of the critical point at θ= 0 is in very good agreement with the exact result Fc= log(1 + √2)/2≈0.4407, and the free and internal energy densities at θ=π agree also well with the analytical prediction, at least in the high-temperature regime, thus giving reliability to our results in this region. Therefore this model could be a good laboratory to check proposals to simulate physical systems afflicted by a SSP. Moreover, we have confirmed that the system under discussion has a more involved phase diagram than in the expected θQCD case, which would have, at most, a single critical point in θ=π. In this sense, the reconstruction method of [13], which when applied in [21] was capable of detecting the presence of a critical region, is strongly supported. 56 CHAPTER 3. ISING MODEL WITH A θTERM Chapter 4 Massive 1-flavor Schwinger model with a θterm We analyze here the massive 1-flavor Schwinger model with a θterm and a quantized topological charge. Our work, published in [29], relies on the approach of Azcoiti et al in [13]. We are able to calculate the full dependence of the order parameter with θin a system that includes dynamical fermions. Moreover, our results at θ=πare compatible with Coleman’s conjecture [22] on the phase diagram of this model. This chapter is organized as follows: after motivating the topic in the first section, we summarize some relevant features of the Schwinger model with a topological term in Sec. 4.2. Since the proposal [13] to analyze physical systems with a topological term in the action has been found to be particularly well suited to bypass the sign problem in asymptotically free gauge theories, we decided to apply it, and Sec. 4.3 contains a brief review of the method. In Sec. 4.4 we give some technical details concerning the lattice setup and the computer simulations performed. Sec. 4.5 shows our results for the topological charge density as a function of θat several fermion masses and gauge couplings, and finally we end this chapter by reporting our conclusions. 4.1 Motivation The nature and origin of dark matter constitute one of the most wide open questions in modern physics. To gain some insight in this puzzling problem, it is highly desirable to elucidate the existence of new low-mass, weakly interacting particles from a theoretical, phenomenological and experimental point of view. As was outlined in Chapter 2, the light particle that has gathered the most attention is the axion, predicted by Weinberg and Wilczek [99], and Wilczek [100] in the Peccei and 57 58 CHAPTER 4. SCHWINGER MODEL WITH A θTERM Quinn mechanism [98] to explain the absence of parity and temporal invariance violations induced by the QCD vacuum. The axion is one of the more interesting candidates to make the dark matter of the universe, and the axion potential, that determines the dynamics of the axion field, plays a fundamental role in this context. The QCD axion model relates the topological susceptibility χTwith the axion mass maand decay constant fathrough the relation χT=m2 af2 a. The axion mass is, on the other hand, an essential ingredient in the calculation of the axion abundance in the Universe. Therefore, a precise computation of the topological properties of QCD and of their temperature dependence becomes of primordial interest in this context. Understanding the role of the θparameter in QCD and its connection with the strong CP problem is one of the major challenges for high energy theorists [122]. The calculation of the topological susceptibility in QCD is already a challenge, but calculating the complete potential requires a strategy to deal with the presence of a highly oscillating term in the path integral; in other words, one needs to circumvent a severe sign problem. In fact euclidean lattice gauge theory, our main non-perturbative tool for studying QCD from first principles, has not been able to help us much because of the imaginary contribution to the action coming from the θterm, that prevents the applicability of the importance sampling method [102]. This is the main reason why the only progress in the analysis of the finite temperature θdependence of the vacuum energy density in pure gauge QCD, outside of approximations, reduces to the computation of the first few coefficients in the expansion of the free energy density in powers of θ[123], and the situation in full QCD with dynamical fermions is, on the other hand, even worse [124–129]. Much experience has been developed in the last years concerning the strengths and weaknesses of the approaches [12,13], which aspire to eventually overcome the sign problem in θQCD. As a matter of fact, it has been applied successfully to the computation of the vacuum energy density and the topological charge density in a handful of interesting physical systems [21,109–111,130]. Our purpose in the present chapter is to take advantage of this experience to perform a first step in the ambitious program of computing the θdependence of the QCD vacuum energy density. Thereby, we analyze the θdependence of a toy model for QCD, the Schwinger model, on the lattice. Strictly speaking, the Schwinger model in the continuum is not asymptotically free, as QCD, since it is super-renormalizable and the CallanSymanzik β-function vanishes. However, in the lattice version, since the continuum coupling is dimensionful, the continuum theory is reached at infinite inverse square gauge coupling β= 1/e2a2, much in the same way as four-dimensional asymptotically free gauge theories such as QCD. Furthermore the model is confining [131], 4.2. THE MASSIVE SCHWINGER MODEL WITH A θTERM 59 exactly solvable at zero fermion mass, has non-trivial topology and shows explicitly the UA(1) axial anomaly [132] through a non-vanishing value of the chiral condensate in the chiral limit, in the one-flavor case. These are basically the reasons why this model has been extensively used as a toy model for QCD. It should be noted that, for two dimensional systems such as the Schwinger model with a θterm, there exist numerical methods such as Hamiltonian methods [133–135] and the Grassmann tensor renormalization group method [136] that have been applied successfully. However, such approaches are currently only applicable to two-dimensional systems, whereas our aim is to test a method that should, in principle, be applicable also to four-dimensional theories such as QCD. 4.2 The massive Schwinger model with a θterm The Schwinger model is Quantum Electrodynamics in 1+1-dimensions [137]. The euclidean continuum action reads S=Zd2x¯ ψ(x)γµ(∂µ+ieAµ(x)) ψ(x) + m¯ ψ(x)ψ(x) + 1 4F2 µν(x),(4.1) where mis the fermion mass and eis the electric charge or gauge coupling, which has the same dimension as m. After a simple rescaling of the fields the action can be written as S=Zd2x¯ ψ(x)γµ(∂µ+iAµ(x)) ψ(x) + m¯ ψ(x)ψ(x) + 1 4e2F2 µν(x),(4.2) where Fµν(x) = ∂µAν(x)−∂νAµ(x) and γµare 2×2 matrices satisfying the algebra {γµ, γν}= 2gµν.(4.3) where gµν stands for the Euclidean metric tensor. At the classical level this action is invariant in the chiral limit under the UA(1) global transformations ψ→eiαγ5ψ, (4.4) ¯ ψ→¯ ψeiαγ5,(4.5) leading to the conservation of the axial current JA µ(x) = ¯ ψ(x)γµγ5ψ(x).(4.6) 60 CHAPTER 4. SCHWINGER MODEL WITH A θTERM However the axial symmetry is broken at the quantum level because of the axial anomaly, as was discussed in detail in Section 1.3.3. The divergence of the axial current is ∂µJA µ(x) = 1 2πµνFµν(x),(4.7) with µν the antisymmetric tensor, and therefore does not vanish. The axial anomaly induces a topological θterm in the action of the form Stop =iθ 4πZd2xµνFµν(x),(4.8) where the topological charge Q=1 4πRd2xµνFµν(x) is an integer. Our purpose is then to analyze the θdependence of the model described by the action (4.2)+(4.8) S=Zd2x¯ ψγµ(∂µ+iAµ)ψ+m¯ ψψ +1 4e2F2 µν +iθ 4πµνFµν.(4.9) A simple analysis of this model on the lattice suggests that it should undergo a phase transition at some intermediate fermion mass mand θ=π, even at finite lattice spacing. Indeed the lattice model is analytically solvable in the infinite fermion mass limit (pure gauge two-dimensional electrodynamics with topological term) [138, 139], and it is well known that the density of topological charge approaches a non-vanishing vacuum expectation value at θ=πfor any value of the inverse square gauge coupling β, exhibiting spontaneous symmetry breaking. On the other hand by expanding the vacuum energy density in powers of m, treating the fermion mass as a perturbation [140], one gets for the vacuum expectation value of the density of topological charge the following θdependence: h−iqi=mΣsinθ +1 2m2sin (2θ) (χP−χS) + ··· ,(4.10) with Σ the vacuum expectation value of the chiral condensate in the chiral limit and at θ= 0 (Σ = eγee/2π3/2in the continuum limit), and χPand χSthe pseudoscalar and scalar susceptibilities respectively. Equation (4.10) shows how the Z2 symmetry at θ=πis realized order by order in the perturbative expansion of the topological charge in powers of the fermion mass m, and therefore a critical point separating the large and small fermion mass phases is expected. Indeed the model was analyzed in the continuum by Coleman in [22], where he conjectured the existence of a phase transition at θ=π, and some intermediate fermion mass mseparating a ”weak coupling” phase ( e m<< 1), where the Z2 symmetry of the model at θ=πis spontaneously broken, from a ”strong coupling” phase ( e m>> 1) where the Z2symmetry is realized in the vacuum. This 4.3. COMPUTING THE ORDER PARAMETER AS A FUNCTION OF θ61 conjecture was corroborated in [133, 134] using the lattice Hamiltonian approach with staggered fermions, and more recently in [136] using the Grassmann tensor renormalization group and Wilson fermions. 4.3 Computing the order parameter as a function of θ To compute the θdependence of the density of topological charge we use the approach proposed in reference [13]. The only assumption in this approach is the absence of phase transitions at real values of θexcept at most at θ=π. The method is based in extrapolating a suitably defined function to the origin. This function turns out to be very smooth in all the cases considered up to now [21, 109–111], and this makes us confident on the whole procedure. Here we summarize the main steps. From numerical simulations of our physical system at imaginary values of θ= −ih (real values of h), which are free from the severe sign problem, we compute the density of topological charge q(−ih) as a function of h, and introduce the following functions: z= cosh h 2,(4.11) y(z) = q(−ih) tanh h 2 .(4.12) The procedure to find out the density of topological charge at real values of θ relies on scaling transformations [13]. We define the function yλ(z) as yλ(z) = yeλ 2z.(4.13) For negative values of λ, the function yλ(z) allows us to calculate the order parameter tanh h 2y(z)below the threshold z= 1. If y(z) is non-vanishing for any positive z,1then we can plot yλ/y against y. Furthermore, in the case that yλ/y is a smooth function of yclose to the origin, then we can rely on a simple extrapolation to y= 0. Of course, a smooth behavior of yλ/y cannot be taken for granted; however no violations of this rule have been found in the exactly solvable models. 1Even though the possibility of a vanishing y(z) for some value z > 0 cannot be completely excluded, it does not happen for any of the analytically solvable models we know. 68 CHAPTER 4. SCHWINGER MODEL WITH A θTERM 0 0.5 1 1.5 2 0 0.01 0.02 0.03 0.04 0.05 0.06 γ yλ β = 2 β = 3 β = 4 m = 0.0 Figure 4.4: Exponent γfor m= 0.0 and various coupling constants. The shaded areas give an estimation of the ambiguity in the extrapolation to yλ= 0. The continuous red line is the analytic result in the pure gauge theory, corresponding to infinite fermion mass. We plot in Fig. 4.6 the results of each of the independent analysis for an interval of θ. As can be seen, the errors we would obtain by averaging the independent points are fully consistent with the synthetic-data estimation. In Fig. 4.6 we present q(θ) at β= 3 for two masses in the symmetry restored phase, as well as at β= 2 and m= 0.5, in the symmetry broken phase (and also the corresponding analytic results for the pure gauge case at both values of βfor comparison). In Fig. 4.7 we show the results for m= 0 and the three different values of the coupling constant we have simulated. We can clearly see the restoration of the symmetry as we approach θ=π. In Fig. 4.8 we show, for β= 3.0 and m= 0, the order parameter q(θ) in the vicinity of θ=π. Fitting q(θ) near θ=πin the symmetry restored phase allows us to extract the exponent (π−θ), which is related to γby =γ−1.6We present in Table 4.1 our results for . 6The numerical procedure used to extract the two exponents is different, and therefore the 4.5. RESULTS 69 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 0.05 y yλ β = 3.0 , m = 0.0 Figure 4.5: Fit of yversus yλ. Table 4.1: β m  2.0 0.0 0.67(4) 2.0 0.05 0.43(5) 3.0 0.0 0.92(7) 3.0 0.05 0.70(21) 4.0 0.0 0.94(19) To finish this Sec. we want to discuss a little bit more on the results for the massless Schwinger model reported in Fig. 4.7. It is well known that the continuum formulation of the massless Schwinger model shows no θdependence, because the θterm in the action can be canceled by an anomalous chiral transformation which results, although compatible within errors, will also be different. 70 CHAPTER 4. SCHWINGER MODEL WITH A θTERM 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0 0.5 1 1.5 2 2.5 3 q θ β = 3.0 β = 2.0 m = ∞ m = 0.5 m = 0.05 m = 0.0 Figure 4.6: Order parameter as a function of θ. The data at m= 0.0 and m= 0.05 correspond to β= 3.0, whereas the points at m= 0.5 correspond to β= 2. Blue points, corresponding to the results of four independent runs, are also shown, to provide a different estimate of the error. The continuous line labeled m=∞is the pure gauge analytic result for β= 3.0, whereas the dotted line is the corresponding analytic result for β= 2.0. does not change the fermion-gauge action if the fermion mass vanishes. Hence the non-trivial θdependence of the density of topological charge shown in Fig. 4.7 may seem surprising. However, the massless staggered Dirac operator does not have exact zero-modes, and therefore, for a given gauge configuration, a nonzero value of the quantized topological charge Qdoes not imply the existence of a corresponding number of zero-modes in the staggered Dirac operator, as would be the case, for example, with the overlap Dirac operator. What we should expect instead is that, as we approach the continuum limit, the topological charge density vanishes. This is indeed what seems to happen, as is suggested by Fig. 4.9. 4.6. CONCLUSIONS AND OUTLOOK 71 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.018 0 0.5 1 1.5 2 2.5 3 q θ m = 0.0 β = 2.0 β = 4.0 β = 3.0 Figure 4.7: Order parameter as a function of θ, at m= 0.0 and different coupling constants. 4.6 Conclusions and outlook All our results are compatible with the standard lore on this model, and in particular with Coleman’s conjecture on the existence of two distinct phases at θ=π, a symmetry breaking phase at large mass, and a symmetry restored phase at small mass. Our simulations are a proof of concept, and are not extensive enough to determine precisely the position of the critical mass at θ=πor its properties in detail. But the important point is that we have succeeded in calculating the full dependence of the order parameter in θin a gauge theory with fermions and a quantized topological charge, using a method that should, in principle, work also in higher dimensional theories. 72 CHAPTER 4. SCHWINGER MODEL WITH A θTERM 0 0.0005 0.001 0.0015 0.002 0.0025 0.003 2.8 2.85 2.9 2.95 3 3.05 3.1 q θ β = 3.0 , m = 0.0 Figure 4.8: Order parameter as a function of θnear θ=π. 4.6. CONCLUSIONS AND OUTLOOK 73 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0 0.5 1 1.5 2 2.5 3 β2q θ β = 2 β = 3 β = 4 m = 0.0 Figure 4.9: Rescaled topological charge density at m= 0.0 and different coupling constants. 74 CHAPTER 4. SCHWINGER MODEL WITH A θTERM Chapter 5 Preliminary results on two-flavor Schwinger—a pseudofermionic approach Hereafter we present our work with the massive Schwinger model with a θterm and two distinct fermionic species. Although the starting point of the study consists in the application of the same Monte Carlo algorithm that was developed for Chapter 4—a standard implementation of Kogut-Susskind fermions with a Metropolis update—a search for more efficient algorithms proves to be necessary in order to fully test the capabilities of the q(θ) reconstruction approach of [13]. 5.1 Motivation The study of the Schwinger model with a single flavor of massive fermions and the inclusion of a θterm was carried out in the previous chapter with considerable success, since the dependence of the topological charge on θwas determined on the lattice for the whole domain of the vacuum angle, up to θ=π. This was achieved thanks to standard Monte Carlo simulations performed at purely imaginary values of θ, which are the basic input for the reconstruction method of [13]. In these real-action computations, staggered (Kogut-Susskind) fermions were used, and the determinant of the fermionic matrix was calculated at every Metropolis step by extracting numerically all its eigenvalues—a technique as reliable as inefficient. In order to adjust the computation to the one-flavor case, the square root of this determinant needs to be taken; since the actual weight employed by the algorithm depends on the logarithm of the determinant, this rooting procedure amounts to multiplying by a factor of 1/2. If instead we consider a factor of Nf/2, the discussion is valid for the Schwinger model with Nffermion species—in fact, this 75 76 CHAPTER 5. PRELIMINARY RESULTS ON NF= 2 SCHWINGER is the only change required in the algorithm. Having developed an algorithm that can be trivially extended to the multiflavored case, it seems natural to apply it, at least, to the cases that are closely related to one of our main interests during this thesis: the study of topological objects on lattice gauge theories, and its implications in QCD. In fact, this is the case for the two-flavored version of the model: its action holds a U(2) symmetry in the chiral limit, of which its axial U(1) subgroup is broken by the anomaly, much in the same way as QCD. The remaining SU(2) group constitutes a true symmetry of the theory that, contrary to what occurs in low temperature QCD, is exactly preserved—as granted by a Theorem due to Coleman,1a continuous symmetry cannot be spontaneously broken in a two-dimensional system, as long as interactions are kept sufficiently local. This exactly preserved symmetry results to be an interesting property, as long as it is shared by the high temperature phase of QCD. In other words, QCD at high temperatures has in the chiral limit an exactly preserved chiral symmetry (contrary to the less exotic low temperature phase, as was discussed in Section 1.3.3). This fact favors the study the two-flavor Schwinger model as a mechanism to gain insight about the topological properties of finite temperature QCD. Beyond its interest as a toy model of QCD, the Nf= 2 Schwinger model presents a more involved θbehaviour than its single flavored counterpart. The later presents, at θ=π, two distinct phases depending on the coupling e/m. While Psymmetry is spontaneously broken at weak coupling, in the strong coupling (or light mass) limit the vacuum energy density can be expanded in terms of the fermion mass m, its leading contribution being E(θ)∼me cos θ. (5.1) As a consequence, the symmetry is exactly preserved and the topological susceptibility remains finite. By the contrary, the two-flavor version of the model, which has a similar behavior in the weak coupling region, has a more involved θdependence on the vacuum angle. As it was shown by Coleman [22], a strong coupling approximation allows to write the energy density θdependence as E(θ)∼m4 3e2 3cos4 3θ 2,(5.2) which eventually leads to Pexact conservation at θ=π, but with a divergent topological susceptibility—the characteristic of a continuous phase transition. This approximation, valid in principle when e/m >> 1, implies a value of δ= 1/3 for the 1Although this result is proven by Coleman in the context of quantum field theories [147], it is commonly known as Mermin-Wagner theorm, since they arrived to the same conclusions in statistical physics [148]. 5.2. THE MODEL 77 associated critical exponent, which describes how the topological charge density vanishes as θapproaches π. Additionally, the lightest bosons of the spectrum are predicted to be an isotriplet and an isosinglet, the quotient of its masses being √3. It is worth noting that precisely this mass ratio has been recently the subject of some controversy, since a recent work by Azcoiti [149] found a subtantial discrepancy with respect to the original computation of Coleman [22]. Furthermore, Georgi [150] has drown even more attention to this model, by proposing a solution to the three questions posed by Coleman in [22] that could entail the existence of a novel mechanism, capable of generating the appearance of fine-tuning in lowenergy effective theories and, consequently, with promising potential concerning any of the hierarchy problems that afflict the Standard Model. In any case, since a critical point is expected in this model at θ=π, this system poses a relevant challenge for the reconstruction method that was applied during Chapter 4; to this effect, it serves as an additional motivation to this work—a particularly pragmatic one, arguably. 5.2 The model The action of the one-flavor massive Schwinger model with a θterm can be easily generalized to its multi-flavor version by adding an index f, running from 1 to Nf, to the original expression (4.9). In this manner, the action for Nfflavors of equal charge eand mass myields SNf=Zd2x   Nf X 1 ¯ ψf[γµ(∂µ+iAµ) + m]ψf+1 4e2F2 µν +iθ 4πµνFµν   .(5.3) Following the reasoning of the previous chapter, it is possible to discretize the above continuum action by using Kogut-Susskind fermions and the standard Wilson action for the gauge part, as in (4.15). At this point we recall that a procedure commonly known as rooting was needed to get the one-flavor theory from the corresponding lattice action, since staggered fermions are not completely free of the doubling problem—they describe two degenerate species of fermions, in two dimensions. But, as long as we are interested in the two-flavor version, it suffices to consider the exact same action (4.15) and dismiss the rooting step. The next step would imply performing a Monte Carlo simulation, much in the same way as in Chapter 4. However, while in our one-flavor study it was enough to perform a proof-of-concept calculation, our aim with the Nf= 2 case is to determine more involved quantities, such as the critical exponent of the expected θ=πphase transition. Even if the brute force approach of the previous chapter was able to deliver results in, roughly speaking, a few months of computer 84 CHAPTER 5. PRELIMINARY RESULTS ON NF= 2 SCHWINGER 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 hqi h Method [29] Pseudofermions, n= 10 Figure 5.5: Results for the two-flavor Schwinger model in a 16 ×16 lattice, with β= 3 and mf= 0.12. δmax = 0.25 and n= 10. The pseudofermion approach is confronted with the method of Chapter 4, with great agreement for small values of h, although systematic deviations are otherwise observed. Chapter 6 Exploratory ghost-gluon study on αswith HISQ fermions In this last chapter we put aside the study of topological effects to face a nonperturbative quantity of the utmost relevance in QCD—its running coupling constant αs. Our intention is to explore the dependence of the coupling on the momentum transfer q, by means of a purely gluonic method, following Sternbeck et al [26]. In this way, it is possible to compute gluon and ghost propagators, G(q2) and D(q2), in the lattice, eventually leading to the determination of αs(q2). In the following sections we motivate the topic and review briefly the most relevant mathematical relations, describing with some detail the technical issues involved in the propagators determination. After that, our results—obtained from large sets of configurations generated by the MILC collaboration—are presented. The chapter is finished with some considerations concerning how the present study could be extended. 6.1 Motivation The determination of the αscoupling constant constitutes a very active field of research. In fact, the Particle Data Group periodically provide a global average of this quantity, including both theoretical calculations of lattice simulations and experimental determinations, such as from hadronic τdecays or e+e−annihilation processes, that allow to give an estimate—at a given scale, usually that of the Z boson—of the strong coupling αs[2]. Remarkably enough, the last decade has been marked by a setback of the precision achieved when calculating the global average of this quantity, in part due to the existence of previously underestimated sources of systematic error that were present in lattice computations. In this situation, there cohabit several independent approaches within the lat85 86 CHAPTER 6. EXPLORATORY GHOST-GLUON STUDY ON αS tice formalism, developed by a number of international collaborations such as HPQCD [4, 151, 152], PACS-CS [153], ETM [154], and other research groups [155, 156]. These approaches differ in a number of technicalities, including how fermions are implemented on the lattice. Moreover, different observables can be studied in order to extract αs, an example being heavy quark correlators in [152], or the determination of the static potential in [156], which in fact gives a value for αs(MZ) that exhibits some tension with the rest of lattice predictions. Our intention in the current chapter is to explore the potential of combining the ghost-gluon vertex technique of Sternbeck et al [26], which was already applied to twisted mass fermions in [154], with the highly improved staggered action (HISQ), used by the HPQCD collaboration [78]. To this aim, we will analyze ghost and gluon fields in a collection of ensembles generated by the MILC collaboration with the HISQ action [157]. 6.2 αsand the ghost-gluon vertex Before going any further, it should be noted that the running coupling αsis not a physically observable quantity. Instead, it acquires a precise meaning only in the context of perturbation theory. Then, in order to compare two given results for the coupling it is necessary to take into account the renormalization scheme. Typical choices in the literature are momentum subtraction schemes, the more common including MS, MS and MOM. While the standard computed value of αs(MZ) is tipically given in the literature in the MS scheme [2], the present approach is defined in a MOM scheme, in which renormalization constants are defined by requiring twoand three-point functions to equal their tree level expressions at a given energy scale µ[158]. There exist a number of ways that allow to calculate the value of the running coupling αson the lattice. As it has been just mentioned, the present work follows the approach of Sternbeck et al [26], although with the focus set in a different region (since we are not particularly interested in the infrared limit of αs). Hereafter, we review the essentials of the method. The computation of αsin the lattice starts by realizing how ghost and gluon propagators can be exploited. Its dressing functions can be used to determine the running coupling as a renormalization group invariant, in a momentum subtraction scheme [25], as αs(q2) = g2 0 4πZD(q2)Z2 G(q2),(6.1) where ZDand ZGare respectively the bare dressing functions of gluon and ghost propagators; we discuss how to compute these functions in the next sections. 6.3. THE GLUON PROPAGATOR 87 Although in usual lattice computations, as we reviewed in Chapter 1, fixing the gauge is not necessary, the mathematical expressions for ghost and gluon propagators adopt a simpler form in the Landau gauge—which in fact makes the whole computation feasible. Consequently, in what follows all expressions will be understood to be valid within the Landau gauge. 6.3 The gluon propagator To begin with, the starting point is the standard 4-dimensional lattice of Vsites, and 4Vlink variables Ux,µ ∈SU(Nc= 3). The lattice gluon fields, which live in the mid-point of each link, Ax,µ ≡Aµ(x+ ˆµ/2), are defined as Ax,µ := 1 2i(Ux,µ −U† x,µ)−1 6iTr(Ux,µ −U† x,µ).(6.2) Additionally, we recall that the color components Aa x,µ of the gluon field can be computed as Aa x,µ := 2 Tr(TaAx,µ)=2·Im Tr(TaUx,µ),(6.3) where the definition (6.2) has been used. With these expressions, we can compute the bare gluon propagator on the lattice as Dab µν(k) = D˜ Aa µ(k)˜ Ab ν(−k)EU,(6.4) where ˜ Aµ=˜ Aa µTaare the Fourier transformed gluon fields. In other words, Dab µν(q(k)) = 1 V*X x,y Aa x,µAb y,νeik·(x+ˆµ/2)e−ik·(y+ˆν/2)+U ,(6.5) where the momentum q(k) is given by qµ(kµ) = 2 asin πkµ Lµ.(6.6) If we assume now that Dab µν(q(k)) has the same tensor structure than its continuum counterpart, Dab µν(q) = δab δµν −qµqν q2D(q2),(6.7) it suffices a bit of algebra to obtain an expression for the scalar part of the propagator D(q2), D(q2) = 1 (D−1)(N2 c−1) X aµ Daa µµ(q),(6.8) 88 CHAPTER 6. EXPLORATORY GHOST-GLUON STUDY ON αS which is related with the dressed propagator simply by ZD(q2)≡q2D(q2).(6.9) In practical terms, the key step in this computation—in terms of computational complexity—is the Fourier transformation of the gauge fields. Fortunately, Fast Fourier Transform algorithms allow to compute very efficiently expression (6.5), especially taking into account that a single application of the algorithm delivers the propagator for every lattice value of the momenta qµat once. 6.4 The ghost propagator Following [27], the ghost propagator in the lattice is defined in the Landau gauge as Gab(k) = a2*X xy M−1ab xy eik·(x−y)+=δabG(q),(6.10) where the real symmetric matrix Mis the Fadeev-Popov operator, defined by Mab xy =X µAab x,µδx,y −Bab x,µδx+ˆµ,y −Cab x,µδx−ˆµ,y(6.11) with Aab x,µ =Re Tr {Ta, Tb}(Ux,µ +Ux−ˆµ,µ),(6.12) Bab x,µ = 2 ·Re Tr TbTaUx,µ,(6.13) Cab x,µ = 2 ·Re Tr TaTbUx−ˆµ,µ.(6.14) In order to compute (6.10), the following system of equations needs to be solved    Max,bycby c=δac cos (k·x), Max,bysby c=δac sin (k·x).(6.15) The 8V-component vectors cc,scare computed with the conjugate gradient method and can be used to determine the inverse of M. Together with (6.10), and assuming the tensor structure of the continuum, Gab(qµ) = δabG(q2), we have G(q2) = 1 (N2 c−1) X ax [cos (k·x)cax a+ sin (k·x)sax a],(6.16) with the corresponding dressed propagator being given by ZG(q2)≡q2G(q2).(6.17) 6.5. RESULTS 89 β m0 l/m0 sm0 sm0 cN3 s×Nta(fm) # of cnfgs 6.00 1/5 0.0509 0.0635 243×64 0.1218(7) 1053 6.30 1/5 0.0370 0.0440 323×96 0.0879(5) 1008 6.72 1/5 0.0240 0.0286 483×144 0.0573(4) 1017 Table 6.1: Details of the three ensembles studied. In contrast with the gluon determination of the previous section, the steps here described are much more expensive in computational terms. In particular, numerically solving the system of equations (6.15) is a very demanding task which, at large lattice sizes—as the ones studied in this chapter are—requires large resources, both in terms of memory and processing time. Furthermore, lattice artifacts are expected to be more intense both at large and at off-diagonal momenta, due to the lack of rotational symmetry on the lattice in the latter case [159]. For this reasons, the ghost propagator, and as a consequence also αs, have been computed only for a handful of selected diagonal momenta. 6.5 Results We have analyzed three sets of configurations, made available by the MILC collaboration [157]. The details of their parameters are summarized in Table 6.1. Prior to the propagators determination, we fixed every configuration to Landau gauge. To this end, it is necessary to make use of an iterative optimization algorithm. As is well known, this type of algorithms could suffer from a severe critical slowing down problem. In the case of large lattices, as some of the ones analyzed here, this obstacle can turn insurmountable. However, we have evaded this difficulty by applying a Fourier-accelerated algorithm, originally proposed by Davies et al [160], which allows to alleviate the computational overhead, making the gauge fixing procedure feasible. In this process, the fixing procedure was stopped only when every local gauge field verified the transversality condition—the lattice version of ∂µAµ= 0—up to Θ <10−14, with the same definition of [160]. Such level of precision was proved to be necessary, since the propagators are quite sensitive to the gauge condition. Once the whole sets were fixed to Landau gauge, we have computed both propagators, G(q2) and D(q2), for a total of seven diagonal momenta, k= (n, n, n, n) for n= 1,...,7.(6.18) As we mentioned earlier, the ghost computation, and in a lesser way the gauge 90 CHAPTER 6. EXPLORATORY GHOST-GLUON STUDY ON αS # of configurations 1 1053 Landau gauge-fixing time (core-h) 7.1 7.1k Propagators for 7 momenta (core-h) 19.7 20.8k Total time (7 momenta) (core-h) 26.928.3k Each additional momentum (core-h) +3.0k Table 6.2: Distribution of computing times involved in the 243×64 ensemble. fixing procedure, are very demanding in terms of computing resources. This being the case, actual calculations have required to be performed in large cluster facilities, which provide both the computing power and the memory needed to store the largest configurations. To this end, the University of Cambridge computing services have been used, including the cluster Darwin, and its 2017 update CSD3. For the first set, of volume 243×64, 28k core-hours were used in Darwin; to exemplify how these are distributed, see Table 6.2. For the 323×96 set, a total of 135k core-hours, also in Darwin, were used. Finally, the ensemble corresponding to the bigger lattice size, 483×144, was analyzed with a total cost of 452k core-hours of the newer cluster CSD3. Our results for both bare unrenormalized propagators are shown in Figure 6.1, as a function of the squared lattice momenta in physical units. By means of Eq. 6.1, the previous results can be applied to determine αs. Taking into account that, according with the particular implementation of the HISQ action in the analysed MILC configurations, g2 0≡5 3 2Nc β=10 β,(6.19) our final results for the unrenormalized version of the running coupling are shown in Fig. 6.2. As a preliminary conclusion, these results can be checked to be qualitatively in agreement with those of the literature, see e.g. [25,26,154]. In any case, in order to provide an estimate for the αs(MZ) in the MS renormalization scheme, the corresponding βfunction should be carefully studied and integrated, together with a thorough study of present volume effects and lattice artifacts, since the knowledge of its dependence would allow to estimate the systematic errors of the method. Finally, the relatively low value of the estimated statistical errors in the studied ensembles drive us to conclude that completing the present analysis would be worthwhile, potentially leading to a valuable contribution to the world average of the strong running coupling. 6.5. RESULTS 91 0 1 2 3 4 1 10 0 1 2 3 4 5 6 1 10 q2[GeV 2] G(q)−243×64 323×96 483×144 q2[GeV 2] D(q)−243×64 323×96 483×144 Figure 6.1: Ghost and gluon bare propagators, in physical units. 0 0.5 1 1.5 2 1 10 q2[GeV2] 243×64 323×96 483×144 Figure 6.2: The running coupling αs(q2) as defined in Eq. 6.1. Shaded areas connect 1σintervals. 92 CHAPTER 6. EXPLORATORY GHOST-GLUON STUDY ON αS Conclusions and outlook In the beginning of this thesis we emphasized the key role of lattice QCD as the main tool capable of dealing with non-perturbative phenomena, in QCD and beyond. Their successes are many, as it has been documented in Chapter 1. Up to the present day, there exist several international collaborations that work intensively to provide theoretical estimates that match the ever-increasing precision of the experimental measures. But, notwithstanding the achievements in phenomenology grounds, there exist open questions that seem to evade any attempt of resolution. Indeed, this is the case of QCD when the topological θterm is included into its action, as it has been reviewed in Chapter 2. Deeply connected with the strong CP problem, and of great importance in axion physics, the study of the θdependence of QCD has been the long-sighted goal of this thesis. Certainly, this is an ambitious objective, since there exist many difficulties when considering θQCD on the lattice, one of the most fundamental being the presence of a severe sign problem that makes the system unattainable to standard Monte Carlo simulations. In this scenario, our efforts have been first directed towards the testing of a promising method that would in principle allow to reconstruct the full θdependence of observables as the topological charge [13]. Thereby, in Chapter 3, following [28], we have crosschecked a previous work that, making use of the reconstruction method of Azcoiti et al, found signs of a rich phase structure in the antiferromagnetic Ising model within an imaginary magnetic field—which can be thought of as a θterm [21]. By making use of a combination of analytical and numerical—but exact—techniques, our results have supported the previous qualitative picture, thus confirming the potential of the reconstruction method in a theory that, when compared to θQCD, holds a more intricate phase diagram. In the same spirit, Chapter 4 has been devoted to the study of the massive one-flavor Schwinger model with a θterm in the action, which frequently serves as a toy model of QCD, as presented in [29]. To apply the reconstruction method [13] in this system implies a stringent exam concerning its possibilities of being successfully applied to full QCD, as it shares many of its potential obstacles— but not all, since it holds a very elementary definition of the topological charge, in contrast with the difficulties that appear in the SU(3) 4d theory. Actually, we have 93 100 CHAPTER 6. EXPLORATORY GHOST-GLUON STUDY ON αS Appendix A Computation of the cumulants κn In order to use expressions (3.17) and (3.20), we need to compute the cumulants κn. The nth cumulant can be calculated in terms of the first nnoncentral moments µ0 n, µ0 n≡* X <ij> sisj!n+m1,m2 ,(A.1) by means of the recursion relation (3.16). The summation over < ij > runs over each couple of neighboring spins, or in other words, over each link. Two neighboring spins always belong to different sublattices. Before going further, let us comment on two intermediate results. First, we consider a lattice of Nspins, the magnetization of which is the sum m=Pisi, and ask about the expected value of the product of nof these spins at fixed m(or fixed N+, the number of positive spins), that is, hs1s2···snim. One can perform this calculation by means of the microcanonical formalism, arriving at hs1s2···snim=1 N N+ n X k=0 (−1)kn k N−n N+−n+k.(A.2) In the above expression, kcan be read as the number of negative spins in the product s1s2···sn. In this way, the first summand, k= 0, counts the number of states with zero negative spins in the product s1s2···snand multiplies it by the expectation value of the product in this case, (−1)0= 1. The second one, k= 1, does the same for one negative spin in s1···sn, and so on. Dividing the sum by the total number of configurations with magnetization m= 2N+/N −1, one obtains the previous expectation value at fixed m. Secondly, consider an observable O(m1, m2) in our two sublattice system, with a dependence on m1and m2such as we can write it as O1(m1)O2(m2). In this case, from the definition 101 102 APPENDIX A. COMPUTATION OF THE CUMULANTS κN (3.12) of the expectation value at fixed m1and m2, we have hO1(m1)O2(m2)im1,m2=hO1(m1)im1hO2(m2)im2.(A.3) This immediately applies to the spin product s1s2···sn. We can always divide it into two products sa···sband sα···sβ, each one containing the spins of one of the sublattices, and then hs1s2···snim1,m2=hsa···sbim1hsα···sβim2.(A.4) With the previous couple of results, we come back to Eq. (A.1), and apply the linearity of the expectation value, arriving at µ0 n=X <ij>,<kl>,··· ,<pq> hsisjsksl···spsqim1,m2,(A.5) which is the sum of the expectation values of the product of nlinks, running over all permutations with repetitions of these links. Then, in every summand we have the product of 2nspins, in some cases with some of them identical. Taking into account that s2 i= 1 ∀i, each summand can be reduced to the expectation value of the product of n1+n2different spins, n1and n2being the number of spins in each sublattice. Since by means of Eq. (A.2) we already have an expression that computes hs1···snim, the problem is reduced to count how many summands in Eq. (A.5) have (n1, n2) spins. We call these numbers geometrical factors, and denote them by G(n1, n2). Following this convention, we can write the nth central moment as µ0 n=X {n1,n2}G(n1, n2)hsa···sb | {z } n1spins im1hsα···sβ | {z } n2spins im2,(A.6) where the sum runs over the couples of integers (n1, n2) the sum of which is even and less than or equal to n. The computation of the geometrical factors G(n1, n2) can be done by hand for the first few cumulants. As an example, for the second noncentral moment µ0 2we have to compute four cases: the two links being the same (sharing both spins), sharing only one spin belonging to the first or the second sublattice, and finally not sharing any spin at all. That is, in terms of the previous notation, {(n1, n2)}={(0,0),(2,0),(0,2),(2,2)}.(A.7) The factors G(n1, n2) can be computed easily in this case, even for an hypercubic lattice of arbitrary dimension d, arriving at the following expression for the second moment µ0 2=Ndh1i+Nd(d−1)(hs1s2im1+hs1s2im2) +Nd(Nd −2(d−1) −1)hs1s2im1hs1s2im2. (A.8) A.1. TRANSLATIONAL SYMMETRY 103 We can use this expression to calculate the second cumulant κ2, κ2=µ0 2−µ02 1 N→∞ −−−→ Nd(m2 1−1)(m2 2−1),(A.9) where we have taken the thermodynamic limit, keeping only the terms of order O(N), which is the leading order for all cumulants. Subleading orders can be preserved if needed, but they are not relevant for our paper. The difficulty of the previous computation escalates quickly with the order nof the cumulant, and it is quite cumbersome for just n≥4. In order to get beyond this limitation, we have developed a program which computes the geometrical factors G(n1, n2) numerically for a finite L×Lbidimensional lattice. Since these factors G(n1, n2) are polynomials in Nof order ≤n(and with integer coefficients), we can run the program for lattices of n+ 1 different sizes, obtaining a set of (N,G(N)) points, which we can use to recover the exact integer coefficients of each geometrical factor, by means of the Lagrange interpolation formula. The basic idea of the program is very simple. We just construct a periodic rectangular L×Mlattice, with L, M > n,nbeing the order of the cumulant we want to compute. With this restriction we avoid products of links crossing the entire lattice, that would not appear in the thermodynamic limit for any finite cumulant. Once we have this, we start a loop running over all the permutations with repetitions of nlinks, and perform the following steps, •We have a product of nlinks, or equivalently 2nspins, s1···s2n. •Recursively, we remove couples of equal spins from this product. •We classify the remaining product by the number of spins in each sublattice, (n1, n2). •We add one to the geometric factor G(n1, n2) and proceed to the next iteration. When the algorithm finishes, we obtain all the G(n1, n2) values for a given N= LM. The computational cost is associated to the number of iterations of the main loop, which grows as (LM)n, that is, exponentially with the order of the cumulant. In practice, we have only reached the computation of the fourth cumulant with this program. However, a number of optimizations can be implemented in order to reach higher order cumulants, which we summarize in what follows. A.1 Translational symmetry Our lattice is symmetric under translations, implying that all geometrical factors are proportional to Nd, the number of links. Fixing, e.g., the first link of the 104 APPENDIX A. COMPUTATION OF THE CUMULANTS κN product, one obtains the same G(n1, n2), but divided by a common factor Nd. The same factor is gained in the overall speed of the program. In addition to that, the degree of the polynomials G(n1, n2) is also reduced by one, and it suffices with n(instead of n+ 1) different sizes in order to recover the Ndependence. One can go even further by realizing that the geometrical factor corresponding to non-neighboring links, G(n, n), is the only one with maximum degree Nn−1. This allows us to express it in terms of the remaining factors, 1 NdG(n, n) = (Nd)n−1 −1 Nd X {(n1,n2)}\(n,n)G(n1, n2),(A.10) which are only of order n−2 or less. This means that it is enough to run the program for n−1 lattice sizes, compute all the geometrical factors but G(n, n) via the Lagrange interpolator, and then with the previous expression find the N dependence of this last factor. A.2 From permutations to combinations The product of links commutes, so its contribution to the geometrical factors is the same regardless of the order. Then, we can change the main loop over permutations with repetition to a loop over combinations with repetition, by taking into account the multiplicity of each combination. Schematically, we perform Pi,j,...,k contrib(lilj···lk) →X i≤j≤···≤k mult ×contrib(lilj···lk),(A.11) where contrib represents a function in our program that takes a product of links and returns the contribution to the geometrical factors. If there are rdifferent links, each one appearing k1, . . . , krtimes, the multiplicity of the combination is given by mult = n! k1!···kr!.(A.12) A.3 Blocks - Grouping links together Many of the link products have few, if any, repeated spins, and their contributions to the geometrical factors can be counted without having to analyze one by one each of them. This is possible by grouping them in sets of links that we will call in A.4. CLUSTERS OF BLOCKS 105 what follows blocks, and replacing the loop over link products by a loop over block products. When the blocks in a product are not neighbors (i.e., they do not have any common spin), we do not need to perform the computation link by link and the contribution can be summed up trivially. Let b1and b3be two non-neighboring blocks, each one composed by Nblinks, and let us denote the contributions to the geometrical factors by λ(n1, n2), where λis an integer counting how many products of links have n1(n2) spins in the first (second) sublattice. Then we have contrib(b1b3) = N2 b(2,2),(A.13) or in general, for the product of knon-neighboring blocks, Nk b(k, k). Following this strategy, we divide our lattice into unidimensional blocks of 2Mlinks, in a way that the jth block, bj, contains all links the first spin of which belongs to the jth column. As a consequence, bjis a neighbor of blocks j−1 and j+ 1, and, taking into account the boundary conditions, b0and bL−1are neighbors too. When we have a product of neighboring blocks, we proceed as before, analyzing the link products one by one, and there is no computational saving. But when the nblocks are not neighbors, we move from (Nd)niterations to a single one. A.4 Clusters of blocks The block method, as defined above, fails to save any computation time if two or more blocks are neighbors in a given block product. However, we can extend the method by dividing each block product into several subproducts, which we will denote as clusters. In each cluster, one can always connect one block to another by the equivalence relation of being neighbors (sharing spins). And in the same way, in each product different clusters never share any spin. This allows us to compute the contributions of each cluster separately, and then compose them with the following law, λ1(a, b)⊕λ2(c, d) = λ1λ2(a+c, b +d).(A.14) If the contributions of the clusters involve more than one geometrical factor, linearity applies, X ab λab(a, b)⊕X cd λcd(c, d) = X ab,cd λabλcd(a+c, b +d).(A.15) Processing one cluster with kblocks takes a computing time proportional to (Nd)k. So dividing the whole block product in smaller clusters implies for almost every 106 APPENDIX A. COMPUTATION OF THE CUMULANTS κN block product a significant amount of time saved. Only when all the blocks are part of the same cluster there is no speed up. Another major optimization can be performed by realizing that translational invariance can also be applied here, since a given cluster, say b0b1b1, and any of its translations, b0+tb1+tb1+t, have the same contribution to the geometrical factors. Then, when a cluster is going to be computed, we can express it in terms of its equivalence class, compute its contribution, and store it in memory. Every time one of its translations appears, we just take the value from the memory, saving a lot of computing time. In addition to that, once we have computed the factors G(n1, n2) for the first size L×M, we know in advance all the cluster contributions for any L0×Mlattice (the blocks keep its size constant). Since almost all the computing time is spent in figuring out the cluster contributions, we reduce in this way the full problem of computing the geometrical factors in lattices of n−1 different sizes to only one size, the smallest one, M×M. In practice, the time spent by the rest of the sizes needed is barely the 1 −2% of that of the first size. A.5 Computation of a cluster The last optimization concerns the computation of the clusters themselves. Until now it is done simply by performing a loop over each possible permutation of links belonging to each of the blocks in the cluster. However, one can go one step further and divide the blocks composing the cluster into smaller sets, that we will call sites. A site is simply the set of two links the first spin of which lies in the site i, j, that is, site(i, j)≡ {sijsi+1,j, sijsi,j+1}.(A.16) With this new subdivision, we can apply in the same way the techniques described above. In order to compute the cluster b1. . . bk, we start a loop over every permutation of sites s1. . . sk, with si∈bi. Each site product is divided into clusters, the contributions of which can be summed with Eq. (A.15) and are calculated by performing another loop over each link product (2kiterations for a site product of kelements). Finally, by summing up each site product contribution, we obtain the whole cluster contribution. All the described optimizations do not remove the exponential dependence on nof the algorithm. However, they allow us to reach the eighth cumulant, which takes about three days of computing time in a modern laptop. Bibliography [1] J. C. Collins, D. E. Soper, and G. F. Sterman, “Factorization of Hard Processes in QCD,” Adv. Ser. Direct. High Energy Phys. 5(1989) 1–91, arXiv:hep-ph/0409313 [hep-ph]. [2] Particle Data Group, M. Tanabashi et al., “Review of Particle Physics,” Phys. Rev. D 98 (Aug, 2018) 030001. [3] K. G. Wilson, “Confinement of quarks,” Phys. Rev. D 10 (Oct, 1974) 2445–2459. [4] HPQCD and UKQCD Collaborations and MILC Collaboration and HPQCD and Fermilab Lattice Collaborations, C. T. H. Davies, E. Follana, A. Gray, G. P. Lepage, Q. Mason, M. Nobes, J. Shigemitsu, H. D. Trottier, M. Wingate, C. Aubin, C. Bernard, T. Burch, C. DeTar, S. Gottlieb, E. B. Gregory, U. M. Heller, J. E. Hetrick, J. Osborn, R. Sugar, D. Toussaint, M. D. Pierro, A. El-Khadra, A. S. Kronfeld, P. B. Mackenzie, D. Menscher, and J. Simone, “High-precision lattice qcd confronts experiment,” Phys. Rev. Lett. 92 (Jan, 2004) 022001. [5] S. D¨urr, Z. Fodor, J. Frison, C. Hoelbling, R. Hoffmann, S. D. Katz, S. Krieg, T. Kurth, L. Lellouch, T. Lippert, K. K. Szabo, and G. Vulvert, “Ab initio determination of light hadron masses,” Science 322 no. 5905, (2008) 1224–1227, http://science.sciencemag.org/content/322/5905/1224.full.pdf. [6] S. Borsanyi, S. Durr, Z. Fodor, C. Hoelbling, S. D. Katz, S. Krieg, L. Lellouch, T. Lippert, A. Portelli, K. K. Szabo, and B. C. Toth, “Ab initio calculation of the neutron-proton mass difference,” Science 347 no. 6229, (2015) 1452–1455, http://science.sciencemag.org/content/347/6229/1452.full.pdf. [7] K. G. Wilson, “Ab initio quantum chemistry: A source of ideas for lattice gauge theorists,” Nuclear Physics B - Proceedings Supplements 17 (1990) 82 – 92. 107 108 BIBLIOGRAPHY [8] L. Fortnow, “The status of the p versus np problem,” Commun. ACM 52 no. 9, (Sept., 2009) 78–86. [9] G. Parisi, “On complex probabilities,” Physics Letters B 131 no. 4, (1983) 393 – 395. [10] J. R. Klauder, “STOCHASTIC QUANTIZATION,” Acta Phys. Austriaca Suppl. 25 (1983) 251–281. [11] J. Berges and I.-O. Stamatescu, “Simulating nonequilibrium quantum fields with stochastic quantization techniques,” Phys. Rev. Lett. 95 (Nov, 2005) 202003. [12] V. Azcoiti, G. Di Carlo, A. Galante, and V. Laliena, “New proposal for numerical simulations of thetavacuum - like systems,” Phys. Rev. Lett. 89 (2002) 141601, arXiv:hep-lat/0203017 [hep-lat]. [13] V. Azcoiti, G. Di Carlo, A. Galante, and V. Laliena, “theta vacuum systems via real action simulations,” Phys. Lett. B563 (2003) 117, arXiv:hep-lat/0305005 [hep-lat]. [14] E. Witten, “Analytic Continuation Of Chern-Simons Theory,” AMS/IP Stud. Adv. Math. 50 (2011) 347–446, arXiv:1001.2933 [hep-th]. (Aug, 2010). [15] E. Witten, “A New Look At The Path Integral Of Quantum Mechanics,” arXiv:1009.6032 [hep-th]. (Sep, 2010). [16] AuroraScience Collaboration, M. Cristoforetti, F. Di Renzo, and L. Scorzato, “New approach to the sign problem in quantum field theories: High density qcd on a lefschetz thimble,” Phys. Rev. D 86 (Oct, 2012) 074506. [17] K. Langfeld, B. Lucini, and A. Rago, “Density of states in gauge theories,” Phys. Rev. Lett. 109 (Sep, 2012) 111601. [18] K. Langfeld, B. Lucini, R. Pellegrini, and A. Rago, “An efficient algorithm for numerical computations of continuous densities of states,” Eur. Phys. J. C76 no. 6, (2016) 306, arXiv:1509.08391 [hep-lat]. [19] K. Langfeld, “Density-of-states,” PoS LATTICE2016 (2017) 010, arXiv:1610.09856 [hep-lat]. BIBLIOGRAPHY 109 [20] V. Matveev and R. Shrock, “On properties of the ising model for complex energy/temperature and magnetic field,” Journal of Physics A: Mathematical and Theoretical 41 no. 13, (Mar, 2008) 135002. [21] V. Azcoiti, E. Follana, and A. Vaquero, “Progress in numerical simulations of systems with a θ−vacuum like term: The two and three-dimensional Ising model within an imaginary magnetic field,” Nucl. Phys. B851 (2011) 420, arXiv:1105.1020 [hep-lat]. [22] S. R. Coleman, “More About the Massive Schwinger Model,” Annals Phys. 101 (1976) 239. [23] F. Fucito, E. Marinari, G. Parisi, and C. Rebbi, “A proposal for monte carlo simulations of fermionic systems,” Nuclear Physics B 180 no. 3, (1981) 369 – 377. [24] A. Deur, S. J. Brodsky, and G. F. de T´eramond, “The qcd running coupling,” Progress in Particle and Nuclear Physics 90 (2016) 1 – 74. [25] A. Sternbeck, E.-M. Ilgenfritz, M. M¨uller-Preussker, and A. Schiller, “Towards the infrared limit in su(3) landau gauge lattice gluodynamics,” Phys. Rev. D 72 (Jul, 2005) 014507. [26] A. Sternbeck, K. Maltman, L. von Smekal, A. G. Williams, E. M. Ilgenfritz, and M. Muller-Preussker, “Running alpha(s) from Landau-gauge gluon and ghost correlations,” PoS LATTICE2007 (2007) 256, arXiv:0710.2965 [hep-lat]. [27] R. Aouane, V. G. Bornyakov, E.-M. Ilgenfritz, V. K. Mitrjushkin, M. M¨uller-Preussker, and A. Sternbeck, “Landau gauge gluon and ghost propagators at finite temperature from quenched lattice qcd,” Phys. Rev. D 85 (Feb, 2012) 034501. [28] V. Azcoiti, G. Di Carlo, E. Follana, and E. Royo-Amondarain, “Antiferromagnetic ising model in an imaginary magnetic field,” Phys. Rev. E96 (Sep, 2017) 032114. [29] V. Azcoiti, E. Follana, E. Royo-Amondarain, G. Di Carlo, and A. Vaquero Avil´es-Casco, “Massive schwinger model at finite θ,” Phys. Rev. D 97 (Jan, 2018) 014507. [30] C. M. G. Lattes, H. Muirhead, G. P. S. Occhialini, and C. F. Powell, “PROCESSES INVOLVING CHARGED MESONS,” Nature 159 (1947) 694–697. [,42(1947)]. 116 BIBLIOGRAPHY [107] G. Aarts, F. A. James, E. Seiler, and I.-O. Stamatescu, “Adaptive stepsize and instabilities in complex Langevin dynamics,” Phys. Lett. B 687 (2010) 154–159, arXiv:0912.0617 [hep-lat]. [108] E. Seiler, “Status of Complex Langevin,” EPJ Web Conf. 175 (2018) 01019, arXiv:1708.08254 [hep-lat]. [109] V. Azcoiti, G. Di Carlo, A. Galante, and V. Laliena, “θdependence of the CP9model,” Phys. Rev. D69 (2004) 056006, arXiv:hep-lat/0305022 [hep-lat]. [110] V. Azcoiti, G. Di Carlo, and A. Galante, “Critical Behaviour of CP1at θ=π, Haldane’s Conjecture, and the Relevant Universality Class,” Phys. Rev. Lett. 98 (2007) 257203, arXiv:0710.1507 [hep-lat]. [111] V. Azcoiti, G. Di Carlo, E. Follana, and M. Giordano, “Critical behaviour of the O(3) nonlinear sigma model with topological term at θ=πfrom numerical simulations,” Phys. Rev. D86 (2012) 096009, arXiv:1207.4905 [hep-lat]. [112] C. Gattringer and K. Langfeld, “Approaches to the sign problem in lattice field theory,” Int. J. Mod. Phys. A31 no. 22, (2016) 1643007, arXiv:1603.09517 [hep-lat]. [113] S.-Y. Kim, “Yang-lee zeros of the antiferromagnetic ising model,” Phys. Rev. Lett. 93 (Sep, 2004) 130604. [114] L. Onsager, “Crystal statistics. i. a two-dimensional model with an order-disorder transition,” Phys. Rev. 65 (Feb, 1944) 117–149. [115] T. D. Lee and C. N. Yang, “Statistical theory of equations of state and phase transitions. ii. lattice gas and ising model,” Phys. Rev. 87 (Aug, 1952) 410–419. [116] B. M. McCoy and T. T. Wu, “Theory of Toeplitz Determinants and the Spin Correlations of the Two-Dimensional Ising Model. II,” Phys. Rev. 155 (1967) 438. [117] E. Ising, “Beitrag zur theorie des ferromagnetismus,” Zeitschrift f¨ur Physik 31 no. 1, (Feb, 1925) 253–258. [118] V. Matveev and R. Shrock, “Complex temperature properties of the 2-D Ising model with Beta H = (+- i pi/2),” J. Phys. A28 (1995) 4859–4882, arXiv:hep-lat/9412105 [hep-lat]. BIBLIOGRAPHY 117 [119] B. M. McCoy and T. T. Wu, The Two-Dimensional Ising Model. Harvard University Press, Cambridge, MA, 1973. [120] P. de Forcrand and T. Rindlisbacher, “The density of states method applied to the ising model with an imaginary magnetic field,” 8, 2016. https:// conference.ippp.dur.ac.uk/event/530/session/3/contribution/58. [121] V. Azcoiti and A. Galante, “Parity and ct realization in qcd,” Phys. Rev. Lett. 83 (Aug, 1999) 1518–1520. [122] R. D. Peccei, “Why PQ?,” AIP Conf. Proc. 1274 (2010) 7, arXiv:1005.0643 [hep-ph]. [123] C. Bonati, M. D’Elia, H. Panagopoulos, and E. Vicari, “Change of θ dependence in 4d SU(n) gauge theories across the deconfinement transition,” Phys. Rev. Lett. 110 (Jun, 2013) 252003. [124] C. Bonati, M. D’Elia, M. Mariti, G. Martinelli, M. Mesiti, F. Negro, F. Sanfilippo, and G. Villadoro, “Axion phenomenology and θ-dependence from Nf= 2 + 1 lattice QCD,” JHEP 03 (2016) 155, arXiv:1512.06746 [hep-lat]. [125] P. Petreczky, H.-P. Schadler, and S. Sharma, “The topological susceptibility in finite temperature QCD and axion cosmology,” Phys. Lett. B762 (2016) 498, arXiv:1606.03145 [hep-lat]. [126] S. Borsanyi et al., “Calculation of the axion mass based on high-temperature lattice quantum chromodynamics,” Nature 539 no. 7627, (2016) 69, arXiv:1606.07494 [hep-lat]. [127] V. Azcoiti, “Topology in the SU(Nf) chiral symmetry restored phase of unquenched QCD and axion cosmology,” Phys. Rev. D94 no. 9, (2016) 094505, arXiv:1609.01230 [hep-lat]. [128] W. Bietenholz, K. Cichy, P. de Forcrand, A. Dromard, and U. Gerber, “The Slab Method to Measure the Topological Susceptibility,” PoS LATTICE2016 (2016) 321, arXiv:1610.00685 [hep-lat]. [129] V. Azcoiti, “Topology in the SU(Nf) chiral symmetry restored phase of unquenched QCD and axion cosmology. II.,” Phys. Rev. D96 no. 1, (2017) 014505, arXiv:1704.04906 [hep-lat]. [130] V. Azcoiti, G. Cortese, E. Follana, and M. Giordano, “A geometric Monte Carlo algorithm for the antiferromagnetic Ising model with ”topological” 118 BIBLIOGRAPHY term at Θ = π,” Nucl. Phys. B883 (2014) 656, arXiv:1312.6416 [hep-lat]. [131] A. Casher, J. B. Kogut, and L. Susskind, “Vacuum polarization and the absence of free quarks,” Phys. Rev. D10 (1974) 732. [132] J. B. Kogut and L. Susskind, “How to Solve the eta –¿ 3 pi Problem by Seizing the Vacuum,” Phys. Rev. D11 (1975) 3594. [133] C. J. Hamer, J. B. Kogut, D. P. Crewther, and M. M. Mazzolini, “The Massive Schwinger Model on a Lattice: Background Field, Chiral Symmetry and the String Tension,” Nucl. Phys. B208 (1982) 413. [134] T. Byrnes, P. Sriganesh, R. J. Bursill, and C. J. Hamer, “Density matrix renormalization group approach to the massive Schwinger model,” Phys. Rev. D66 (2002) 013002, arXiv:hep-lat/0202014 [hep-lat]. [135] B. Buyens, S. Montangero, J. Haegeman, F. Verstraete, and K. Van Acoleyen, “Finite-representation approximation of lattice gauge theories at the continuum limit with tensor networks,” Phys. Rev. D95 no. 9, (2017) 094509, arXiv:1702.08838 [hep-lat]. [136] Y. Shimizu and Y. Kuramashi, “Critical behavior of the lattice Schwinger model with a topological term at θ=πusing the Grassmann tensor renormalization group,” Phys. Rev. D90 no. 7, (2014) 074503, arXiv:1408.0897 [hep-lat]. [137] J. Schwinger, “Gauge invariance and mass. ii,” Phys. Rev. 128 (Dec, 1962) 2425. [138] N. Seiberg, “Topology in strong coupling,” Phys. Rev. Lett. 53 (Aug, 1984) 637–640. [139] U. J. Wiese, “Numerical Simulation of Lattice θVacua: The 2-dU(1) Gauge Theory as a Test Case,” Nucl. Phys. B318 (1989) 153. [140] H. Leutwyler and A. V. Smilga, “Spectrum of Dirac operator and role of winding number in QCD,” Phys. Rev. D46 (1992) 5607. [141] F. D. M. Haldane, “Continuum dynamics of the 1-D Heisenberg antiferromagnetic identification with the O(3) nonlinear sigma model,” Phys. Lett. A93 (1983) 464–468. BIBLIOGRAPHY 119 [142] D. G¨oschl, C. Gattringer, A. Lehmann, and C. Weis, “Simulation strategies for the massless lattice Schwinger model in the dual formulation,” Nucl. Phys. B924 (2017) 63, arXiv:1708.00649 [hep-lat]. [143] V. Azcoiti, G. di Carlo, and A. F. Grillo, “A New proposal for including dynamical fermions in lattice gauge theories: The Compact QED case,” Phys. Rev. Lett. 65 (1990) 2239. [144] V. Azcoiti, G. Di Carlo, A. Galante, A. F. Grillo, and V. Laliena, “The Schwinger model on the lattice in the microcanonical fermionic average approach,” Phys. Rev. D50 (1994) 6994, arXiv:hep-lat/9401032 [hep-lat]. [145] S. D¨urr, “Physics of η0with rooted staggered quarks,” Phys. Rev. D 85 (Jun, 2012) 114503. [146] S. D¨urr and C. Hoelbling, “Staggered versus overlap fermions: A study in the schwinger model with Nf= 0,1,2,” Phys. Rev. D 69 (Feb, 2004) 034503. [147] S. R. Coleman, “There are no Goldstone bosons in two-dimensions,” Commun. Math. Phys. 31 (1973) 259–264. [148] N. Mermin and H. Wagner, “Absence of ferromagnetism or antiferromagnetism in one-dimensional or two-dimensional isotropic Heisenberg models,” Phys. Rev. Lett. 17 (1966) 1133–1136. [149] V. Azcoiti, “Interplay between SU(Nf) chiral symmetry, U(1)Aaxial anomaly, and massless bosons,” Phys. Rev. D 100 no. 7, (2019) 074511, arXiv:1907.01872 [hep-lat]. [150] H. Georgi, “Automatic fine-tuning in the 2-flavor schwinger model,” arXiv:2007.15965 [hep-th]. [151] HPQCD Collaboration, C. McNeile, C. T. H. Davies, E. Follana, K. Hornbostel, and G. P. Lepage, “High-precision cand bmasses, and qcd coupling from current-current correlators in lattice and continuum qcd,” Phys. Rev. D 82 (Aug, 2010) 034512, arXiv:1004.4285. [152] HPQCD Collaboration, B. Chakraborty, C. T. H. Davies, B. Galloway, P. Knecht, J. Koponen, G. C. Donald, R. J. Dowdall, G. P. Lepage, and C. McNeile, “High-precision quark masses and qcd coupling from nf= 4 lattice qcd,” Phys. Rev. D 91 (Mar, 2015) 054508, arXiv:1408.4169 [hep-lat]. 120 BIBLIOGRAPHY [153] PACS-CS, S. Aoki et al., “Precise determination of the strong coupling constant in Nf= 2+1 lattice QCD with the Schrodinger functional scheme,” JHEP 10 (2009) 053, arXiv:0906.3906 [hep-lat]. [154] ETM, B. Blossier, P. Boucaud, M. Brinet, F. De Soto, V. Morenas, O. Pene, K. Petrov, and J. Rodriguez-Quintero, “High statistics determination of the strong coupling constant in Taylor scheme and its OPE Wilson coefficient from lattice QCD with a dynamical charm,” Phys. Rev. D 89 no. 1, (2014) 014507, arXiv:1310.3763 [hep-ph]. [155] K. Maltman, D. Leinweber, P. Moran, and A. Sternbeck, “The Realistic Lattice Determination of alpha(s)(M(Z)) Revisited,” Phys. Rev. D 78 (2008) 114504, arXiv:0807.2020 [hep-lat]. [156] A. Bazavov, N. Brambilla, I. Tormo, Xavier Garcia, P. Petreczky, J. Soto, and A. Vairo, “Determination of αsfrom the QCD static energy: An update,” Phys. Rev. D 90 no. 7, (2014) 074038, arXiv:1407.8437 [hep-ph]. [Erratum: Phys.Rev.D 101, 119902 (2020)]. [157] MILC Collaboration, A. Bazavov, C. Bernard, N. Brown, J. Komijani, C. DeTar, J. Foley, L. Levkova, S. Gottlieb, U. M. Heller, J. Laiho, R. L. Sugar, D. Toussaint, and R. S. Van de Water, “Gradient flow and scale setting on milc hisq ensembles,” Phys. Rev. D 93 (May, 2016) 094510, arXiv:1503.02769 [hep-lat]. [158] A. Sternbeck, The Infrared behavior of lattice QCD Green’s functions. PhD thesis, 9, 2006. arXiv:hep-lat/0609016. [159] R. Aouane, F. Burger, E.-M. Ilgenfritz, M. M¨uller-Preussker, and A. Sternbeck, “Landau gauge gluon and ghost propagators from lattice QCD with Nf=2 twisted mass fermions at finite temperature,” Phys. Rev. D87 no. 11, (2013) 114502, arXiv:1212.1102 [hep-lat]. [160] C. T. H. Davies, G. G. Batrouni, G. R. Katz, A. S. Kronfeld, G. P. Lepage, K. G. Wilson, P. Rossi, and B. Svetitsky, “Fourier acceleration in lattice gauge theories. i. landau gauge fixing,” Phys. Rev. D 37 (Mar, 1988) 1581–1588.