TL;DR
This paper develops a 3+1D numerical simulation method for the Glasma, a gluonic state in heavy ion collisions, capturing rapidity dependence and improving numerical stability over previous boost-invariant models.
Contribution
It introduces a novel 3+1D numerical approach for non-boost-invariant Glasma simulations, including a new gauge-covariant scheme with enhanced stability.
Findings
Energy density varies with rapidity, aligning with experimental data.
Pressure anisotropy remains large in the Glasma.
Rapidity profiles depend on collision energy.
Abstract
The Glasma is a gluonic state of matter which can be created in collisions of relativistic heavy ions and is a precursor to the quark-gluon plasma. The existence of this state is a prediction of the color glass condensate (CGC) effective theory. In many applications of the CGC framework, the boost invariant approximation is employed. It assumes that the longitudinal extent of the nuclei can be approximated as infinitesimally thin. Consequently, the Glasma produced from such a collision is boost invariant and can be effectively described in 2+1D. Therefore, observables of the boost invariant Glasma are by construction independent of rapidity. The main goal of this thesis is to develop a numerical method for the non-boost-invariant setting where nuclei are assumed to be thin, but of finite longitudinal extent. This is in conflict with a number of simplifications that are used in the…
| [GeV] | [fm] | [GeV] | [GeV] | |
|---|---|---|---|---|
| 200 | 0.0364 | 1.94 | 2.08 | 0.595 |
| 130 | 0.0560 | 1.84 | 2.11 | 0.551 |
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Code & Models
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
DISSERTATION
**Simulations of the Glasma in 3+1D
**
ausgeführt zum Zwecke der Erlangung des akademischen Grades eines Doktors der Technischen Wissenschaften unter der Leitung von
Priv.Doz. Dr.techn. Andreas Ipp
Institut für Theoretische Physik,
Technische Universität Wien, Österreich
eingereicht an der Technischen Universität Wien,
Fakultät für Physik
von
Dipl.-Ing. David Müller, BSc
Matrikelnummer 0825271
Antonie-Alt-Gasse 12/8/11
1100 Wien
[TABLE]
[TABLE]
Abstract
The Glasma is a gluonic state of matter which can be created in collisions of relativistic heavy ions. It only exists for a short period of time before it evolves into the quark-gluon plasma. The existence of the Glasma is a prediction of a first-principles classical effective theory of high energy quantum chromodynamics called the color glass condensate (CGC). In many analytical and numerical calculations within the CGC framework, the boost invariant approximation is employed. It assumes that the Lorentz-contracted longitudinal extent of the nuclei can be approximated as infinitesimally thin. Consequently, the Glasma produced from such a collision is boost invariant and can be effectively described in 2+1 dimensions. Observables of interest such as energy density, pressure or gluon occupation number of the boost invariant Glasma are by construction independent of rapidity.
The main topic of this thesis is to study how the assumption of infinitesimally thin nuclei can be relaxed. First, we discuss the properties of the CGC and Glasma by starting with the boost invariant case. The McLerran-Venugopalan (MV) model is used as a simple model for large heavy nuclei. The Yang-Mills equations, which govern the dynamics of the Glasma, generally cannot be solved analytically. Numerical solutions to these equations are therefore often the only reliable approach to studying the Glasma. We discuss the methods of real time lattice gauge theory, which are the usual approach to numerically solving the Yang-Mills equations in a gauge-covariant manner. Having established the standard tools used to describe the boost invariant Glasma, we focus on developing a numerical method for the non-boost-invariant setting where nuclei are assumed to be thin, but of finite longitudinal extent. This small change is in conflict with a number of simplifications and assumptions that are used in the boost invariant case. In particular, one has to describe the collisions in 3+1 dimensions in the laboratory or center-of-mass frame, compared to the co-moving frame of the traditional method. The change of frame forces the explicit inclusion of the color charges of nuclei. In numerical simulations this is achieved using the colored particle-in-cell method.
The new method is tested using a version of the MV model which includes a parameter for longitudinal thickness. It reproduces the boost invariant setting as a limiting case. Studying the pressure components of the Glasma, one finds that the Glasma in 3+1 dimensions does not differ much from the boost invariant case and that the pressure anisotropy remains large. On the other hand, one finds that the energy density of the Glasma depends on rapidity due to the explicit breaking of boost invariance. The width of the observed rapidity profiles is controlled by the thickness of the colliding nuclei. Using only a very simple model for nuclei, the profiles can be shown to agree with experimental data. If simulation parameters are not carefully chosen, the numerical scheme employed in the 3+1D method suffers from a numerical instability. To eliminate this instability, a completely new numerical scheme for real-time lattice gauge theory is developed. This new scheme is shown to be gauge-covariant and conserves the Gauss constraint even for large time steps.
Zusammenfassung
Das Glasma ist ein gluonischer Zustand, welcher in relativistischen Schwerionenkollisionen erzeugt werden kann und nur für sehr kurze Zeit existiert, bevor er in das Quark-Gluon-Plasma zerfällt. Die Existenz des Glasmas ist eine Vorhersage des Farbglaskondensats (engl. “color glass condensate” (CGC)). Das CGC ist eine klassische effektive Theorie, welche direkt aus der fundamentaleren Theorie der Quantenchromodynamik abgeleitet werden kann. In vielen analytischen und numerischen Rechnungen im Rahmen des CGCs kommt die boost-invariante Näherung zur Anwendung. In dieser Näherung nimmt man an, dass die dünne longitudinale Ausdehnung von Atomkernen (also der Lorentz-kontrahierte Durchmesser entlang der Strahlachse bzw. Bewegungsrichtung) infinitesimal ist. Eine Konsequenz dieser Näherung ist, dass das erzeugte Glasma invariant unter Lorentz-Boosts wird und somit effektiv in 2+1 Dimensionen beschrieben werden kann. Observablen, also im Prinzip beobachtbare Größen wie die Energiedichte, die Druckkomponenten und die Gluonenbesetzungszahl, sind dadurch per constructionem unabhängig von der Rapiditätskoordinate.
Das Thema dieser Dissertation ist eine neue Methode zu entwickeln, mit der man die Annahme der Boost-Invarianz lockern und umgehen kann. Ich beginne mit einer Diskussion über die physikalischen Eigenschaften des Glasmas und des CGCs im boost-invarianten Fall. Als einfaches Modell für große, schwere Atomkerne kommt das McLerran-Venugopalan-Modell (MV) zum Einsatz. Die Yang-Mills-Gleichungen, welche die Dynamik des Glasmas bestimmen, können im Allgemeinen nicht analytisch gelöst werden. Numerische Lösungsmethoden sind somit oft der einzige verlässliche Weg, um das Glasma zu untersuchen. Daher wird Echtzeit-Gittereichtheorie benötigt, welche die Standardmethode zur numerischen Lösung der Yang-Mills-Gleichungen darstellt. Nach dieser Einführung in die Standardwerkzeuge, welche verwendet werden, um das boost-invariante Glasma zu beschreiben, liegt der Fokus auf der Entwicklung einer numerischen Methode für den Fall, dass Boost-Invarianz nicht mehr gilt, also wenn man für relativistische Atomkerne eine kleine, aber endliche Ausdehnung entlang der Bewegungsrichtung annimmt. Diese kleine Änderung führt dazu, dass viele Annahmen und Vereinfachungen, die noch im boost-invarianten Fall verwendet werden konnten, nicht mehr gültig sind. Inbesondere muss die Kollision im Labor- bzw. Schwerpunktsystem in drei räumlichen Dimensionen beschrieben werden, anstelle des sich mit dem Glasma mitbewegenden Koordinatensystems im boost-invarianten Szenario. Dieser Koordinatensystemwechsel erfordert unter anderem, dass die Farbladungen der Atomkerne explizit berücksichtigt werden müssen. In numerischen Simulationen gelingt das mit der Particle-in-Cell-Methode, verallgemeinert auf Farbladungen.
Die neue numerische Methode wird getestet, indem Kollisionen von Kernen mit endlicher Dicke simuliert werden. Als Anfangsbedingung für diese Simulationen dient ein erweitertes MV-Modell, welches einen neuen Parameter für longitudinale Ausdehnung besitzt. Es wird gezeigt, dass die neue Methode das boost-invariante Szenario als Grenzfall beschreiben kann. Weiters wird auch die Anisotropie der Druckkomponenten des dreidimensionalen Glasmas untersucht, wobei nur wenige Unterschiede zum boost-invarianten, zweidimensionalen Glasma festgestellt werden können. Anders verhält es sich mit der Brechung der Boost-Invarianz: betrachtet man die Energiedichte des Glasmas im lokalen Ruhesystem, kann eine starke Rapiditätsabhängigkeit festgestellt werden, welche durch die Dicke der kollidierenden Kerne beeinflusst wird. Im Vergleich mit experimentellen Resultaten von echten Kollisionsexperimenten zeigt sich, dass mit diesem sehr einfachen Modell realistische Rapiditätsprofile erzeugt werden können. Die numerische Methode, welche für dreidimensionale Kollisionssimulationen entwickelt wurde, ist auf die Wahl der Simulationsparameter sensibel und kann in gewissen Fällen instabil werden. Die Ursache dieser numerischen Instabilität wurde identifiziert und eine Erweiterung der ursprünglichen Methode entwickelt, welche sich als stabil erweist. Es wird gezeigt, dass diese neue Methode eichkovariant ist und das Gaußsche Gesetz während der Simulation auch für große Zeitschritte erfüllt bleibt.
Preface
The results and methods presented in this thesis are largely based on the following published articles:
- •
D. Gelfand, A. Ipp and D. Müller, Phys. Rev. D 94, no. 1, 014020 (2016) [1605.07184] [1]
- •
A. Ipp and D. Müller, Phys. Lett. B 771, 74 (2017) [1703.00017] [2]
- •
A. Ipp and D. Müller, PoS EPS HEP2017, 176 (2017) [1710.01732] [3]
- •
A. Ipp and D. Müller, Eur. Phys. J. C 78, no. 11, 884 (2018) [1804.01995] [4]
Acknowledgements
There are a lot of people that I owe gratitude toward for their direct or indirect involvement with my PhD studies and this thesis.
First and foremost, I would like to thank my supervisor Andreas Ipp for giving me the opportunity to work on this project and providing guidance during my PhD. I thank Daniil Gelfand for collaboration during the early development of the project. I also wish to thank Jean-Paul Blaizot, François Gelis, Edmond Iancu, Aleksi Kurkela and Sören Schlichting for interesting discussions about the CGC, the Glasma and heavy-ion collisions in general. My colleagues Patrick Kappl and Axel Polaczek deserve special thanks for many interesting discussions about technical details regarding numerical algorithms, programming and lattice gauge theory. I thank Kayran Schmidt for carefully reading this thesis. I thank Tuomas Lappi for two very productive research stays at the University of Jyväskylä. I also thank Elena Petreska and Carlos A. Salgado for inviting me to the University of Santiago de Compostela and enlightening discussions.
Furthermore, the scientific and personal support of my friend and colleague Alexander Haber has been invaluable throughout the past ten years of studying physics. The same applies to Georg Harrer, Sebastian Schönhuber, Alexander Soloviev, Wolfgang Steiger, Franz-Stephan Strobl and many other friends I’ve made in the past years. I would also like to thank Gerald Hattensauer for many supportive discussions. I’d like to give special thanks to Frederic Brünner and Christian Ecker, who I’ve had the pleasure to share an office with in the last five years. My research stays in Finland have been very enjoyable due to Kirill Boguslavski, Risto Paatelainen, Jarkko Peuron and Andrecia Ramnath. Over the past years I’ve also received valuable guidance and advice from my former supervisor Andreas Schmitt.
Finally, I would like to thank my family, in particular my parents Albert and Irene Müller, my uncle Paul Müller, and my partner Alina G. Dragomir.
Contents
-
2.5.3 Rapidity dependence and finite longitudinal extent of nuclei
-
5.1 Numerical Cherenkov instability and numerical dispersion
Chapter 1 Introduction
Heavy-ion collision experiments such as the ones at the Relativistic Heavy Ion Collider (RHIC) and the Large Hadron Collider (LHC) provide fascinating insights into the properties of strongly interacting nuclear matter under extreme conditions. The main phenomenon of interest in these collisions is the creation and evolution of the quark-gluon plasma (QGP), an exotic state of matter, where the constituents of nuclei, neutrons and protons, come apart and break into their fundamental building blocks, namely quarks and gluons. The study of the QGP is of crucial importance for testing the theory of quantum chromodynamics (QCD) by explaining and predicting experimental data. Moreover, the QGP created in relativistic nuclear collisions can be used as a model for the early universe immediately after the big bang and therefore also has implications for cosmology [5].
Although the basic equations of QCD, which describe the interactions between quarks and gluons, are well established, the theoretical description of heavy-ion collisions in terms of the full theory is only partially tractable. This is due to the phenomenon of asymptotic freedom: the strong coupling constant becomes large at low momenta and consequently perturbative techniques are bound to fail. The evolution of the “fireball” created in heavy-ion collisions is therefore split into various stages with different appropriate models used for each stage which approximate the underlying QCD processes [6]. Roughly speaking, the three main stages are the pre-equilibrium, the equilibrium and the freeze-out. The pre-equilibrium describes the earliest stage, starting directly after the collision until the evolving matter is in thermal equilibrium, which is a process known as thermalization. The fireball then evolves as a QGP in thermal equilibrium until the freeze-out, where quarks and gluons recombine and form a gas of hadrons. This gas stops interacting after some time and the free streaming hadrons travel towards the detectors.
One of the most striking results to come out of heavy-ion collision experiments is that the QGP in the equilibrium stage almost behaves like an ideal fluid [7, 8, 9, 10]. The time evolution of the QGP can be described in large parts by relativistic viscous hydrodynamics [11, 12, 13], but this successful description only applies to the evolution of the QGP fireball itself and not to the earliest stages of the collision. Relativistic hydrodynamical simulations therefore require initial conditions from the pre-equilibrium stage for which there are different types of models. Popular choices are variants of the phenomenological MC-Glauber model [14], but a more sophisticated approach to the initial state of nuclear collisions, developed over the past two and a half decades, is the color glass condensate (CGC) [15, 16, 17]. It is a classical effective theory for high energy QCD and provides a first-principles model of the early stages of relativistic heavy-ion collisions.
The main idea behind the CGC is a separation of scales: the hard constituents of a relativistic nucleus, i.e. partons which carry most of total momentum such as valence quarks, are described as highly Lorentz-contracted, thin sheets of classical color charge. These fast color charges generate a highly occupied color field, which represents the soft partons of the nucleus, namely mostly gluons at lower momenta. The longitudinal momentum cutoff, at which one performs the separation into soft and hard partons, is entirely arbitrary and by requiring that observables do not depend on this artificial cutoff one can obtain a set of renormalization group equations known as the JIMWLK equations [18, 19]. The CGC is therefore an effective description of high energy nuclei in terms of classical color fields and color currents, whose dynamics are governed by classical Yang-Mills theory. In the CGC model the result of a collision of two nuclei is a state called Glasma [20] (a combination of glass and plasma), which is a precursor to the QGP and can also be treated in terms of classical field theory. The combination of relativistic hydrodynamics with CGC/Glasma initial conditions from Yang-Mills simulations (CGC+Hydro) using models such as MC-KLN [21, 22] or IP-Glasma [23, 24] has been highly successful in explaining observed phenomena in high energy heavy-ion and proton-nucleus collisions: total particle multiplicity, azimuthal anisotropy [25] and higher flow coefficients [26] observed at the LHC can all be understood in terms of CGC/Glasma and hydrodynamics. Simulations of proton ion collisions even reveal details about the sub-nucleonic structure of the proton [27, 28]. The CGC also provides an explanation for long-range rapidity correlations (the ridge phenomenon) [29, 30].
Despite the obvious success of the CGC+Hydro approach, the description of the initial state is incomplete with regard to rapidity dependence and longitudinal dynamics. In the limit of very high collision energies the color field of a nucleus is usually assumed to be an infinitesimally thin shock wave, which makes it possible to derive analytic expressions for the Glasma fields directly after the collision [31]. As a result of this approximation, the Glasma behaves in a boost invariant fashion, which means that observables are by construction independent of space-time rapidity and reduce the system from 3+1 dimensions to effectively only 2+1 dimensions. Any non-trivial rapidity dependence of observables is usually obtained from independent boost invariant simulations at different values of rapidity [32, 33, 34]. Although 3+1 dimensional Yang-Mills simulations of the Glasma exist, they are mostly performed in the context of studying non-Abelian plasma instabilities [35, 36, 37, 38, 39, 40, 41, 42] and not yet for obtaining realistic initial conditions for hydrodynamic simulations with correct rapidity dependence.
In this thesis we develop our new approach for simulations of the pre-equilibrium stage of 3+1 dimensional collisions in the CGC framework, which is able to go beyond the boost-invariant approximation. This approach is distinguished from the conventional boost invariant scenario by working directly in the laboratory (or center-of-mass) frame of the colliding nuclei using Cartesian coordinates. Instead of treating the nuclei as infinitesimally thin shock waves, we allow for a finite Lorentz-contracted longitudinal extent along the beam axis, which is proportional to , where is the nuclear radius and is the Lorentz factor. The numerical method that is used to solve the classical Yang-Mills equations in this frame takes inspiration from the colored particle-in-cell approach (CPIC) [43, 44] for simulating classical non-Abelian plasmas. This new approach to simulating collisions not only covers the evolution of the Glasma but also the collision event itself and therefore gives a genuine 3+1 dimensional picture of the early stages of heavy-ion collisions in the CGC/Glasma framework.
This thesis is organized as follows: chapter 2 is a short introduction to the CGC and the Glasma in the boost invariant scenario. We discuss the McLerran-Venugopalan model of large nuclei as a simple model for relativistic heavy ions and develop the methods used to perform the classical time evolution of the Glasma numerically on a lattice. The main physical phenomena and some observables in the Glasma are also mentioned. In chapter 3 the numerical method for simulating 3+1 dimensional collisions in the laboratory frame is developed. In chapter 4 this numerical method is applied to a specific extension of the McLerran-Venugopalan model in 3+1 dimensions and some numerical results are presented. In particular, in section 4.4 the rapidity profile of the Glasma energy density is computed, which shows significant deviations from the boost-invariant case. Finally, in chapter 5 a numerical improvement is developed that allows for more stable and accurate simulations by employing a semi-implicit solving method. Longer, more tedious derivations and calculations are summarized in the appendix. For physical units, sign conventions and other special notation used in this thesis, see appendix A.
Chapter 2 Aspects of the boost invariant Glasma
This chapter is an introduction to the methods and key phenomena associated with the boost invariant Glasma, i.e. the Glasma created in an ultrarelativistic heavy-ion collision where the colliding nuclei are assumed to be infinitely thin. In the following, we will focus on the necessary fundamentals in order to formulate the 3+1 dimensional model of the Glasma presented in the later chapters. For a more general overview see e.g. [15, 16, 17, 20, 45].
First, we introduce the classical field theory description of high energy nuclei using a very simple model, namely the McLerran-Venugopalan (MV) model. Then the boost invariant collision of two nuclei is discussed within this model. We derive the initial state and the equations of motion governing the Glasma after the collision. In order to investigate the time evolution of the Glasma, the equations of motion have to be solved numerically, which can be accomplished using real-time lattice gauge theory. Using these methods it is possible to study quantities like the energy density and the pressure components of the Glasma. Finally, state-of-the-art extensions to the MV model are discussed which allow for correct modeling of realistic nuclei and we discuss how to go beyond the boost invariant approximation.
2.1 The McLerran-Venugopalan model
As discussed in chapter 1, nuclei at very high energies can be described within the CGC framework, which is a high energy effective theory for QCD and enables an effectively classical description of nuclei in terms of classical color fields and color charges. However, even before the full development of the CGC framework as a genuine effective theory, McLerran and Venugopalan already proposed a classical model for high energy nuclei [46, 47].
They argued that the valence quarks, which carry most of the total momentum of the nucleus (referred to as “hard” degrees of freedom), are to be treated as static, recoilless color charges whose dynamics are “frozen” due to time dilation. Furthermore they proposed that the color charge density of the valence quarks provides an energy scale much larger than the QCD scale . This implies that the Yang-Mills coupling constant can be considered to be weak, which justifies a classical (or eikonal) approximation. A priori the exact positions and charges of the valence quarks are unknown and therefore the classical color charges are considered to be fluctuating random variables. The valence quarks are therefore represented by a random classical color charge density (or more generally a classical color current ) whose probability distribution is specified by a probability functional .
On the other hand, the gluons, which carry only a fraction of the total momentum (the “soft” degrees of freedom), are considered to be dynamic. McLerran and Venugopalan realized that because of the large number of gluons in high energy nuclei, quantum mechanical effects should be, as a first approximation, neglected111More precisely, the gluon field is highly occupied and forms a coherent state, which is essentially a classical field state.. The soft gluons are therefore represented by a classical color field . Since the classical Yang-Mills equations must hold, the color field of the gluons is fixed by the color charge density of the valence quarks, such that the Gauss constraint is fulfilled.
Using these assumptions we can start solving the classical problem. The equations of motion are most easily solved by employing light cone coordinates , where and are the laboratory frame coordinates. When using light cone coordinates, Latin indices as in are reserved for transverse coordinates, i.e. . We then consider (without loss of generality) the color current of a nucleus moving in the positive direction (i.e. ) at the speed of light. The only relevant component is then which we associate with the color charge density . We can drop the dependency of the current because we assume the valence quarks to be static in . We are then left with
[TABLE]
where are the coordinates in the transverse plane spanned by and , and are the generators of the gauge group in the fundamental representation. The fast moving nucleus that we are describing is highly Lorentz contracted in the direction which tells us that the support along the direction must be very thin. In the ultrarelativistic limit the longitudinal support becomes infinitesimal and the color current is proportional to :
[TABLE]
For the present discussion however, we will keep the more general form eq. 2.1 with the color current being strongly peaked around . A schematic diagram is shown in fig. 2.1.
The next step is to solve the classical Yang-Mills equations (see section A.2 for the conventions used in this thesis)
[TABLE]
in the presence of the external color current . In order to make progress we must choose a gauge condition and an appropriate ansatz for the color field. A useful starting point is the covariant (or Lorenz) gauge condition
[TABLE]
with the ansatz
[TABLE]
The only remaining field component is and the gauge condition forces us to drop the dependency of the color field due to
[TABLE]
We also see that the ansatz is compatible with non-Abelian charge conservation, i.e.
[TABLE]
due to . Inserting the ansatz into the field strength tensor
[TABLE]
we find that the only non-zero components are given simply by
[TABLE]
Plugging this into the Yang-Mills equations yields
[TABLE]
which reduces to a Poisson equation in the transverse plane
[TABLE]
This equation can be readily solved by inverting the two-dimensional Laplace operator , for instance using the Fourier transform. Defining the partial transform
[TABLE]
and its inverse
[TABLE]
we find the solution
[TABLE]
Note that the color field inherits its longitudinal support and shape from the color charge density . It is also instructive to analyze the field strength of the color field: switching back to laboratory frame coordinates and we find
[TABLE]
because . The only non-zero field strength components are then and with , which means that the nucleus only has transverse color-electric and color-magnetic fields
[TABLE]
which are orthogonal to each other and have the same magnitude. Note that refers to the three-dimensional Levi-Civita symbol, while is the two-dimensional Levi-Civita symbol in the transverse plane. It turns out that the solution for a single propagating nucleus is completely analogous to the electromagnetic field of an electric charge moving at the speed of light, i.e. ultrarelativistic Liénard-Wiechert potentials. Due to the ansatz and a clever choice of gauge given by eq. 2.4, all non-linear terms of the Yang-Mills equations can be ignored and the resulting solution ends up being remarkably simple.
The choice of covariant gauge becomes inconvenient when discussing collisions of nuclei, where it turns out that light cone (LC) gauge is much better suited. Our goal therefore is to find the gauge transformation acting on the gauge field via
[TABLE]
such that vanishes. This requirement leads us to
[TABLE]
The solution to this equation is given by the path-ordered exponential
[TABLE]
Note that we use the “left means later” convention for path ordering, i.e. if , i.e. comes after along the path. Equation 2.20 implies that the Wilson line at the asymptotic boundary is a unit matrix.
We identify the gauge transformation with the lightlike Wilson line starting at and ending at . Due to we still have like in covariant gauge. On the other hand, the transverse components of the gauge field are given by
[TABLE]
The color current also has to be transformed accordingly and is given by
[TABLE]
where the subscript “LC” is used to differentiate the LC gauge current from the covariant gauge current . For ultrarelativistic nuclei we can derive a simpler relationship between the transverse gauge fields and the LC gauge current : a -shaped charge density as in eq. 2.2 implies that the transverse gauge field has the form of a Heaviside step function given by
[TABLE]
where is the asymptotic Wilson line
[TABLE]
The ultrarelativistic LC current is given by
[TABLE]
Inserting the above current and eq. 2.23 into the Yang-Mills equations (2.3) yields
[TABLE]
which results in the relation
[TABLE]
We see that the space-time picture of the LC gauge solution eq. 2.23 is very different compared to the covariant gauge case: the transverse gauge fields are non-zero for and extend to . However, this is merely an artifact of the LC gauge condition as the transverse gauge fields are pure gauge, i.e. there exists a gauge condition where the field becomes zero (namely the covariant gauge). In contrast, the actual color-electric and color-magnetic field strengths and the color current are still concentrated around .
Now that the relationship between the color charge density and the color field is established, one has to specify what the color charge distribution of a large nucleus looks like. In their original formulation McLerran and Venugopalan assumed that the charge density is -shaped as in eq. 2.2 and proposed a simple Gaussian probability distribution for the color charge density of the valence quarks. The distribution is defined by the charge density one- and two-point functions
[TABLE]
The one-point function eq. 2.29 guarantees that the nucleus is on average color neutral. The two-point function eq. 2.30 fixes the average color charge fluctuation around zero, where is the phenomenological MV parameter given in units of energy or inverse length. For a large nucleus with nucleons they estimated from the average density of valence quarks that (see section A.1 for the unit conventions used in this thesis)
[TABLE]
For a gold nucleus with this estimate gives . The MV model does not assume a finite transverse extent of the nucleus. Instead, it approximates very large nuclei as infinitely thin, but transversely infinite walls of color charge. Within the transverse plane charges at different points are completely uncorrelated due to the term in eq. 2.30. On average, the MV model exhibits translational and rotational invariance in the transverse plane. The MV model thus gives only a very crude approximation of a realistic nucleus, but due to its simplicity with only one dimensionful parameter and its high symmetry, many otherwise complicated calculations can be performed analytically.
Using the one- and two-point functions eqs. (2.29) and (2.30), together with the assumption that the random color charges obey a Gaussian distribution, one can define the probability functional as
[TABLE]
where is a normalization constant. The probability functional is used to define expectation values of arbitrary observables via a functional integral over all charge density configurations
[TABLE]
where it is implied that is the color field associated with the charge density . The probability functional , eqs. (2.29) and (2.30) are invariant under arbitrary gauge transformations :
[TABLE]
It was later realized [48] that for certain derivations the -approximation of eq. 2.2 can be problematic. A more rigorous approach is to first regularize the -peak using eq. 2.1 for intermittent calculation steps and then only performing the ultrarelativistic limit in the final results. A generalization of the original MV model (see eqs. (2.29) and (2.30)) with finite longitudinal support is given by [48]
[TABLE]
or equivalently
[TABLE]
Here the function defines the average color charge fluctuation in the nucleus and the longitudinal shape along the direction.
The physical picture of the generalized MV model is slightly different compared to the original MV model, see fig. 2.2: for a nucleus with thin but finite longitudinal support, one can think of the nucleus as a stack of uncorrelated, infinitely thin sheets of color charge due to the additional term in the charge density two-point function eq. 2.36. On the other hand, the original MV model collapses this stack of color sheets into a single, infinitesimal sheet of color charge. Differentiating between the two formulations of the MV model is important because in [49] it was found that the ultrarelativistic limit of the generalized MV model is actually not just simply obtained from the replacement due to subtleties of the path ordering of color charges. Even in the case of infinitesimal width along , the asymptotic Wilson line eq. 2.25 “remembers” the ordering of the uncorrelated sheets of color charge.
In [49] a regularization of eq. 2.25 is presented that can be used for both numerical and analytical calculations: as described above, one imagines the nucleus as a stack of separate, infinitesimally thin, uncorrelated sheets of color charge. The charge density two-point function can then be approximated as
[TABLE]
where the indices are introduced to denote different sheets. This regularization effectively discretizes the term in eq. 2.36. For each sheet one then solves the Poisson equation
[TABLE]
The lightlike Wilson line eq. 2.25 is given by
[TABLE]
Taking the limit one obtains the ultrarelativistic limit of the generalized MV model. For one is left with the original formulation of the MV model. We will therefore refer to this limiting case as the “single color sheet approximation”. Obviously, results generally depend on due to the nature of path ordering. This fact is relevant when one tries to relate the MV parameter to another central quantity, the saturation momentum .
One way to define the saturation momentum within the CGC framework is via the inverse correlation length of the Wilson line two-point function in the fundamental representation
[TABLE]
where the fundamental representation Wilson lines are given by eq. 2.25. Translational and rotational invariance in the transverse plane implies that the two-point functions only depend on the distance . An intuitive physical picture associated with eq. 2.41 is a recoilless quark-antiquark pair “probing” the nucleus at different transverse coordinates and at the speed of light. Due to the quarks being recoilless (or eikonal), they pass through the nucleus without changing their lightlike trajectories along . However, the quarks are not unaffected by the nucleus: as they pass the color field, their color charges rotate in accordance with non-Abelian charge conservation, i.e. the continuity equation has to hold. After passing through, one can compare the rotated color charges of the quarks: if they are very close, meaning that the transverse separation is much smaller than any transverse length scale of the nucleus (e.g. ), the quarks experience the same color rotation and thus . At very large separation, the color fields and Wilson lines and are uncorrelated: the color charges of the quarks are completely random and thus . The Wilson line two-point functions therefore define a characteristic length scale at which an initially correlated probing recoilless quark pair becomes de-correlated. The inverse of this length is associated with the characteristic transverse momentum scale . Specifically, one defines
[TABLE]
There exist multiple, slightly different definitions of the saturation momentum (see [50] for details), involving for instance the adjoint representation Wilson line two-point function instead of eq. 2.41. Moreover, the above definition involves the two-point function in coordinate space, but it is also possible to define using the momentum space representation of eq. 2.41. On dimensional and parametrical grounds should be roughly . One factor of is due to eq. 2.36 and the second one is introduced in the Wilson line eq. 2.25. In the color sheet regularization the ratio depends on and has to be determined numerically [50]. The limiting case , however, can be performed analytically (see e.g. [48] and [45] for a detailed derivation). Ignoring all of these complications, it will be sufficient for the purpose of this thesis to assume that the ratio is roughly close to .
Another interesting two-point function to characterize the MV model is the correlator of the gauge field given by
[TABLE]
which can be directly related to the charge density correlator eq. 2.36 via eq. 2.14. We then find that
[TABLE]
The momentum integral in the gauge field two-point function of the MV model turns out to be infrared divergent and has to be regularized. A popular way of curing this divergence is to introduce an infrared regulator by replacing eq. 2.14 with
[TABLE]
where has units of energy. Effectively, suppresses the long-range behavior of the gauge field , which can be illustrated by choosing to be a point charge
[TABLE]
which yields the gauge field
[TABLE]
where are modified Bessel functions of the second kind. At large distances the gauge field falls of exponentially
[TABLE]
The regulator thus plays the role of a screening mass. A possible physical interpretation one can assign to this infrared regulator is that its inverse mimics the confinement radius (roughly ) by imposing color neutrality at distances larger than the nucleon size. Thus, a more careful way of introducing infrared regulation is to change the color charge density via
[TABLE]
without modifying eq. 2.14. This circumvents the problem of adding a gauge-invariance breaking term in order to obtain eq. 2.45. Moreover, it is now also clear that the color charge density is globally color neutral because the zero mode of vanishes. Finally, we can compute the (now finite) two-point function of the gauge field
[TABLE]
We will use this two-point function in the introduction of chapter 4. For now it, is important to know that the MV model contains an infrared divergence that has to be manually regularized. There are other ways to cure the divergence by imposing color neutrality: one can require only global color neutrality by eliminating the zero mode of , which is equivalent to subtracting the monopole contribution of the charge density. However, results will then depend on the size of the system. Since the MV model has no notion of finite size in the transverse plane, one is forced to manually choose some kind of system size. In practice the system size is fixed by the transverse size of the lattice in numerical simulations. Optionally one can also choose to not only eliminate the monopole, but also subtract the dipole contribution of , which again depends on system size and changes the infrared behavior of in a slightly different manner. For an extended discussion on how to impose color neutrality see [51]. In this thesis we use eq. 2.49 not only for convenience, but also because this way of regularization has been established in more elaborate models of nuclei such as IP-Glasma [23, 24].
Even though the results discussed in this section apply to the MV model for nuclei moving in the positive direction (along ), completely analogous calculations can be done for the opposite direction along . In this case the color current and the gauge field in covariant gauge are given by
[TABLE]
In LC gauge, now referring to the gauge in which , we find
[TABLE]
with the lightlike Wilson line along given by
[TABLE]
The charge density correlator of the MV model reads
[TABLE]
2.2 The Glasma initial conditions
Equipped with the single nucleus solutions of the classical field equations from the last section, we can start investigating ultrarelativistic collisions and the kinds of color fields produced in these collisions. We denote the two nuclei by “A” and “B”. The combined color current of both nuclei is given by
[TABLE]
where nucleus “A” moves along the axis and “B” along . In the ultrarelativistic limit this color current defines the sharp boundary of the light cone in fig. 2.3, which separates the Minkowski diagram into four distinct regions I - IV. Due to causality the regions I, II and III are unaffected by the collision and therefore we can use the single nuclei solutions from the last section. However, non-linear interaction terms in the Yang-Mills equations, result in a non-trivial color field in region IV, i.e. the future light cone. This color field is the Glasma. Our goal is to find this solution to the Yang-Mills equations in the presence of these moving color charges. Formally, the solution in region IV is a functional of the two charge densities, i.e. . Using this color field we can compute arbitrary observables via functional integration
[TABLE]
where naturally one has to integrate over both color charge densities.
An important assumption is that the color charges do not suffer any recoil during the collision, i.e. they do not lose any longitudinal momentum and stay on their fixed, lightlike trajectories. Consequently, both color currents must be separately conserved:
[TABLE]
Even if the trajectories are unaffected, non-Abelian charge conservation still leaves the possibility that the color charges experience color rotation due to the presence of the field, leading to non-trivial () dependence of (). In order to avoid this complication we choose a gauge condition such that () vanishes along the () axis. A suitable choice is the Fock-Schwinger gauge condition
[TABLE]
which we use in the future light cone . In regions I - III we use the LC gauge condition , which can be smoothly connected to Fock-Schwinger gauge at the boundaries , and , .
The combined color current eq. 2.56 features an important symmetry, namely invariance under longitudinal boosts, i.e. boost invariance. Specifically, this means that under the Lorentz transformation
[TABLE]
the current remains unchanged in its functional form:
[TABLE]
Since the color current is the only input that we have for solving the Yang-Mills equations (apart from boundary conditions and choice of gauge), the solution and observables should also reflect this symmetry. One therefore introduces a coordinate system in which boost invariance becomes manifest, namely the , coordinate set (see section B.1). We define
[TABLE]
with proper time and space-time rapidity to parametrize the future light cone (, ) with a boundary formed by . Inverting the above definitions yields
[TABLE]
Applying the coordinate transformation to the gauge field we find (see section B.1)
[TABLE]
We can therefore make a manifestly boost invariant ansatz for the color field by “stiching” together separate functions in the four regions of the Minkowski diagram [52]:
[TABLE]
where we suppressed all dependencies on , rendering the ansatz boost invariant. In principle it would be possible to perform an -dependent gauge transformation without breaking boost invariance in gauge invariant observables, but for the sake of simplicity we also require this symmetry at the level of gauge fields. The transverse gauge fields of the single nuclei fulfill
[TABLE]
and are pure gauge in region (II) (in the case of nucleus A) and (III) (nucleus B). Due to eq. 2.67 the Fock-Schwinger gauge condition renders the temporal component zero. Plugging this ansatz into the Yang-Mills equations one obtains a few worrying terms such as and due to partial derivatives acting on the Heaviside functions. Setting all coefficients of these problematic terms to zero one finds matching conditions at the boundary [52]:
[TABLE]
together with . These matching conditions form the starting point for the time evolution of the color field in the future light cone, i.e. the Glasma field.
The longitudinal electric and magnetic field in the laboratory frame is related to the field strength tensor in the co-moving frame via (see eqs. B.17 to B.20)
[TABLE]
At this yields
[TABLE]
where and are the Kronecker delta and the Levi-Civita symbol in two dimensions. The magnetic contribution of the single nucleus fields vanishes due to them being pure gauge as specified by eq. 2.24. On the other hand, the transverse electric and magnetic field components vanish for . The collision of two nuclei therefore produces fields that are initially purely longitudinal. This leads to the following physical picture directly after the collision [20]: as the purely transverse color fields of the nuclei move away from the collision region, longitudinal electric and magnetic flux tubes span between the two nuclei. The size of these Glasma flux tubes corresponds to the size of correlated domains in the color fields of the nuclei, i.e. roughly [53]. After the initial formation, the flux tubes quickly evolve and start to expand in the transverse directions. In the next section we will introduce methods to study the time evolution of this system numerically.
2.3 Numerical time evolution
2.3.1 Boost invariant equations of motion
The gluonic matter produced in an ultrarelativistic heavy-ion collision is described as a boost invariant classical color field in the forward light cone. In the last section we sketched the derivation of the initial conditions, i.e. how this color field is determined by the color fields of the nuclei at the boundary . For the time evolution is governed by the source-free Yang-Mills equations
[TABLE]
However, since it is our goal to work in coordinates, it is reasonable to start from the Yang-Mills action in flat coordinates
[TABLE]
and perform the coordinate transformation there. Furthermore, we implement temporal (or Fock-Schwinger) gauge in order to be consistent with the initial conditions at and restrict ourselves to manifestly boost invariant gauge fields and . The boost invariant action then reads (see section B.1 for a derivation)
[TABLE]
supplemented by the Gauss constraint
[TABLE]
The canonical momenta and are given by
[TABLE]
Varying the action with respect to the gauge fields, one obtains the equations of motion (EOM)
[TABLE]
In these variables the initial conditions are given by
[TABLE]
In general, there exist no closed-form solutions to these equations due to the non-linear self-interaction terms. A number of approximate solution approaches have been proposed: for example, one can solve the boost invariant Yang-Mills equations by expanding in terms of the color fields of the nuclei [31]. This is only strictly correct if the color fields can be considered weak. Within the color glass condensate framework it is predicted that nuclei have non-perturbatively large fields . Thus, a series expansion is bound to fail at least for nucleus-nucleus collisions, but is valid for proton-nucleus collisions. Remarkably, one still obtains good agreement with numerical results by approximating the Glasma as weak and non-interacting, but keeping all orders of the color field in the Glasma initial conditions at [54, 55]. A very intriguing semi-analytic approach is to perform a series expansion of the gauge fields in powers of proper time [56, 57, 58, 59]. A downside is that due to slow convergence of the series, very high orders in are required to obtain reasonable agreement with known numerical results. A more general approach, which is also easily extendable to models of nuclei more complicated than the MV model and which in principle works up to arbitrary , is to solve the equations numerically on a lattice.
2.3.2 Real-time lattice gauge theory
Real-time lattice gauge theory is a numerical method used to solve the classical Yang-Mills equations on the lattice in a gauge covariant manner. It is based on lattice gauge theory (see [60, 61, 62] for introductions to the topic), but instead of solving path integrals in Euclidean space, one formulates a discretized, gauge invariant action in Minkowski space to obtain (discretized) equations of motion. The main strength of this method is that it retains a notion of gauge covariance even for the discretized system. We will quickly discuss the ideas behind this method.
A naive approach to solving the Yang-Mills equations could be the following: starting from the continuous action
[TABLE]
we make a choice on the degrees of freedom to use. An obvious choice are the gauge fields and their canonical momenta . Then, we replace Minkowski space by a regular hypercubic lattice with lattice spacings (see section A.3):
[TABLE]
where with unit vectors . In this step one replaces all continuous derivatives (such as in ) with finite difference expressions accurate up to some order in . In principle, this procedure yields a discretized action, and upon variation, discrete equations that have a correct continuum limit. Unfortunately, doing so one destroys the most important property of a gauge theory, namely local gauge invariance. The main issue is that local gauge transformations
[TABLE]
involve derivatives acting on the gauge transformation . Consequently, the discretized action and the equations that follow from variation are only gauge invariant up to some order in and gauge invariance too is only approximate. This is particularly troublesome when solving path integrals, which involve integrating over non-smooth field configurations, but breaking gauge invariance is also undesirable when merely solving classical equations.
However, there is a way to remedy this problem. Instead of making the obvious choice of , we use a different set of degrees of freedom. In section 2.1, we introduced Wilson lines: given a path starting at and ending at , the Wilson line is defined by the path ordered exponential
[TABLE]
Applying a local gauge transformation to , one can show that the Wilson line transforms non-locally at the start and end points
[TABLE]
In particular, if the path is closed (i.e. a “loop”), then one talks about Wilson loops. In this special case transforms locally, i.e.
[TABLE]
Taking the trace of a Wilson loop yields a gauge invariant quantity
[TABLE]
On a lattice, for instance defined by eq. 2.91, the shortest possible Wilson lines are called gauge links, which connect nearest neighbors on the lattice. For example connects the points and . In lattice gauge theory, is actually the anti-path-ordered Wilson line (denoted by )
[TABLE]
which transforms as
[TABLE]
or shorter (see section A.3 for an introduction to this shorthand notation)
[TABLE]
Close to the continuum limit we can use the approximation
[TABLE]
where the gauge field is placed at the mid-point of the path. The smallest possible Wilson loops that we can form on the lattice are the rectangular plaquettes ( Wilson loops)
[TABLE]
where we define gauge links with negative directions via . Taking the continuum limit of the plaquette one finds
[TABLE]
and consequently
[TABLE]
which is a gauge invariant expression. Note that the error term is one order higher than one would expect just from the continuum limit of the plaquette. In this expansion the terms proportional to cancel, because the term amounts to taking the real value of the plaquette, which eliminates precisely all terms (see subsection 2.3.2 of [61], or [60, 63]). Apart from the correct continuum limit, the most important property of eq. 2.103 is its exact invariance under gauge transformations eq. 2.99. Therefore, by reducing the space of local gauge transformations to so-called lattice gauge transformations only defined at lattice points , one can retain a notion of exact gauge invariance for the discretized system. This allows us to formulate a lattice gauge invariant action in Minkowski space. We start by splitting the action into its electric and magnetic part
[TABLE]
with
[TABLE]
where . We then use the gauge invariant combination of plaquettes eq. 2.103 to approximate
[TABLE]
where is the discrete space-time volume. The standard Wilson gauge action [64] is therefore given by
[TABLE]
Since we have discretized all of Minkowski space as a lattice, including the time direction , the above action is merely a function of gauge links (instead of a functional). Extremizing with respect to all gauge links yields discretized equations of motion and constraints.
However, coming back to solving the time evolution of the boost invariant Glasma, the Wilson action section 2.3.2 in its original formulation is not particularly useful. It is formulated in the frame, instead of the more appropriate frame in which boost invariance becomes manifest. A rigorous approach to derive the correct action (see e.g. [65, 66]) is to start from the discretized Wilson action section 2.3.2 and take the partial continuum limit . One then ends up with an action where the transverse plane is discretized as a lattice, while and are continuous. It is then possible to perform the coordinate transformation to the system and implement boost invariance by assuming all fields to be independent of .
Alternatively, one can also obtain the boost-invariant discretized action in temporal gauge by taking a few shortcuts (see section B.2 for more detailed derivation): the idea is to start directly from the continuous boost invariant case given by eq. 2.80 and replace the transverse fields with transverse gauge links . The term can be discretized using eq. 2.103 in a straightforward manner. Introducing the gauge-covariant forward and backward finite differences
[TABLE]
we can formulate the discretized action for the boost invariant system:
[TABLE]
with the canonical momenta given by
[TABLE]
The relation for guarantees that the gauge links remain special unitary matrices throughout the time evolution.
As we have already implemented temporal gauge, it is not possible to find the Gauss constraint through variation of the above action. Nevertheless, an educated guess yields
[TABLE]
If we would have started from an action continuous in without fixing any gauge condition (as detailed in [66]), the above equation would be the result of varying with respect to . Performing the variation with respect to is straightforward and yields
[TABLE]
where (no sum implied) is the second order gauge-covariant finite difference. Varying with respect to using (see section C.2)
[TABLE]
yields the equations of motion for the transverse components
[TABLE]
with the parallel transported field and the anti-hermitian, traceless part (see eq. A.20 of section A.2)
[TABLE]
Even though a bit cumbersome, it is straight-forward to prove that the equations of motion conserve the Gauss constraint eq. 2.115.
The final step is to discretize the proper time coordinate as well. Using a finite time step one can formulate a leapfrog scheme, where momenta and fields are evaluated at fractional and whole time steps (i.e. and with ) respectively. The complete set of equations then reads
[TABLE]
An advantage of using the leapfrog scheme is its accuracy up to second order in the time step and its symmetry under time reversal. Furthermore, the use of exponentiation in the evolution equation for guarantees that gauge links remain special unitary matrices throughout the numerical time evolution, also for finite . Probably the most important aspect is that the leapfrog equations are gauge covariant, i.e. they transform consistently under lattice gauge transformations. Gauge covariance also implies conservation of the discrete Gauss constraint given by
[TABLE]
In the expression for the Gauss constraint momenta are evaluated at while gauge fields and links are evaluated at . In this form the Gauss constraint is actually exactly conserved for finite time steps . It is a non-trivial task to find different update equations that are consistent with a discretized Gauss constraint, in particular if one starts from an action that is continuous in instead of discretizing already at the level of the action. This issue will be addressed in more detail in chapter 5.
2.3.3 Initial conditions on the lattice
The last missing piece in the numerical description of the boost invariant Glasma is to formulate the initial conditions (see eqs. 2.86, 2.87, 2.88 and 2.89) on the lattice in terms of transverse gauge links and the canonical momenta .
We have to find the pure gauge color field of a single nucleus in LC gauge. Starting in covariant gauge we discretize the charge density two-point function eq. 2.38 of the MV model in the transverse plane. The correlator then reads
[TABLE]
where with are transverse lattice spacings and are discretized transverse coordinates. The indices refer to the set of independent color sheets. The probability functional of the MV model is Gaussian, which makes it easy to generate random numbers on a computer according to the above two-point function: at each point in the transverse plane, each color component (of which there are ) and each sheet index , the color charge density is an independently sampled Gaussian random number with zero mean and variance . Since the MV model describes nuclei with infinite extent in the transverse plane and computers provide only a finite amount of memory, we have to use a finite transverse lattice with periodic boundary conditions. For simplicity, we use a square lattice with cells, lattice spacing and transverse length .
After generating a random color charge density on the lattice, one has to solve the discretized Poisson equation
[TABLE]
which is a lattice regularization of eq. 2.39. A fast and elegant way of numerically solving the Poisson equation above is to do so in momentum space with the help of the discrete Fourier transformation. We define the Fourier series
[TABLE]
where . The vector can be interpreted as the discretized momentum on the lattice. The above transformation can be readily calculated on a computer using the fast Fourier transformation (FFT). Using this definition the Poisson equation simply reads
[TABLE]
with the squared transverse lattice momentum
[TABLE]
The momentum representation also makes it possible to directly perform infrared and (optional) ultraviolet regulation. The regularized solution is then given by
[TABLE]
with the additional condition to eliminate the zero mode of , i.e. , in order to guarantee global color neutrality. In addition to the screening mass , we also introduce an ultraviolet regulator to regulate high momentum modes that are poorly resolved on the lattice. Thus, the discretized Poisson equation can be solved in the following way:
Generate a color charge configuration according to eq. 2.125 using a Gaussian random number generator. 2. 2.
Perform an FFT to obtain . 3. 3.
Solve the regularized Poisson equation according to eq. 2.130. 4. 4.
Perform an inverse FFT to obtain the solution on the transverse lattice .
The lightlike Wilson line is then given by (see eq. 2.40)
[TABLE]
where the subscript “A” denotes the right-moving nucleus. The transverse color field of nucleus “A” is
[TABLE]
which follows from performing the gauge transformation from covariant gauge to LC gauge on the lattice, see eq. 2.99. For the left-moving nucleus “B” we perform the analogous procedure and find
[TABLE]
The lattice formulation of eqs. (2.72) and (2.73) has been derived in [65]. The discretization of eq. 2.72 reads
[TABLE]
which is an implicit equation for . For there is an explicit solution given by
[TABLE]
The inverse can be readily computed, see section C.1.4. For eq. 2.134 needs to be solved numerically. The lattice version of eqs. (2.73) and (2.86), i.e. the discretized momentum , is given by
[TABLE]
where . The rest of the components are zero:
[TABLE]
To initialize the leapfrog integrator one has to define the gauge fields at , while the momenta are defined at :
[TABLE]
2.4 Glasma observables
In the following section we discuss a few standard results and phenomena using the numerical procedure outlined in the previous sections. We focus on the energy density and the pressure components of the Glasma, i.e. components of the energy momentum tensor .
Having already worked out the action eq. 2.80 and the canonical momenta eq. 2.82, eq. 2.83 it is a simple task to perform the Legendre transformation and obtain the Hamiltonian density
[TABLE]
which can be related to the energy density in the frame at mid-rapidity via , where the expectation value refers to averaging over the color charge densities of both nuclei. In order to identify the various terms in the energy density, we use the field strength components of the laboratory frame at mid-rapidity
[TABLE]
where and refer to the longitudinal field components. We also introduce the transverse components and analogously . Using these definitions we can split the energy density into four contributions:
[TABLE]
From the expression of the energy momentum tensor
[TABLE]
we can read off the other diagonal components: the longitudinal pressure and the transverse pressure given by
[TABLE]
Due to being traceless we have . The four energy density contributions, the energy density and the pressure components are the main observables that we are interested in.
In this section we consider head-on (i.e. central) collisions of gold nuclei at “RHIC-like” collision energies of roughly . Due to the limitations of the MV model and the ultrarelativistic approximation, the only energy dependence comes from the saturation momentum . However, the MV model by itself does not predict any dependence . A from-first-principles determination of collision energy dependence of the initial conditions, or rather dependence on the momentum rapidity of the colliding nuclei, is possible through the JIMWLK renormalization group equations. A proper discussion and application of the JIMWLK equations can be found in [18, 19, 67]. For the purposes of this thesis, we simply choose a value for by hand, such that it more or less corresponds to the case of Gold nuclei colliding at (see [68, 69]).
Ideally, fixing the saturation momentum would be done indirectly by choosing a value for the MV model parameter and the coupling constant , since naively one expects (see discussion below eq. 2.42). However, the number of color sheets and the screening mass , introduced as regularizations of the model, also play a considerable role in the relation between and [50]. Even though the MV model is comparatively simple, its parameter uncertainties are not insignificant. For the purposes of this thesis it will therefore be enough to choose parameters that are phenomenologically somewhat realistic, but it should be stressed that comparisons to experimental results should not be taken too seriously.
With this disclaimer in mind, we proceed as follows: at collision energies of one expects a saturation momentum (see e.g. [68, 69]). Since sets the relevant QCD energy scale, one can use the one-loop beta function to find a coupling constant or , which is the popular choice for in CGC literature. Assuming and ignoring all other dependencies, these considerations fix the MV model parameter to be roughly . Alternatively, one can use the phenomenological estimate eq. 2.31 and which approximately yields a similar value .
As we are interested in central collisions, the size of the circular overlap region of the two gold nuclei is , where is the nuclear radius with . Approximating the overlap region as a square of identical area, we set the transverse length of the simulation box to . Due to the use of periodic boundary conditions, boundary effects are completely neglected and the system size in the transverse directions is effectively infinite. The transverse size has to be large enough in order to accommodate the largest wavelengths in the color field of the nucleus. If color neutrality is realized at the size of the nucleus, the area of the transverse plane must be sufficiently large, which is the reason for choosing .
In order to regulate the infrared modes of the color fields of the single nuclei we use two different prescriptions: either (I) and global color neutrality by eliminating the zero mode of the charge density, or (II) which implements color neutrality at the size of nucleons. For prescription (II), the zero mode of the color charge density is also removed. Ultraviolet modes are cut off at , which should only insignificantly affect the results, but eliminates high momentum modes for which the lattice treatment of the system might be less than optimal.
For the number of color sheets we either choose , i.e. the single color sheet approximation, or up to in order to approach the continuum limit of the MV model.
The simulations are performed using as the color gauge group. We choose two instead of three colors to make use of a simple and elegant parametrization for (see section C.1) which makes the implementation of a real-time lattice Yang-Mills code comparatively easy. Furthermore, since there are only three color components in (compared to eight in ), performance and memory requirements of the code are much better, which will be particularly important in the three-dimensional setup. Most importantly however, the qualitative physics of the Glasma is essentially the same for any number of colors: phenomena such as the pressure anisotropy and the decay of Glasma flux tubes are similar for all groups . It is also possible to rescale results from simulations with dependent factors to give reasonable agreement with simulations of gauge fields [70]. Analytical calculations (see e.g. [71]) show that the initial energy density scales with the Casimir of the fundamental representation . This yields a rescaling factor of , which we will use to extrapolate all energy density contributions.
The code used to solve the lattice equations of motion and produce the results in this section was developed in Python and Cython [72] and is hosted on GitLab222The code is freely available at https://gitlab.com/dmueller/curraun_cy..
Figure 2.4 shows the energy density of a single collision event as a function of proper time and reveals the longitudinal flux tube structure of the Glasma, which, at least in the ultrarelativistic limit, is formed instantaneously at . Even though the MV model is homogeneous on average, single events sampled from the probability functional are highly random: the statistical nature of the color charge densities leads to random hotspots of energy. The expansion of flux tubes after , which generates transverse chromo-electric and -magnetic fields, leads to characteristic circular patterns in the transverse plane seen in the center plot. As the evolution continues, the energy density distributes across the transverse plane and the system becomes more and more dilute. At later times the four energy density contributions begin to converge towards the same value, as can be seen in fig. 2.5 on the left. Afterwards, the system shows free streaming behavior signaled by shown in fig. 2.5 on the right. The reason for this decrease of energy density is simple: as the colliding nuclei recede, the Glasma expands not only in the transverse plane, but also in the longitudinal direction and energy is transported away from the mid-rapidity region. This expansion occurs in a boost invariant manner and therefore the energy per unit rapidity, i.e. , becomes constant. Since the energy density and the field amplitudes decrease at later times, non-linear interactions become less important and eventually the system can be treated as effectively Abelian.
Let us now focus on some quantitative results for the Glasma energy density. In [69] it was shown that the initial energy density at of the MV model is not a particularly useful quantity to look at as it suffers from a UV divergence. However, this divergence dies down quickly for and at results become independent of UV regularization and the lattice spacing [54, 55, 69]. At we find an energy density of (including the color scaling factor ) for and for . Similar simulations [69] with as the gauge group and yield a value of . Performing the numerical evolution with and using the color scaling factors seems to underestimate the Glasma energy density. In [73] a different color scaling factor of was used, where is the dimension of the adjoint representation for . Using this other factor we find for and for . We see that this time the energy density is overestimated compared to . Obviously, playing such extrapolation games can never be exact, but when considering all other parameter uncertainties of the MV model, it is clear that still provides phenomenologically relevant “ballpark” results.
For fixed values of , and , the energy density also strongly depends on the number of color sheets. Plotting as a function of for various values of we find the results presented in fig. 2.6. In general, increases with , but this growth stops around with convergence being faster for larger values of . In the continuum limit one finds an approximate logarithmic dependency [54, 55], which is well reproduced in fig. 2.6. On the other hand, the single color sheet approximation shows more complicated behavior with a peculiar minimum around . In [49] it was shown that the case inadequately describes the generalized MV model of eq. 2.36 as the path ordering of the lightlike Wilson lines is completely neglected. In light of these results we can expect that the form of this curve depends on the number of colors as well and the minimum should probably not be taken too seriously.
The plot shows that for a given the energy density is significantly higher in the continuum limit compared to the naive single color sheet approximation. This is something to keep in mind and shows that for phenomenological applications one has to be careful when choosing . It is always possible to simply view as a phenomenological parameter that can be fitted to known results [49]. Nevertheless, it is still interesting to study the interplay of parameters and their consequences on observables.
The formation of purely longitudinal Glasma flux tubes has a peculiar feature: when computing the two pressure components and , one finds that, because transverse fields only build up later in the evolution, the longitudinal pressure is initially negative. In fact, due to we have and consequently the Glasma is maximally anisotropic directly after the collision with and . The time evolution of the pressure components can be seen in fig. 2.7. For the flux tubes expand and generate and , which decreases the anisotropy in the system until is slightly positive. In the purely boost-invariant Glasma one never encounters isotropization in the sense that or . After roughly (around the same time when ) the system settles into a free streaming expansion characterized by and . This happens because the expansion of flux tubes leads to equilibration of the four energy density components. It would be necessary to generate more transverse fields than possible through such an expansion to reach an isotropized state. It turns out that the free streaming limit is a generic feature of the boost invariant expansion of longitudinal flux tubes even in the Abelian limit: in [74] it was shown that the Abelian evolution of electric or magnetic flux tubes of size qualitatively agrees with fully non-Abelian Glasma simulations with free streaming setting in after .
The pressure anisotropy of the Glasma poses a serious problem for heavy-ion collision modeling, since a transition from the Glasma state into the QGP somewhere around is necessary. Famously, the QGP is well described by relativistic viscous hydrodynamics [13], which is able to support some pressure anisotropy , but not as large as predicted by the classical Yang-Mills evolution of the Glasma where . An ad-hoc solution to this problem is to simply assume some kind of mechanism that cures the anisotropy when switching from the Glasma to hydrodynamics. Even though this approach is successful (see e.g. [25]), it is unsatisfying from a theoretical viewpoint. A more recent and elaborate approach is to introduce an intermediate kinetic theory stage between the Glasma and the QGP [75, 76], which drives the system towards isotropization within . On the other hand, there exist possible mechanisms for reaching an isotropized system from the Glasma within classical Yang-Mills theory, which will be discussed in the next section.
Up until now we have only discussed observables derived from the energy-momentum tensor , but it is also possible to analyze the gluon spectrum of the Glasma by computing gluon occupation numbers. Although this quantity is not studied in this section, it is still important to mention the procedure, as it allows one to make more direct connections to experiments. We saw earlier that the expansion of the Glasma at later times is characterized by . As the field amplitudes die down, non-linear interactions are less relevant and consequently one can treat the system as an ensemble of weakly interacting gluons [68, 77]. Furthermore, using the assumption that the electric and magnetic contribution is equal at late times (at least on average), the Hamiltonian of the system (2.140) can be approximated by
[TABLE]
Performing a Fourier transformation one finds
[TABLE]
Identifying this with the Hamiltonian describing free particles
[TABLE]
where is the free dispersion relation and is the gluon occupation number in momentum space, one finds the relation
[TABLE]
However, due to the system being described by gauge theory, this identification is not unique. The residual gauge freedom, i.e. time-independent gauge transformations, has to be fixed. The standard practice [68, 77, 78, 79, 80] is to choose transverse Coulomb gauge
[TABLE]
which minimizes the functional . Computing the gluon occupation number , one can determine the total number of gluons per unit rapidity produced in a collision by integrating over all momenta. Assuming that all gluons decay into hadrons, it is possible to make direct comparisons to the number of charged particles per rapidity measured in heavy-ion collision experiments [24, 34, 68, 81]. This can be used to fix the MV model parameter directly using only classical Yang-Mills simulations of the Glasma. Other quantities such as the energy density at can not be measured directly, making the gluon spectrum a valuable tool for phenomenological applications.
2.5 Beyond the MV model
The goal of this chapter was to give a short presentation on the standard techniques and main phenomena of the boost invariant Glasma using the MV model as a simple approximation of nuclei. The MV model is highly restrictive due to its simplicity: we approximate the colliding nuclei to be infinitesimally thin in order to justify boost invariance and thus neglect any rapidity dependence of observables. We can only describe a crude approximation to central, head-on collisions, because the MV model is infinitely extended in the transverse plane. As such there is no notion of impact parameter dependence. Finally, the model only depends on one dimensionful parameter (apart from regularization parameters). Consequently, the transverse structure of the model is very simple: even though color charge densities are random on an event-by-event basis, the average density is completely homogeneous. A realistic model of a nucleus should also include its nucleonic and possibly even sub-nucleonic structure.
In this section we discuss a few extensions to this simple picture, which will naturally lead us to the main topic of this thesis.
2.5.1 Glasma instability and isotropization
As we saw in the previous section, the boost invariant expansion of the Glasma leads to a highly anisotropic system with vanishing longitudinal pressure. Even at very large proper times the system will remain in this state. This asymptotic behavior of the boost-invariant system stands in contrast to experimental indications of a nearly isotropized quark-gluon plasma, where after a few [10, 82, 83]. This scenario of fast isotropization is necessary for traditional relativistic viscous hydrodynamics to be applicable at early times.
Since boost invariance was put in by hand, it is a reasonable question to ask whether a system like this is actually stable against small variations which break this symmetry. This question led to a large number of publications (see e.g. [35, 36, 39, 41, 42, 84, 85, 86]) studying a similar setup: an initially boost invariant Glasma at is perturbed with rapidity dependent fluctuations and then evolved using 3+1 dimensional Yang-Mills equations in the frame. These studies found that the initially boost invariant Glasma is indeed unstable and that the pressure anisotropy becomes less severe at late times . Rapidity-dependent fluctuations (whose origin is thought to be quantum) therefore provide a possible mechanism to explain the transition from the pre-equilibrium Glasma to the equilibrated QGP. However, the general consensus is that this mechanism is not efficient enough to fully isotropize the system within the short time frame before the system becomes hydrodynamic at around . Without an intermediate stage described by kinetic theory [75, 76], this either means that the Glasma description (even including fluctuations) is insufficient or that the system might not be isotropic at early times. It should be noted that the scenario of fast isotropization ( – ) is not undisputed. The development of anisotropic relativistic hydrodynamics (aHydro) [87, 88] over the past few years has shown that the QGP can evolve hydrodynamically while still maintaining a large pressure anisotropy throughout its lifetime (see [89] and [90] for recent reviews and pedagogical introductions to this topic). In any case, it is clear that the boost invariant approximation neglects more than just simple rapidity dependence of observables and that there are interesting dynamics to be found in the 3+1 dimensional Glasma.
2.5.2 Realistic models of nuclei
In the beginning of this section we have outlined the shortcomings of the MV model and how it is only a very crude approximation of high energy nuclei. Fortunately, the MV model (see eqs. (2.35) and (2.36)) can be easily generalized by promoting the parameter to a function of the transverse coordinate via
[TABLE]
The color charge densities are still considered to be Gaussian, but now the variance of the fluctuation varies locally. This extension allows us to consider nuclei with finite radius by choosing to have e.g. the form of a Wood-Saxon distribution with appropriate parameters. In turn, finite nuclei allow us to describe off-central collisions. This is in fact how the earliest numerical studies of the Glasma with impact parameter dependence have been performed [68, 91].
The above model assumes nuclei to be homogeneous disks of diameter with possibly some smooth boundary, thus neglecting any structure inside the nucleus. This can be improved upon, by also considering the nucleonic structure of nuclei: since nuclei are comprised of nucleons it is reasonable to assume that the color charge density inside is not homogeneous (not even on average) and only becomes large in the vicinity of nucleons [91]. One can formulate even more accurate models for nuclei from its building blocks (i.e. nucleons) if one has an accurate description for single nucleons. The IP-Glasma [23, 24] is a state-of-the-art model for high energy nuclei that is constructed in the above described manner. In the IP-Glasma, the saturation momentum of a single proton is considered as a function of relative to the center of the proton, which can be related to the color charge density fluctuation . The question then arises how one should choose in order to arrive at an accurate representation of a full nucleus. This leads to the most important feature of the IP-Glasma model: the exact form of is extracted from experimental deep inelastic scattering (DIS) data [23, 24]. The color charge density of neutrons is determined in the same manner. Using this procedure the IP-Glasma has become the most successful model when applying CGC/Glasma to phenomenology. For instance, it successfully describes the charged particle multiplicity [33] computed from the gluon occupation numbers at both RHIC and LHC not just for the nucleus-nucleus case but also for proton-nucleus and proton-proton collisions. When coupled to relativistic viscous hydrodynamics after the Yang-Mills evolution, it correctly describes anisotropic hydrodynamic flow at RHIC and LHC [25].
2.5.3 Rapidity dependence and finite longitudinal extent of nuclei
Even though the Glasma initial conditions were derived under the assumption of boost invariance, it is still possible to simulate collisions with rapidity dependence, at least in an approximate manner. The main idea, as outlined in [32], is to introduce a momentum rapidity dependent MV model parameter for both nuclei “A” and “B” via
[TABLE]
with a parameter (see also [92] for the rapidity dependence of ). This relation can be justified by recalling that the saturation momentum depends on the longitudinal momentum fraction of the gluons and that the relevant can be related to the momentum rapidity of the colliding nuclei [17]. The IP-Glasma uses a very similar procedure to obtain rapidity dependent initial conditions [33].
The rapidity is then viewed as an external parameter and one performs independent, boost invariant simulations for a range of values of . The computed observables for each , such as the energy density or gluon numbers, are identified with the observables at that rapidity.
The same simulation procedure is used in [34], but the rapidity dependent initial conditions are not obtained as outlined above. Instead, the initial conditions are computed from a JIMWLK evolution, which takes into account quantum effects that lead to the saturation momentum growing with rapidity.
Treating each rapidity direction independently obviously neglects any interactions between fields at different rapidities, which is a reasonable approximation at later times in the evolution, but can not be completely accurate at the earliest stages of the Glasma. Furthermore, the procedure relies on using the boost invariant initial conditions eqs. (2.72) and (2.73), as there are no analytical results for rapidity-dependent Glasma initial conditions. One could interpret this procedure as using a different effective (boost invariant) description of the nucleus for each rapidity slice and simply “stitching” together the results for each value of . On the other hand, if one assumes that there should be a single classical field configuration valid for a single nucleus associated with rapidity dependent fields, then it is clear that such a field configuration can not have a -shape as assumed in the ultrarelativistic limit. This leads us to the main goal of this thesis. In the next chapters we develop a numerical method to solve the classical Yang-Mills equations for collisions of nuclei with finite longitudinal extent. This method forces us to work in the frame instead of using comoving coordinates in order to correctly incorporate initial conditions for nuclei with finite thickness along or . By virtue of breaking boost invariance explicitly, one finds rapidity dependent observables such as the rapidity profile of the energy density of the Glasma (see section 4.4). This rapidity dependence is a consequence of the finite longitudinal extent and the longitudinal structure of the color fields of the nucleus before the collision. It should be interpreted as a purely classical effect as there is a priori no notion of growing with rapidity. Since the problem will be formulated in a different frame and previously important assumptions like boost invariance will be loosened, we will have to construct the method from the ground up. This will be the main task in the next chapter.
Chapter 3 The 3+1D Glasma on the lattice
In this chapter we develop the numerical method we use to simulate heavy ion collisions in the CGC framework without relying on the ultrarelativistic limit and the symmetry under longitudinal boosts. In section 3.1 we argue that the laboratory (or center of mass) frame is a more suitable choice of frame when considering nuclei with finite extent along the longitudinal axis. We also mention some assumptions that render the numerical problem more tractable. In section 3.2 we introduce the colored particle-in-cell (CPIC) method that we use to solve the problem numerically on a lattice and derive the discretized equations of motion for the color fields and color currents used in a simulation. As the numerical description differs from the boost invariant case, the initial conditions for single nuclei on the lattice have to be formulated in laboratory frame temporal gauge (see section 3.3). In section 3.4 we explain how we compute the energy-momentum tensor on the lattice, which will be the most important observable that we study. Finally, in section 3.5, we discuss some details regarding the actual implementation of the code such as optimizations, some subtleties regarding the initial and boundary conditions, and how to parallelize the simulation.
3.1 Collisions of nuclei with finite longitudinal extent
Let us quickly recapitulate how to derive the color fields of a single nucleus (for details see section 2.1). Starting from a color charge density moving along the axis and fixing the covariant gauge , the only non-vanishing component of the color field is found by solving the Poisson equation
[TABLE]
The subscripts “A” and “B” are used to differentiate the two nuclei. This solution is exact for arbitrary color charges if the charge densities only depend on one of the light cone coordinates . In the last chapter we considered the longitudinal support along the longitudinal axes finite, but only in order to account for path ordering along when switching to LC gauge. Here, we would like to associate finite longitudinal support with the actual physical Lorentz contracted extent of relativistic nuclei proportional to , where is the nuclear radius and is the Lorentz factor. Therefore, at no step in the calculation can we take this longitudinal extent to be infinitesimal as we did in the previous chapter. This immediately poses a problem when considering collisions of two such color fields. A schematic picture of this is shown in fig. 3.1. The charge densities (which are basically the only input for the calculation) are not invariant under longitudinal boosts. The Glasma fields produced in a collision will depend on space-time rapidity and the system cannot be reduced to 2+1 dimensions. With finite longitudinal support the collision and creation of the Glasma is not instantaneous, although the formation still happens very fast with the formation time being proportional to . Without employing any approximations (such as expanding in ) one has to solve the fully non-perturbative Yang-Mills equations in order to compute the rapidity-dependent color field after the collision. Furthermore, without a point-like origin, the Minkowski diagram cannot be separated into distinct, causally independent regions in an unambiguous way. The separation ansatz, eqs. (2.69) and (2.70), is not appropriate in this more general case. In conclusion, it seems that many assumptions that apply to the boost invariant Glasma simply cannot be used for the collisions we are considering and a different approach is necessary.
This brings us to the main novel feature of our 3+1 dimensional setup. Since the separation of space-time into distinct regions is ambiguous and a future light cone cannot be identified, we choose a “brute force” approach and simply stay in the center of mass frame. We therefore use the laboratory time and the beam axis coordinate as our coordinate set. This means that we also have to consider the dynamics of the color current. Recall that in the boost invariant case, the color charges of the nuclei are only relevant for finding the initial Glasma fields. The infinitesimally thin color currents only form the boundary of the future light cone and do not affect the time evolution of the Glasma after its formation. On the other hand, with finite longitudinal extent the color currents, in addition to the color fields of the nuclei, have to be included in the description and treated as dynamic degrees of freedom.
There is one restriction we can put on the color currents: we can assume them to be eikonal, i.e. the color charges do not recoil or lose any momentum. Thus, the color charges are only affected by parallel transport, meaning that they can rotate in color space as they undergo the collision, but their lightlike trajectories remain fixed. This color rotation is a consequence of non-Abelian charge conservation via the continuity equation (see eq. 2.7). Writing this conservation law down for the nucleus “A” we find
[TABLE]
In the asymptotic past of the nucleus the density is indeed independent of , but as soon as the two nuclei collide the color current is parallel transported due to the color field which is comprised of the color field of the second nucleus “B” and the “” component of the Glasma fields. Since the color charges are only rotated, we can use the ansatz
[TABLE]
where is the unaffected current at . The unitary matrix fulfills the relation
[TABLE]
with the general solution given by the Wilson line
[TABLE]
An analogous solution can be found for the current of the left-moving nucleus “B”. In the previous chapter the currents are unaffected due to the choice of Fock-Schwinger gauge which eliminates the ( resp.) component at the ( resp.) boundary of the future light cone, but in collisions with finite extent we cannot exploit gauge freedom to set both light cone components and to zero simultaneously everywhere. Parallel transport therefore needs to be accounted for in the 3+1 dimensional setting.
Now that we have understood the scope of the problem at hand, let us formulate it in detail. First we have to pick an appropriate gauge. The obvious choice for numerical time evolution is the laboratory temporal gauge . Transforming from the covariant gauge to the temporal gauge we have to solve
[TABLE]
where the subscript “cov” is used to explicitly denote covariant gauge fields. The general solution of this equation is given by the time-like Wilson line
[TABLE]
where the spatial coordinate vector is given by . If we apply temporal gauge to the color field of nucleus “A”, the Wilson line reads
[TABLE]
where we used and . Using the fact that only depends on and via we can reparametrize the integral and find that is in fact the same transformation matrix associated with switching from covariant gauge to LC gauge eq. 2.20. Writing the path-ordered exponential first as a product integral
[TABLE]
and performing the coordinate transformation from to , while keeping fixed, yields
[TABLE]
This fact can also be intuitively understood in a diagram, see fig. 3.2. Using an analogous reparametrization one finds that can be identified with the spatial Wilson line starting at :
[TABLE]
where . As a consequence, it is easy to show that the longitudinal component vanishes in temporal gauge. Since the spatial Wilson line fulfills
[TABLE]
we find
[TABLE]
The only remaining components are the transverse gauge fields given by
[TABLE]
Analogous arguments can be made for the left-moving nucleus. The color fields of single nuclei in temporal gauge are therefore completely equivalent to the LC gauge solutions.
Before the collision, at some initial time , when the nuclei are sufficiently well separated, it is possible to superimpose the transverse gauge fields
[TABLE]
The total color current is also a superposition
[TABLE]
where only the longitudinal component is non-zero:
[TABLE]
This initial condition is then evolved to by solving the Yang-Mills equations (see eq. 2.3), which read
[TABLE]
where the color currents are given by
[TABLE]
Parallel transport is accomplished by
[TABLE]
The parallel transport of the color currents (and the signs used in the exponential) can be better understood if one thinks of the color charge density of the nuclei to be comprised of an ensemble of point-like particles. Consider a single point-like particle in 1+1 in a background color field in temporal gauge. Assuming that it follows a lightlike trajectory defined by (such that ), the color charge density and longitudinal current is given by
[TABLE]
Due to non-Abelian charge conservation the particle’s color charge is time dependent. Requiring that the continuity equation (see eq. 2.7)
[TABLE]
is fulfilled, yields
[TABLE]
We can generalize this to a smooth distribution of color charge by letting depend on . The color charge density is then given by
[TABLE]
with
[TABLE]
which is exactly eq. 3.20. Repeating the same steps for a particle moving in the opposite direction defined by we find
[TABLE]
which corresponds to eq. 3.21. The negative sign in the exponential can be traced back to the opposite sign occurring in the spatial color current, i.e. .
We see that the equations for the 3+1 dimensional setup are quite complicated. Not only do we need to include the additional longitudinal axis , but due to parallel transport the color fields couple to the color currents in a highly non-linear way. Fortunately, there is a numerical method that can handle dynamic color charges in a gauge covariant manner. This will be the topic of the following sections.
3.2 The colored particle-in-cell method
The colored particle-in-cell (CPIC) method is a numerical approach to simulating an ensemble of colored point-like particles coupled to non-Abelian gauge fields. It combines the real-time lattice description of gauge fields with classical particles obeying Wong’s equations [93], which is the non-Abelian equivalent of the Lorentz force. To name a few applications, CPIC has been successfully used in hard thermal loop simulations [43, 94], where colored particles act as hard (i.e. ultraviolet) degrees of freedom which interact with soft (infrared) classical color fields. It can also be used to study non-Abelian plasma instabilities [44, 95] in the QGP. When including collisions among the colored particles, it is possible to simulate energy loss and momentum broadening of jets in the QGP [96].
The CPIC method originates from the popular particle-in-cell (PIC) method used in (Abelian) plasma physics [97]. Pedagogical introductions can be found in [98, 99, 100]. PIC was developed to accurately simulate electromagnetic plasmas with typically very large numbers of charged particles coupled to the electromagnetic field. In a PIC simulation fields are usually sampled in cells on a numerical grid and a finite difference approximation is used to numerically solve Maxwell’s equations. On top of this lattice, an ensemble of charged particles is able to move through the simulation volume. The particle trajectories are only affected by the fields within their closest (and possibly also neighboring) cells. This approach solves a number of problems encountered in the simulation of plasmas: instead of having to compute the Lorentz force of each particle pair in the plasma individually, which is an operation that scales quadratically with the number of particles , one only has to consider the field in the immediate vicinity of each particle (in the simplest approximations just a single cell), which only involves operations. The PIC method therefore makes it possible to perform large scale simulations of systems with high numbers of particles (usually ). Furthermore, compared to early electrostatic simulations, it allows for a fully dynamic treatment of the electromagnetic field. Consequently, the fields exhibit a finite speed of propagation and can be consistently coupled to relativistic particles, enabling the simulation of fully relativistic systems such as the acceleration of electrons in laser wakefield experiments [101].
An important aspect of PIC is how to relate the field degrees of freedom to the particle degrees of freedom. Fields are defined at lattice sites, whereas particles have continuous positions and velocities and move freely across the lattice. This is accomplished through a procedure called interpolation or weighting. Interpolation is necessary in two directions: interpolating from the particles to the grid in order to compute charge densities and currents on the lattice and interpolation from the grid to particles to obtain the force acting on them. Each particle is given a certain shape characterized by a continuous shape factor or weight . Given the continuous position of a particle, the shape factor determines how its point-like charge is distributed among the cells on the discrete lattice and, conversely, how the fields act on the particle.
The concept of interpolation is best explained by giving a specific example: consider a simple 1+1 dimensional setting with a single particle with charge following a trajectory . How should one relate this charge to the discrete charge density defined at the lattice points , where is the lattice spacing? The solution is to specify the shape of the particle via the shape factor and perform the mapping
[TABLE]
Depending on the exact form of the charge is distributed among multiple neighboring cells or just a single cell. Charge conservation is guaranteed by requiring that the total charge on the lattice is always the same at any time independent of the particle position, i.e.
[TABLE]
This leads to the normalization condition for charge conserving shape factors
[TABLE]
where is an arbitrary position and the sum runs over all lattice sites . In principle, shape factors can be arbitrary functions as long as the above constraint holds, but for simplicity they are usually constructed from piecewise polynomials and categorized by the order of the polynomials. In fig. 3.3 two simple shape factors are presented: the nearest-grid-point (NGP) or zero-order weight, and the cloud-in-cell (CIC) or first-order weight. The NGP scheme is defined by
[TABLE]
while the CIC scheme is given by
[TABLE]
In this one-dimensional example the NGP scheme interpolates the charge to the (as the name implies) nearest neighboring lattice site, whereas the CIC scheme distributes the charge among two neighboring cells (except if the particle is exactly at a lattice site ). Higher order shape factors can be systematically constructed using quadratic, cubic or higher order polynomials [102].
The interpolation step can easily be generalized to multiple particles by summing over each individual particle
[TABLE]
In more than one spatial dimension the shape factor is defined as the product of the one-dimensional shape factors, e.g. for three spatial dimensions and a trajectory we have
[TABLE]
where are the coordinates of lattice site . Shape factors of particles are also relevant when performing the inverse interpolation step, i.e. interpolation from fields to particles. In this simple 1+1 dimensional example using the temporal gauge condition, there is only one gauge field component defined at the lattice sites . The local gauge field acting on a particle at position is computed using the mapping
[TABLE]
The use of the same shape factors for both interpolation directions is crucial, as this eliminates numerical instabilities caused by an unphysical self-force, i.e. a spurious force due to particles generating fields that act on themselves (see [103] and section 5.3.3 of [99]). This elimination of numerical instabilities only works in the electrostatic case, where one assumes that particle velocities are much smaller than the propagation speed of electric and magnetic fields, i.e. the speed of light. In the more general electromagnetic case, new numerical instabilities due to self-interaction can occur. A prominent example is the numerical Cherenkov instability, which will be discussed in chapter 5.
As particles move across the lattice, they generate a charge current that affects the time evolution of the fields. In order to preserve the Gauss constraint, the current must be locally charge conserving and satisfy the continuity equation. In the 1+1 dimensional example, the continuity equation reads
[TABLE]
A straightforward discretization of the above equation is given by
[TABLE]
Here refers to the time step and to the charge current density on the lattice. The current is evaluated at shifted lattice sites and at fractional time steps. This assignment is known as the Yee scheme [104] which is an example of a staggered grid. The main benefit of using the Yee scheme is that simple forward differences can be interpreted as central differences, which have quadratic accuracy:
[TABLE]
We now consider a single particle with charge and trajectory in the NGP scheme. The charge current generated by the moving particle can be computed directly from the discretized continuity equation with a few simple assumptions. We consider a particle moving in the positive direction. Notice that if the particle moves in such a way that its nearest neighbor cell does not change from one time step to the next, then the discretized charge density stays the same. In this case no current is generated. However, once the particle crosses from one cell into the next the time derivative of is non-zero and a non-vanishing current must be generated. Assuming that the NGP of is and the NGP of is we have and . Plugging this into eq. 3.39 evaluated at and we obtain
[TABLE]
Using the boundary condition that no spurious current should be generated far away from the particle position, it is intuitive to induce the current at the mid-point . The charge current that satisfies the continuity equation is therefore
[TABLE]
At all other lattice sites, the current is set to zero. Reversing the above example and considering a particle moving from to yields the same current with opposite sign
[TABLE]
The NGP scheme is peculiar because current is only generated when a boundary crossing happens. In the case of a single particle, this leads to sudden and high spikes in the discrete current density, which introduce unwanted numerical noise when coupled to the fields. This artifact can be avoided with the CIC and higher order schemes which smear the particle charge among multiple cells and induce smoother currents. In section 3.5 we will see how this problem can be mitigated even for the NGP scheme.
The focus of this section is how interpolation is performed in PIC simulations. Although we left out details regarding the computation of forces and the time evolution of fields, a summary about PIC is incomplete without mentioning the typical simulation cycle. A PIC simulation proceeds as follows:
Interpolation from particles to fields: given all particle positions, velocities and the shape factor, determine the lattice charge and current density. 2. 2.
Field evolution: solve Maxwell’s equations on the lattice using the current density and previous fields. 3. 3.
Interpolation from fields to particles: determine the local field acting on each particle with the same shape factor as in the first step. 4. 4.
Particle evolution: use interpolated fields to compute the force acting on each particle and update velocities and positions accordingly.
This cycle is repeated for each time step until the simulation is complete. At the end of each cycle one can compute various observables such the energy momentum tensor of the fields or particle densities.
The CPIC method generalizes the concept of dynamic point-like particles in cells coupled to discretized fields to non-Abelian gauge theory. The main ingredient in PIC, namely the interpolation step, can be readily generalized to reflect non-Abelian gauge covariance. Color currents induced by the movement of colored particles can be computed from the non-Abelian continuity equation. We apply CPIC to the problem of colliding nuclei in the sense that we approximate the smooth color charge densities of the nuclei as an ensemble of colored particles. A schematic overview of this approach is shown in fig. 3.4. For this purpose, we only use CPIC in a very limited way, because we assume the colliding nuclei to be recoilless. We can therefore neglect any deceleration of color charges and only care about parallel transport. This greatly simplifies the problem since particle movement can be treated one-dimensionally. In the following subsections we derive the equations for gauge fields on the lattice and colored particles moving along lightlike trajectories.
3.2.1 Field equations
The starting point of our derivation is the Wilson action section 2.3.2 coupled to an external color current (see section A.3 for some clarifications regarding the shorthand lattice notation)
[TABLE]
where and are the color charge and current density. The current is identified with the spatial components of the four-vector . This action is slightly inconsistent in the sense that gauge links are used to formulate the pure Yang-Mills part, while the interaction term with the external current is written in terms of the gauge fields . Even worse, due to the presence of the term, the action is not even strictly gauge invariant under arbitrary gauge transformations. This is not an artifact of the lattice description, but is also true in the continuum case [105]. Nevertheless, the action still yields gauge covariant equations of motion (EOM) and constraints upon variation. The constrained variation of gauge links is described in section C.2. The variation of the coupling term is simply given by
[TABLE]
which is invariant under and . Performing the variation of the Yang-Mills part (see section C.3) with respect to temporal components and yields the Gauss constraint
[TABLE]
while variation with respect to spatial components and yields the EOM
[TABLE]
Here, denotes the anti-hermitian, traceless part, see eq. A.20. Equivalently, we can write the Gauss constraint as
[TABLE]
and the EOM as
[TABLE]
where we have used
[TABLE]
and the identity (see eq. A.24)
[TABLE]
Both the constraint and EOM are gauge covariant under lattice gauge transformations. In order to perform time evolution we use temporal gauge fixing
[TABLE]
where is the regular hypercubic lattice (see eq. A.26). The temporal plaquettes then reduce to the product of consecutive spatial gauge links (see eq. 2.101)
[TABLE]
The goal is now to express all gauge links in the next time slice in terms of previous gauge links. Once the temporal plaquettes are known from the EOM, the next gauge links can be computed via
[TABLE]
At a first glance, the EOM can only be solved for the anti-hermitian traceless part of the temporal plaquette , but since elements of in the fundamental representation are completely specified by real parameters, there are enough equations to reconstruct the full matrix . However, obtaining from depends on the parametrization and on the specific gauge group. For we use the parametrization (see section C.1)
[TABLE]
where , are the Pauli matrices. The real-valued parameters , fulfill the constraint (see eq. C.8)
[TABLE]
Using
[TABLE]
the EOM can be written as
[TABLE]
To obtain , we take the positive branch of eq. 3.58
[TABLE]
assuming that the changes from one time slice to the next are “small”, i.e. the temporal plaquette will be “closer” to than to . Obviously, if exceeds , the time evolution is ill-defined. This puts a limit on the field amplitudes and lattice spacings used in the simulation. The equations can be cast into a more familiar form by introducing the chromo-electric field via
[TABLE]
The lattice EOM then read
[TABLE]
and the temporal plaquette in is (see eq. C.30)
[TABLE]
which is used to update gauge links via eq. 3.56. The Gauss constraint reads
[TABLE]
or simply
[TABLE]
where is the backward gauge-covariant finite difference.
At this point a remark about the accuracy of these discrete EOM is necessary. Recall that the Wilson action is accurate up to quadratic terms in the lattice spacing and the EOM derived from this action should exhibit the same accuracy. At first glance the EOM eq. 3.63 seem to resemble a forward Euler scheme, due to the electric field and the gauge links “living” in the same time-slice according how they transform under lattice gauge transformations. Forward Euler schemes are only linearly accurate in the time step , which seems to contradict the accuracy of the action these equations were derived from. This apparent contradiction is resolved by interpreting the electric field to be evaluated in-between two time steps. The reasoning for this is as follows: the plaquette best approximates the field strength tensor at the center of the plaquette . Electric fields are defined via eq. 3.62, i.e. they are related to temporal plaquettes , which are naturally centered in-between time slices. The field therefore corresponds to (no sum over is implied) in the continuum. The EOM eq. 3.63 are therefore actually equivalent to a leapfrog scheme, which is quadratically accurate in the time step. We refer to this method as the Lagrangian leapfrog approach from now on.
Although the above presented time evolution follows directly from varying the Wilson action, it is less commonly used [78, 79, 106]. The more popular method is the Hamiltonian leapfrog approach. It based on the Kogut-Susskind Hamiltonian [107], which can be derived from the Wilson action by taking the continuous time limit and performing a Legendre transformation. Since we have already derived the discrete equations eq. 3.49, it easy to take the limit using
[TABLE]
In this limit the lattice reduces to the spatial three-dimensional lattice , where denotes the continuous time coordinate. The EOM in temporal gauge then read
[TABLE]
The quadratic accuracy arises naturally from the time reflection symmetry of the Wilson action. Identifying the time derivative of spatial gauge links with the chromo-electric field
[TABLE]
we find
[TABLE]
and therefore
[TABLE]
In the Hamiltonian approach the chromo-electric field is identified with the conjugate momentum of the gauge links. Repeating similar steps, the Gauss constraint reads
[TABLE]
or simply
[TABLE]
In order to solve the first-order differential equations eqs. (3.69) and (3.71) numerically, the time coordinates is (again) discretized and one obtains the leapfrog equations
[TABLE]
which are accurate up to second order in like the Lagrangian leapfrog approach. The momenta (electric fields) are shifted by half a time step compared to the gauge links, which is again reminiscent of a leapfrog scheme. Note that the use of the exponential in the gauge link update guarantees that gauge links remain in .
Both methods agree up to quadratic order in . The only difference lies in the update of gauge links from one time slice to the next and the definition of the electric field. To see this we first note that in the Hamiltonian leapfrog method, the temporal plaquette is given by
[TABLE]
For the exponential can be computed via
[TABLE]
where and . Expanding up to third order in we find
[TABLE]
In the Lagrangian leapfrog method, the chromo-electric field is identified as (see eq. 3.62)
[TABLE]
which can be thought of as naturally centered around the mid-point between and . The temporal gauge EOM look identical to the Hamiltonian leapfrog method. The temporal plaquette in is given by
[TABLE]
where . With the above definition of the electric field we then find
[TABLE]
The Lagrangian and Hamiltonian leapfrog methods agree up to a cubic term in and are interchangeable since both are only correct up to quadratic order. The Hamiltonian leapfrog equations were used in the first two publications [1, 2], but during the development of the numerical improvements introduced in [4] and presented in chapter 5 it became apparent that the Lagrangian leapfrog approach is necessary for more complicated actions than the standard Wilson action.
3.2.2 Particle equations
In this section we derive a locally charge conserving description of colored particles coupled to lattice Yang-Mills fields. As already outlined earlier, the trajectories of particles are completely fixed and we only need to derive how color charges are parallel transported as they move across the lattice.
Local charge conservation is closely related to the Gauss constraint and gauge invariance. Without external color sources, the lattice EOM eq. 3.49 exactly preserve the Gauss constraint eq. 3.48 (see section C.4) due to lattice gauge invariance. If external currents are included, the constraint is only conserved if and only if the non-Abelian continuity equation holds
[TABLE]
As colored particles move across the lattice, their charge density and color current must satisfy the above equation. Similar to the interpolation step derived previously in the Abelian case, we can use the non-Abelian continuity equation to derive the charge conserving current given a certain interpolation scheme.
Our setup is effectively 1+1 dimensional in terms of particle movement due to the eikonal approximation. It is therefore enough to consider a single point-like particle with color charge moving along the longitudinal axis at the speed of light. In general, is a function of time because of parallel transport and the color current only has a non-vanishing longitudinal component .
We choose the NGP scheme for the interpolation of charges to the lattice due to its simplicity. If the particle moves in such a way that its NGP does not change, there is no current generated, even if the actual particle position varies. This situation would be slightly more complicated if we had not fixed temporal gauge. Then, even if the particle is at rest, its color charge would be parallel transported along the time axis such that the charge would be gauge-covariantly constant
[TABLE]
In our setup, temporal gauge is used to solve the EOM, so we can ignore this complication of the more general case. It is then only necessary to study boundary crossings of particles as in the Abelian case. For this, we consider a particle moving in the positive direction. Assume that at time its NGP is and in the next time slice its NGP is . The charge density on the lattice then reads
[TABLE]
where is the volume of a cell. The charge density is zero everywhere else. We have also accounted for a possible color rotation of the charge. Plugging this into the continuity equation evaluated at and we find
[TABLE]
Analogous to the PIC derivation, we assume that only is non-zero. We then find from the first equation
[TABLE]
Using this in the second equation fixes the parallel transport of the particle charge
[TABLE]
We can repeat the same steps for a particle moving in the opposite direction. Assuming that at its NGP is and in the next time slice the NGP is , the charge density is given by
[TABLE]
The continuity equation yields
[TABLE]
In this case, the color current and the color charge are given by
[TABLE]
and
[TABLE]
Note that the slightly unintuitive signs of the color currents (i.e. a negative sign for the right-moving particle and vice versa) arise because are the spatial components of with lowered indices. Raising the indices yields and one obtains the expected signs.
In principle, one could also consider non-Abelian generalizations of higher order interpolation schemes such as the first-order CIC scheme. The CIC scheme for CPIC simulations has been derived in [95], but it has a slight flaw: the trace of the color charge squared of each particle should be a conserved quantity, because the charge is only affected by rotations. In the CIC scheme this conservation is violated and only restored in the continuum limit (see section C of [95]). On the other hand, the NGP scheme conserves this quantity exactly, i.e. for each individual particle we have
[TABLE]
The main disadvantage of the NGP scheme is that current is only generated at each boundary crossing. This problem can be circumvented by simply using a high number of particles per cell, such that at each time step the color current is non-zero.
3.2.3 Simulation cycle
Having derived all the necessary equations, we can summarize what a single simulation step entails in the Lagrangian leapfrog approach. The Hamiltonian leapfrog approach is completely analogous. Assuming that the ensemble of particles, and are known initially, the simulation cycle reads:
Interpolation from particles to fields: for each particle, determine the NGP from their current position and compute the charge density on the lattice. Assuming that the NGP is , the charge density receives the contribution
[TABLE]
where .
Optionally, check if the Gauss constraint is satisfied
[TABLE]
to make sure that the interpolation procedure works correctly.
Determine if a boundary crossing happens between the current and the next time slice and compute the longitudinal current density on the lattice. For right-moving particles the contribution is
[TABLE]
and for left-moving particles it is
[TABLE] 2. 2.
Update fields: Given , and from the previous step, compute via the EOM eq. 3.63
[TABLE]
Compute from and via eq. 3.64
[TABLE]
and eq. 3.56
[TABLE] 3. 3.
Interpolation from fields to particles: update the color charge of each particle using the gauge links . Right-moving particles are updated via eq. 3.89
[TABLE]
whereas left-moving particle use eq. 3.95
[TABLE]
given that is their NGP at time . 4. 4.
Particle evolution: update the position of each right-moving particle via
[TABLE]
and for each left-moving particle
[TABLE]
where is the unit vector in the longitudinal direction.
This form of a typical simulation cycle closely resembles traditional PIC simulations, but in practice some steps can be simplified. For example, the update to color charges and particle positions can be immediately applied after computing charge densities and currents, thus avoiding multiple loops over each particle. Similarly, the update of the electric fields and the gauge links at a given lattice site can be combined as well.
3.3 Initial conditions on the lattice
In this section we discuss how to put the temporal gauge initial conditions derived in section 3.1 on the lattice. Consider a right-moving nucleus described by the generalized MV model (see eqs. (2.35) and (2.36))
[TABLE]
At some initial time we can express this correlator in laboratory frame coordinates.
[TABLE]
The additional factor of is due to the scaling property of the Dirac distribution . On the three-dimensional lattice this correlator can be readily discretized via
[TABLE]
where is the longitudinal lattice spacing and is the transverse spacing. In practice such a random configuration can be easily generated, because for the charge density is completely uncorrelated. It is therefore sufficient to generate individual Gaussian random numbers for each color component with dependent variance and zero mean at each point on the lattice.
3.3.1 Sampling charge density with particles
In principle, it is a simple task to generate ensembles of particles that mimic a given charge density on the lattice, but there are some subtleties that are discussed in the following.
Starting with a simple example, one could, for each lattice site , generate a single particle with charge , where is the cell volume. This would exactly replicate the original charge density when performing the NGP interpolation step. However, in the derivation of the current in the NGP interpolation scheme, we saw that particles only generate current when they cross cell boundaries such that their NGP changes from one time step to the next. If the time step is , this only happens every second time step. While correct in principle, in the sense that the Gauss constraint is satisfied and the EOM are solved correctly, the time evolution will be adversely affected by numerical noise introduced by the current.
A simple solution to this problem is to use more than one particle per cell such that at each time step at least one cell boundary crossing per cell occurs. This leads to smoother currents and suppresses noise. Assuming that the time step is , where is an integer with , one uses particles per cell and places them equidistantly within each cell. To circumvent rounding errors in the computation of the NGP, we avoid placing particles exactly at the boundaries of cells. The charge within a cell can be distributed equally among the particles with .
This procedure only partially solves the problem of sudden changes in the current. Due to the equal distribution of charge within a cell, the same particle charges cross the boundary for consecutive time steps. As a result, the current generated by the interpolation step stays the same until all particles per cell have moved into the next cell. When all particles have moved, the current jumps a full lattice spacing . This sudden change also negatively impacts the time evolution of the fields.
The smoothness of the current can be further improved because the distribution of color charge among multiple particles per cell is underdetermined. Consider all particle (color) charges for a fixed transverse position lined up along the longitudinal axis. All charges with index belong to the same NGP cell with index . Within each cell we can modify the particle charges via
[TABLE]
without changing the total charge in the cell. We iteratively apply such small changes to smoothen the sub-lattice distribution of charges among particles. For example, one can pick the new first order finite difference of modified charges to be the average of adjacent derivatives via
[TABLE]
The above relation can be solved for :
[TABLE]
In order not to affect the overall charge in an NGP cell, one has to leave out the charges at the right boundary of each cell, i.e. . The above procedure eventually leads to convergence by applying it to each particle (except at boundaries) in each cell multiple times. The result of the above iteration is that the second order finite difference becomes constant in each cell. Since is directly proportional to the third order finite difference, we can write
[TABLE]
Requiring that the procedure converges, we have and thus
[TABLE]
which shows that the second order difference is constant. We therefore refer to the above iteration procedure as second-order charge refinement.
Similarly, we can derive a fourth-order refinement procedure by smoothing out the third order differences. This leads to
[TABLE]
which yields
[TABLE]
Again, since is proportional to the fifth-order finite difference, we can show that this leads to constant fourth-order finite differences. The above procedures work for both Abelian and non-Abelian charges: for color charges one performs the refinement steps for each color component independently.
In practice, we find the fastest convergence if we first apply the second-order refinement and afterwards the fourth-order refinement. The result of using these refinement procedures is shown in fig. 3.5.
3.3.2 Initial fields
The gauge field in covariant gauge is obtained through solving the two-dimensional Poisson equation. In each transverse plane we solve
[TABLE]
This can be done completely analogous to the boost invariant case using the discrete Fourier transformation, see section 2.3.3. In momentum representation the equation reads
[TABLE]
with the squared transverse lattice momentum
[TABLE]
Implementing infrared and ultraviolet regularization, the solution is given by
[TABLE]
where is the infrared regulator (screening mass) and is the ultraviolet regulator. After performing the inverse discrete Fourier transformation, one is left with the covariant gauge field at each lattice point .
In section 3.1 we showed that the solution to the Yang-Mills equations in temporal gauge involves a temporal Wilson line. Through reparametrization of the path-ordered exponential, it is possible to show that the temporal Wilson line is equivalent to the spatial Wilson line along . The longitudinal component in covariant gauge is simply . The spatial Wilson line is computed by constructing the path-ordered exponential starting at and ending at some finite .
[TABLE]
On the lattice, the Wilson line can be computed layer by layer, starting at some large far away from the nucleus where . We then have
[TABLE]
with the initial condition
[TABLE]
The transverse gauge links in temporal gauge are then given by
[TABLE]
which follows from the lattice gauge transformations of gauge links, see eq. 2.99. Longitudinal links are set to unit matrices in accordance with temporal gauge, see eq. 3.13. For left-moving nuclei, the procedure is almost identical, except that the spatial Wilson line starts at . In practice, one simply performs a mirrored computation of the Wilson line starting at some left of the nucleus.
The above procedure allows us to determine the gauge links at some initial time step for arbitrary charge densities given that their longitudinal structure only depends on either or . Since we are solving a second order partial differential equation, we also need to specify the gauge links at the next time step or alternatively the electric field . From the analytical initial conditions we know that the single right-moving nucleus does not depend on . A single time step therefore corresponds to shifting the charge density by to the right along the longitudinal axis. In order to guarantee numerical stability, the time step must be smaller than the smallest lattice spacing, which in our case is usually the longitudinal spacing (see sections 5.3.1 and 5.4.1 for more details on stability conditions). A shift by to the right therefore corresponds to shifting the lattice charge density by a fraction of a longitudinal spacing, which is slightly ambiguous. One solution to this problem is interpolating through the charge density and sampling the interpolated function at the shifted positions. An alternative method, which is more closely related to the CPIC approach, is to sample the charge density with a large number of colored particles as outlined previously. The ensemble of particles (with each particle having a continuous position) can then be easily shifted by to the right. Using the interpolation procedure described in the previous section, we then obtain the charge density in the next time slice. This shifted charge density is used to repeat the computation of the Wilson line . Finally, we determine the transverse gauge links in the next time slice via
[TABLE]
Using and the electric field can be immediately computed from
[TABLE]
The initialization described above has to be performed for each nucleus separately. In order to arrive at the correct field configuration, the nuclei initially have to be far away from each other such that there is no overlap in their color fields. In the implementation, the simulation box is split in half along the longitudinal axis. Each half box is occupied by a single nucleus. We initialize the color fields in each half-box separately. With compact longitudinal support, one simply has to guarantee that the support of each nucleus is contained in their respective half-boxes. In the case of non-compact support (as used in chapter 4, where the longitudinal profile is assumed to be Gaussian) slight errors are inevitable, but strongly suppressed when the initial separation of the nuclei is sufficiently large.
As already mentioned in section 3.1, the single nuclei color fields in temporal gauge are essentially the same as in LC gauge. In particular, they are purely transverse and asymptotically pure gauge behind each nucleus. One therefore arrives at the picture presented in fig. 3.6. The use of temporal gauge forces us to use fixed boundary conditions on the outermost transverse planes, which form the longitudinal boundary.
3.3.3 Temporal gauge charge density
After performing the gauge transformation, the gauge links and color electric fields obey the temporal gauge condition. However, the original charge density used to solve the Poisson equation in each transverse plane is the covariant gauge charge density. To restore gauge covariance one can explicitly perform the lattice gauge transformation of . Alternatively, an even simpler approach is to use the Gauss constraint
[TABLE]
and re-compute the temporal gauge charge density from the electric field and the gauge links. This trivially guarantees that the Gauss constraint is satisfied for the initial conditions. Moreover, any floating point errors introduced by the initialization procedure are absorbed in the temporal gauge charge density. The particle ensemble mimicking the charge density is obtained as described in section 3.3.1.
3.4 Energy-momentum tensor
The most important observables that we want to study (at least for the purposes of this thesis) are the energy density and the pressure components of the Glasma, which are components of the energy-momentum tensor. Since we always consider initial conditions of the MV model type, we do not care about the transverse coordinate dependence of most observables. It is therefore sufficient to consider quantities that are averaged over the transverse plane. We also perform the average over the ensemble of charge densities defined by the probability functional . Due to homogeneity and isotropy (with respect to the transverse plane) of the MV model, the averaged energy-momentum tensor (see eq. 2.149) in the laboratory frame reads
[TABLE]
Here is the energy density, and are the two pressure components, and is the Poynting vector (or energy flux) along the longitudinal direction . The energy-momentum tensor is traceless in the sense of , which yields
[TABLE]
Due to the averaging procedure, all components of are functions of and only. The specific form of is explained by the symmetries of the MV model. For example, there is no net energy flux in the transverse direction and . In the MV model, nuclei have infinite extent in the transverse directions and the average fluctuation of the color charge density is homogeneous. Therefore there is no distinguished direction within the plane. Net energy flux in any of the transverse directions would violate these symmetries. This is not true when considering single events: each collision produces an inhomogeneous distribution of color flux tubes, which expand in the transverse plane. Therefore, we have locally non-vanishing and , i.e. transverse components of the Poynting vector, on an event-by-event basis. However, once we average over all charge densities, these components must vanish. Similar arguments can be made for the shear stress components. These involve expectation values and with . Again, due to rotational symmetry and homogeneity, these terms must vanish.
Analogously to the boost invariant case in section 2.4, the energy density and the pressure can be written in terms of longitudinal and transverse components of the chromo-electric and -magnetic fields:
[TABLE]
In the continuum limit, these four contributions are defined as
[TABLE]
where is the area of the transverse plane. We also use the shorthand and . The Poynting vector or longitudinal energy flux is given by
[TABLE]
On the lattice, we can express these quantities in terms of gauge links. Recall that the plaquette (see eq. 2.101)
[TABLE]
can be used to approximate the squared field strength tensor (see eq. 2.103)
[TABLE]
Here, refers to the center of the plaquette, i.e. . The quadratic error term is of the form , i.e. there are no mixed terms of the form (see e.g. [60] or [63] for a proof). The four energy density components can therefore be written as
[TABLE]
The sum runs over all lattice sites of a given transverse plane. The approximation of field strengths in terms of plaquettes is only accurate up to quadratic order if one considers the field strengths to be evaluated at the centers of plaquettes. For instance, is naturally centered around . If we need to evaluate quantities at positions other than the ones they are naturally defined at, we use averaging. In the case of evaluated at we use
[TABLE]
which is both a spatial and temporal average. Similar averages can be defined for all other components. By using symmetric averages we are able to retain quadratic accuracy for all quantities.
In order to compute the Poynting vector, we need to find approximations of the (un-squared) field strength tensor. This is accomplished by taking the anti-hermitian, traceless part of plaquettes:
[TABLE]
where is the center of the plaquette. This expression follows from inserting eq. 2.102 into eq. 2.119 and expanding for small lattice spacing . Eq. (3.4) has already been used to define the chromo-electric field in the Lagrangian leapfrog method (see eq. 3.62)
[TABLE]
Although transforms locally at the lattice site , the above relation shows that it approximates the electric field in the continuum at up to quadratic order. Similarly, the magnetic field can be approximated via
[TABLE]
where on the LHS are orthogonal directions to . Since and are not defined at the same positions, averaging is also necessary here. In the case of longitudinal energy flux we use
[TABLE]
and
[TABLE]
Here, a time average is used for electric fields and spatial averaging for magnetic fields. The properly discretized version of the longitudinal Poynting vector component from eq. 3.139 then reads
[TABLE]
There are a few things to note: first, the two terms and are not actually defined at the same transverse positions . This poses no problem if we are only interested in quantities averaged over the transverse plane. Secondly, the spatial average of the magnetic field is slightly problematic because it is not gauge-covariant. A better idea might be to use a gauge-covariant spatial average:
[TABLE]
In practice, we found that this does not change results significantly, as long as one does not apply highly discontinuous gauge transformations and the fields can be assumed to vary only slowly from lattice site to lattice site.
A more rigorous expression for the Poynting vector on the lattice can be found the following way: in the continuum, the Poynting vector and the energy density satisfy the Poynting theorem
[TABLE]
if no external currents are considered. This theorem simply follows from the temporal component of the more general energy momentum conservation law
[TABLE]
In a discretized system, one can require that a discrete analogue of eq. 3.155 should hold. Given a discretized version of the energy density , one can then try to find a lattice expression of the Poynting vector, which exactly satisfies the discretized Poynting theorem (see [108] for such a derivation). It is not guaranteed that such an expression for exists for arbitrary discretizations of the local energy conservation law. Furthermore, discretized systems do not even necessarily satisfy local energy conservation exactly. In our numerical studies [1, 2, 3] eq. 3.152 proved to work sufficiently well for our applications.
In some cases, we are interested in components of the energy-momentum tensor without averaging over the transverse plane. For instance, in section 4.2 we study the energy density due to the longitudinal color electric field as a function of the transverse coordinate . We then simply leave out the summation over transverse coordinates. We have to keep in mind where the various energy density components are naturally defined, i.e. in the center of plaquettes.
3.5 Implementation details
During the development of my thesis several implementations have been created based on the numerical procedures presented in this chapter. The first version of the code was built on top of an already existing code base for Abelian particle-in-cell simulations called OpenPixi111The code is publicly available at https://github.com/openpixi/openpixi. The original Abelian PIC simulation can be found at https://github.com/openpixi/openpixi_pic. and co-developed with my colleague Daniil Gelfand and my supervisor Andreas Ipp. This code was written in Java in an object-oriented fashion. This first version was designed to be modular such that components of the simulation (e.g. various solvers for the equations of motion, different implementations of interpolation procedures, initial conditions and gauge groups) can be easily swapped out and replaced. It also features a graphical user interface, which can be used to, for example, monitor observables while the simulation runs or create two- and three-dimensional plots. These features facilitate easier debugging during development and allow users to experiment with the simulation without making changes to the source code. Unfortunately, the original version suffers from performance issues, as the Java programming language is not the best choice for high-performance computing applications. In particular, due to the slightly unpredictable memory management behavior in Java, large simulations are often at risk of running out of available memory. Nevertheless, all results of the first two publications [1, 2] were achieved using the Java version running on the Vienna Scientific Cluster 3 (VSC3). This required great computational resources: for example in [2] the simulation of one single collision took roughly two days to complete. Due to the stochastic nature of charge densities in the CGC framework, one has to perform such simulations multiple times in order to compute expectation values. These long simulation times severely limit the way one can experiment with different initial conditions and parameters. As a result, I wrote a second version222The code is publicly available at https://gitlab.com/monolithu/pyglasma3d. which was implemented in Cython [72] with a focus on a minimal set of features. Using Cython, it was possible to implement non-critical parts of the simulation (such as data in- and output and general program flow) in the Python programming language, while performance critical code (e.g. time evolution, initial conditions) is translated to highly optimized C. This led to a substantial decrease in simulation time by a factor of roughly 7 - 8, while also decreasing memory use by a factor of 6 - 7. The main disadvantage of this code is that it is not object-oriented and tends to be less readable compared to the Java version. Furthermore, due to certain design choices early on in the development, it is also restricted to as the color gauge group. Extension to other gauge groups would require large parts of the code to be rewritten. All results of the most recent publication [4] were obtained using this faster version. A third version333The code is publicly available at https://gitlab.com/openpixi/openpixi_c. written in C++ is currently (at the time of writing) being developed by my colleague Patrick Kappl with comparable performance to the second version, while at the same time maintaining the object-oriented design principle and increased readability of the code [108]. It is also readily extendable to other gauge groups.
In this last section we discuss a few technical details regarding the actual implementation of these codes and some optimizations that were found to improve performance and, in particular memory consumption, of these simulations.
Grouping particles into charge planes
The first thing to notice is that one does not actually have to individually track each particle in the simulation as would be the case in completely general PIC simulations. Since all the trajectories are fixed due to the eikonal approximation (i.e. no recoil) and all particles only move along the longitudinal axis, one can group particles within the same transverse plane into so-called “charge planes”. A single charge plane is defined by its longitudinal position , its velocity or orientation (either left- or right-moving), and the (color) charge density distribution in the transverse plane. Since all particles within a charge plane move collectively, one simply has to check if the NGP associated with the longitudinal position of a charge plane is about to change from one time slice to the next. Computing the current and applying parallel transport is then performed for all charges within the plane. Grouping particles this way avoids iterating over each individual particle and performing unnecessary repeated calculations. Furthermore, as a particle’s position can be trivially inferred from its initial position and the number of already performed simulation cycles, it is also not necessary to update the position as described in section 3.2.3 or actually compute the NGP. It is much easier to keep track of the sub-cell position of each charge plane within each cell, where is the number of particles per cell and its last (integer) lattice position along the longitudinal axis . After each time step the sub-cell position is increased by one for each right-moving plane and conversely decreased for each left-moving plane. At the beginning of each simulation cycle, all right-moving charge planes with and all left-moving planes with will perform a boundary crossing. It is therefore immediately clear that these particular planes contribute to the color current interpolation step and require parallel transport. Accordingly, the lattice positions are updated and the sub-cell positions set to their new values.
Combined current calculation and parallel transport
A further improvement is to combine the current interpolation and the parallel transport into one single step. Recall that current is only generated when particles move from one cell into the next. Parallel transport is also only applied at such boundary crossings. It is therefore natural to perform both calculations simultaneously. If one attempts to do this, it turns out that one has to discriminate between right and left-moving charges: in the computation of the current density of right-moving particles eq. 3.99, one has to use the particle charge of the previous time slice. On the other hand, the current for left-moving particles eq. 3.100 involves the parallel transported charge , which corresponds to the charge in the next time slice , see eq. 3.105. For right-moving particles the combined step therefore involves first computing the color current eq. 3.99 and then applying parallel transport eq. 3.104 immediately after. For left-moving particles the order of steps is reversed: we first apply parallel transport eq. 3.105 and then (using the updated charge ) compute the current density eq. 3.100. This optimization avoids performing parallel transport for left-moving particles twice (which would increase computational effort) or alternatively storing an “old” and “new” version of each particle charge (which would increase memory consumption). This idea can be applied to whole groups of particles as described in the previous paragraph.
Boundary conditions at longitudinal boundaries
As discussed in the previous section, the color fields are non-zero at the longitudinal boundaries of the simulation box due to the choice of temporal gauge. On the lattice, we can implement these Dirichlet-type boundary conditions by adding additional “inactive” transverse planes at the left- and right-most boundary, as shown in fig. 3.7. During the time evolution, these fixed boundaries are simply not updated and therefore the fields within remain constant throughout the simulation. Obviously, the simulation has to be stopped once the nuclei are about to hit these fixed boundaries.
Parallelization with shared memory
A major speed-up is to perform certain calculations in parallel by making use of multiple cores available in modern CPUs. The cores of a typical multi-core CPU can access the same shared memory. In such a shared memory setup, the parallelization of finite difference solvers – such as the real-time lattice gauge theory approach to solving the Yang-Mills equations used in this thesis – are trivially parallelizable. For instance, the computation of the electric field in the next time slice (eq. 3.63 for the Lagrangian leapfrog method or eq. 3.74 for the Hamiltonian leapfrog method) only depends on quantities of previous time slices. Therefore, the calculation of electric fields and at different lattice sites and can be performed independently. The same applies to the update of the gauge links, the interpolation step, parallel transport of charges, and various parts of the initialization procedure. The first version of the code (Java) uses the built-in thread functionality to achieve parallelization. The second version (Cython) and the third (C++) both use OpenMP [109]. All three versions of the code exhibit high degrees of parallelism and scale well with the number of cores used.
Parallelization with distributed memory
Shared memory parallelization is usually restricted to simulations on single computers with multi-core processors. A further speed-up of the calculation can be achieved by employing distributed memory parallelization techniques through a Message Passing Interface (MPI) [110]. In this approach parallelizable tasks are distributed among multiple computers communicating e.g. via network interfaces. Since each individual computer (or node) can only access its own memory directly, such approaches are referred to as parallelization with distributed memory. Communication between nodes is achieved by passing and receiving messages. In the 3+1D Glasma setup a straightforward way of splitting up the computation is to divide the simulation box along the longitudinal coordinate. Each node is assigned one of the smaller simulation boxes as seen in fig. 3.8. Analogous to how fixed boundary conditions were implemented using additional transverse planes, we use communicating boundaries to synchronize the individual smaller boxes. This introduces an additional synchronization step at the end of each simulation cycle. Additionally, when particles cross into communication layers, they are removed from the originating node and reintroduced in the communication layer of the neighboring node. In practice, whole charge planes (i.e. all charges within a transverse plane grouped together) are transfered if such a node crossing happens.444This simple transfer of particles (or charge planes) only works with the NGP scheme. Using higher order schemes (e.g. the CIC scheme), it might be necessary to keep copies of particles in two adjacent simulation boxes simultaneously during an extended transition step. At the time of writing, only the Cython version of the code is able to perform parallel calculations in this manner.
Chapter 4 Single color sheet approximation
Having developed the necessary numerical methods in the previous chapters, we are now finally able to discuss heavy-ion collision simulations in 3+1 dimensions. Much of the material covered in this chapter is based on the first two publications [1, 2] and conference proceedings [3]. In addition, we include extended discussions and some new results.
Our starting point is a simple regularization of the original MV model. In the boost invariant MV model (see section 2.1), the (covariant gauge) color current of a right-moving nucleus is given by eq. 2.2
[TABLE]
with the associated gauge field (2.5)
[TABLE]
The gauge fields inherits the longitudinal profile of the color current (4.1). Explicitly separating the dependence, we can write
[TABLE]
where describes the color structure in the transverse plane. Inserting this into the Yang-Mills equations (2.3), we find
[TABLE]
see eq. 2.11. The transverse charge density obeys a Gaussian probability functional defined by the one- and two-point functions (see eqs. (2.29) and (2.30))
[TABLE]
As already mentioned in section 2.1, this model neglects structure along the longitudinal axis. We therefore refer to this as the single color sheet approximation, as the nucleus resembles only a single infinitesimally thin sheet of color charge defined by the transverse color charge density .
The first step towards three-dimensional collisions using this model is to express everything in terms of laboratory frame coordinates. In this frame we have . Analogously, for the current in the laboratory frame we find
[TABLE]
where we have used for some constant . Similarly, we have
[TABLE]
In order to explicitly break boost invariance, we choose to regularize the delta distribution and replace it with a smooth longitudinal profile function . Our ansatz for the single color sheet model in 3+1 dimensions therefore reads
[TABLE]
with and . We choose to be a Gaussian given by
[TABLE]
where controls the longitudinal width in the laboratory frame. It holds that
[TABLE]
In contrast to the singular color fields of the boost invariant model, the regularized color fields are smooth along and can be resolved on a finite lattice. The boost invariant limit is restored by taking . It is obvious that such a limit cannot be taken in practice, due to the finite resolution of the lattice. In section 4.2 we will see that numerical accuracy is lost for very thin nuclei at fixed lattice spacing.
One can relate to the Lorentz contracted width of the nucleus, where is the nuclear radius and is the Lorentz factor. However, since the Gaussian profile does not have compact support, there is some ambiguity in choosing a proportionality factor for . In [1, 2] we simply set
[TABLE]
such that practically all of the color charge density lies within the longitudinal width . The Lorentz factor can be directly computed from the collision energy per nucleon pair via
[TABLE]
where is the nucleon mass. If we consider the MV model parameter and the Yang-Mills coupling to be constants, then the energy dependence only enters via the thickness parameter .
The regularization of the delta distribution in eqs. (4.9) and (4.10) allows us to compute the total field energy of a single nucleus analytically. From section 2.1 we already know that the color-electric and -magnetic fields of single nuclei moving at the speed of light are purely transverse and, after averaging over color charge configurations, give the same contribution to the energy density. We can therefore write
[TABLE]
In covariant gauge we find for the transverse electric field and thus
[TABLE]
First, we notice that due to the presence of the term, the boost invariant limit, where we take , is ill-defined. Fortunately, this is what one would expect. Integrating the energy density over and and using the Gaussian profile for , we end up with
[TABLE]
The total energy of a nucleus thus diverges with when taking . From our previous considerations we know that . The ill-defined boost invariant limit is therefore related to ultrarelativistic nuclei carrying infinite energy.
The next step is to compute the expectation value , which can be easily related to the two-point function . From the infrared regulated two-point function section 2.1 we can read off
[TABLE]
Using integration by parts and the fact that the MV model is invariant under translations in the transverse plane we can use
[TABLE]
When plugging in the result for the two-point function from above however, one quickly realizes that this does not yield a finite result. Even though the two-point function is well defined through infrared regularization, taking the second derivative does not work at the origin . For the energy density in the laboratory frame, ultraviolet regulation is essential. We therefore compute the required expectation value directly in momentum space:
[TABLE]
In the last line we used . The total field energy then reads
[TABLE]
where the integral over is replaced with the transverse area .
In the single color sheet approximation it is easy to compute the temporal Wilson line for transforming from covariant gauge to temporal gauge as required for the numerical method. The temporal Wilson line (see eq. 3.7) is given by
[TABLE]
where is given in covariant gauge. Using we see that the color structure along the longitudinal direction is trivial, because the profile function simply smears out the color field in . As a consequence, we can drop the path ordering and find
[TABLE]
where the integrated profile function is given by
[TABLE]
For , i.e. far behind the right-moving nucleus, we have and recover the asymptotic pure gauge solution known from the original MV model with a single color sheet (see eq. 2.40 for )
[TABLE]
The fact that the Wilson line can be computed analytically greatly simplifies the initialization procedure of the simulation. Instead of actually constructing the path ordered Wilson line on the lattice as detailed in section 3.3.2, in the single color sheet approximation we can simply evaluate eq. 4.23 at each point of the lattice to compute the temporal gauge initial conditions.
It is noteworthy that the Gaussian profile function does not have compact support along the longitudinal axis as often assumed in the previous chapter. In practice this does not pose a problem, because the Gaussian profile tends to zero exponentially fast. It is sufficient if the two nuclei are well separated at the beginning of the simulation. Any errors can be suppressed below the floating point accuracy used in the code and are therefore negligible.
4.1 Simulation parameters
Before presenting numerical results, we summarize all parameters that define the initial conditions and the simulation itself. First, there are a number of physical parameters related to the initial conditions:
- •
The coupling constant fixes the strength of interactions. The canonical value in CGC simulations is , which can be obtained by evaluating the strong coupling constant at the relevant energy scale .
- •
The MV model parameter , which controls the average color charge density of the nucleus: we use the estimate provided by McLerran and Venugopalan [46, 47] which yields approximately in the case of gold nuclei (see also sections 2.1 and 2.4).
- •
The longitudinal thickness : it can be related to the collision energy per nucleon pair via (see eqs. (4.13) and (4.14)) where is the nucleon mass and is the nuclear radius. For collisions at RHIC with we have .
- •
Infrared modes are regulated by . Usually this is set to , such that we can mimic the effects of confinement. In order to study the dependence of this parameter on collisions, we vary this up to a few .
- •
Ultraviolet modes are cutoff by , which is set to values between and . For some applications we leave out manual UV regulation and let the modes be regulated by the transverse lattice spacing .
On the other hand, there are parameters related to the lattice and size of the simulation.
- •
The longitudinal and transverse system size and fix the finite simulation volume. The transverse area (spanned by and ) is assumed to be square, i.e. . Periodic boundary conditions are used in the transverse directions, which leads to a system size dependent infrared cutoff proportional to the inverse of . One therefore is effectively dealing with two separate infrared regulators. In order to not have any interference between these two regulators, it is necessary to set the screening mass to larger values . The finite length along the longitudinal axis limits the maximum time the simulation can run before the colliding nuclei approach the fixed longitudinal boundaries of the box. If the two nuclei are initially placed at and , they will overlap completely at and and have switched positions at , which is the typical time where we stop the simulation. If the nuclei are very thin compared to , one can keep the simulation running longer while making sure that there is no interference from the boundary conditions.
- •
The lattice size is defined by the number of cells along the longitudinal and transverse directions and . The lattice spacings are given by and . In practice, we work with dimensionless quantities in the simulation. This is done by choosing an energy scale defined by the transverse lattice spacing . All dimensionful quantities are then expressed in terms of . In particular, the transverse lattice spacing can be set to and the longitudinal spacing becomes the length ratio of and . There are two transverse scales that have to be properly resolved on the lattice. First, there is the screening range induced by the infrared regulator . It roughly translates to the size of color neutral regions in the transverse plane. In order to get close to the continuum limit with our lattice calculation, the transverse spacing must be much smaller than . However, there is a second (usually) even smaller length scale provided by which corresponds to the transverse size of Glasma flux tubes produced in the collision. Obviously, the simulation must be able to resolve this scale as well and therefore must be chosen such that it is smaller than the average Glasma flux tube. This leads to a hierarchy of scales on the lattice: . The longitudinal spacing must be small enough in order to resolve the full longitudinal width of the nucleus and possibly also a detailed structure within the longitudinal extent. It also acts as an effective momentum cutoff proportional to . The single color sheet approximation, which has a very smooth color structure along , does not require a highly resolved longitudinal direction.
- •
The time step or must be strictly smaller than the smallest lattice spacing, which is usually . For the methods developed in chapter 3, numerical stability is guaranteed for (see sections 5.3.1 and 5.4.1 for a derivation of the stability criterion). The NGP scheme for particles works best if we set to a fraction of such that the number of particles per cell yields a whole number. This way color current is generated at each time step and no additional numerical noise is introduced into the simulation. In chapter 5 we develop a numerical integrator that allows us to set , but for the results presented in this chapter we use the methods of chapter 3.
4.2 Approaching the boost invariant limit
The result of a typical collision simulated with our code is shown in footnote 2, where we plot the energy density of the color fields. As is evident from the figure, it is possible to study the fully three-dimensional structure of the Glasma.
The first step we take is to validate that the simulation produces meaningful results when compared to the boost invariant case. As discussed earlier, the ultrarelativistic limit should be reached by setting . It is clear that for any fixed lattice resolution this limit cannot be taken directly. If is comparable to , then the nuclei are not properly resolved on the lattice and the numerical solution will suffer from lattice artifacts. Nevertheless, we can check if the Glasma created in a 3+1 dimensional collision resembles the boost invariant limit by choosing to be very small and comparing the fields from our simulation to the boost invariant initial conditions described in section 2.2.
The question then arises what quantities one should actually compare. In the boost invariant case the initial conditions can be formulated as initial values for the longitudinal chromo-electric and magnetic fields as given by eqs. (2.76) and (2.77). These fields are given in the co-moving temporal gauge , which is different from the laboratory frame temporal gauge condition . Even though both gauge conditions agree at mid-rapidity , it is a better idea to compare gauge invariant quantities. We therefore choose the longitudinal electric part of the energy density
[TABLE]
as the main observable. Given the lattice charge densities of two colliding nuclei, we can compute at directly from lattice initial conditions in the boost invariant case, see eq. 2.136. On the other hand, in the 3+1 dimensional setup we use the same charge densities but smeared out along with some finite width and record the longitudinal chromo-electric field on the lattice. One has to decide when and where the two results should be compared. For simplicity, we choose the center of the simulation box at , which corresponds to the mid-rapidity region . This is also the position at which the two nuclei will overlap completely if their initial placement was symmetric around the center . Since the interaction of the nuclei during the collision is not instantaneous anymore, the time at which the fields should be compared is ambiguous. A natural choice would be when the maximum overlap takes place, but for very thick nuclei the collision starts much earlier when the color fields start to overlap. We therefore record the longitudinal electric field in the transverse plane at as a function of time to study this effect. The transverse structure is shown in fig. 4.2, where we see that for small the 3+1 dimensional setup generally agrees with the boost invariant fields (at least visually).
In order to make a more meaningful and quantitative comparison, we compute the correlation coefficient between the boost invariant, “analytical” result (as determined from eq. 2.136) and the numerical result from our 3+1 dimensional simulations. The correlation is given by
[TABLE]
where the covariance across the transverse plane is defined by
[TABLE]
and \overline{\mathrm{tr}\big{[}E^{2}_{L}\big{]}_{\mathrm{ana}}}, \overline{\mathrm{tr}\big{[}E^{2}_{L}\big{]}_{\mathrm{ana}}} are the mean values of the energy density in the transverse plane
[TABLE]
The standard deviations , are given by
[TABLE]
Using has the advantage that it only measures how similar the two energy densities are without taking their absolute values into account. Lower collision energies produce a less energetic Glasma compared to the ultrarelativistic case, but one can still expect similar distributions of flux tubes in the transverse plane, especially for small .
Fixing some initial two-dimensional charge densities and for the two colliding nuclei, we can run the simulation for various values of and compute the correlation in each case. Since the correlation depends on when we make the comparison, we record as a function of . The results of this are shown in in fig. 4.3 (a): Here, the blue curve corresponds to where the two nuclei overlap completely. Starting from very large values of the thickness parameter, the correlation increases as we decrease . However, at some point the correlation is quickly lost and decreases again. In the plot this happens around . As mentioned earlier, the thickness has to be resolved properly on the lattice in order to obtain results similar to the continuum limit. For this is not the case anymore. In the case of these very thin, improperly resolved nuclei we observe a numerical instability. Even before the collision, the color fields of the single nuclei start to disperse and lose their original longitudinal profile. At the same time, unphysical longitudinal fields start to show up. We will come back to this instability and provide a solution in chapter 5. For the time being, we will make sure not to go below a threshold of for the thickness parameter. Unfortunately, we will see that is not the only parameter that affects the development of this instability.
The red curve in fig. 4.3 (a) is the maximum value of the correlation recorded during the collision. Even in the case of large , relatively high values can be achieved, although at earlier times than the overlap . For very thick nuclei the color fields start overlapping sooner, producing the characteristic longitudinal Glasma fields. These initial fields therefore start evolving much earlier and consequently have changed by the time both nuclei overlap.
Furthermore, we also study how infrared and ultraviolet regulation affects our results while keeping fixed. This is shown in fig. 4.3 (b), where we plot the correlation for as a function of and two choices of the UV regulator. The solid lines correspond to which is below the lattice cutoff . The dashed lines correspond to . The same infrared and ultraviolet regulators are used in the boost-invariant scenario. There are a few things to notice here: the manually regulated simulations are consistently above the lattice regulated cases. The use of the ultraviolet regulator helps us achieve more accurate results. The infrared regulator affects the correlation only slightly except when going to very small values. Using small generally leads to large field amplitudes as is evident from the solution of the Poisson equation (3.122). We observe that large field amplitudes lead to similar unstable field configurations as in the case of small thickness . Our results are therefore generally unreliable when we study thin nuclei with large color fields without appropriate lattice resolution. The effects of the instability can be suppressed or postponed to later times when we use very fine lattices, but only at the cost of larger computational effort.
Studying the correlation between our numerical results and the analytic expressions for the boost invariant initial conditions shows that we are able to correctly describe boost invariant collisions in the limit of thin nuclei. It also reveals that one has to be careful in choosing simulation parameters, in particular .
4.3 Pressure and energy density
In this section we study the energy density and the pressure components of the Glasma in more detail. We are particularly interested in deviations from the 2+1 dimensional Glasma and the evolution of the system in the laboratory frame.
4.3.1 Pressure anisotropy
As discussed in sections 2.4 and 2.5, a prominent phenomenon in the early stages of heavy-ion collisions is the pressure anisotropy and a possible subsequent isotropization of the system. The main observables in this context are the transverse and longitudinal pressure components and , averaged over the transverse plane. On the lattice, we compute these observables as described in section 3.4.
Plots of the (averaged) pressure components in the laboratory frame at different times are shown in fig. 4.4. Before the collision, the field strengths of the colliding nuclei are purely transverse and, as such, only contribute to the longitudinal pressure component. Once the collision starts, longitudinal color-electric and -magnetic fields are produced and transverse pressure builds up. This pressure is much smaller than the longitudinal pressure component, because the large color fields of the nuclei dominate the Glasma fields. After the collision, longitudinal pressure remains largely flat. On the other hand, falls off exponentially towards . In this central region we observe a much larger compared to the , which already indicates the characteristic pressure anisotropy of the Glasma.
We now turn towards studying this pressure anisotropy. For the further analysis it will be sufficient to stay in the central region or . Recall from section 2.2 that in the boost invariant scenario the field strengths at are purely longitudinal, which leads to maximally anisotropic initial pressures and . During the expansion of the Glasma flux tubes, transverse color-electric and -magnetic fields are created until and . Once the system settles, it will remain in the free streaming limit. As discussed in section 2.5.1, the free streaming limit is incompatible with the subsequent (isotropic) hydrodynamical evolution of the QGP. Introducing boost invariance breaking fluctuations drive instabilities in the Glasma, which can move the system towards isotropization [39, 41, 42, 86], albeit more slowly than required by the scenario of fast isotropization. In the 3+1 dimensional setup we explicitly violate boost invariance by introducing a finite nucleus thickness. It is therefore interesting to investigate the effects of the thickness parameter on the pressure anisotropy of the Glasma.
For the numerical study of the anisotropy in the Glasma we use the pressure to energy density ratios and with , as introduced in section 2.4 (see fig. 2.7). The free streaming limit then corresponds to and . Isotropization would be signaled by . Both the pressure and energy density components are averaged over the transverse plane and events are used for the statistical sampling. We choose a grid size of and . For collisions of thick nuclei in fig. 4.5 we choose a lattice spacing of , which corresponds to . The transverse area then almost covers the full area of a gold nucleus. For simulations of thin nuclei in fig. 4.6, we are forced333These simulations have been carried out in spring of 2016 with an early build of the Java version. Due to slower performance and high memory requirements, lattice sizes of were at the limit of available computational resources (see section 3.5). to use smaller lattice spacings of (for ), (for ) and (for ). The transverse area then only covers , and of the full area respectively. The time step is set to .
There are several observations we make regarding the results shown in figs. 4.5 and 4.6:
- •
From fig. 4.5 we see that we recover the free streaming limit of the boost invariant case. Isotropization is not reached within possible simulation times due to limitations from the longitudinal length of the simulation box. We observe slight movement of both pressure components towards the desired value of , but not within any realistic time scales.
- •
The initial pressures directly after the collision behave differently compared to the boost-invariant case. In our simulations of thick nuclei in fig. 4.5 we see that in the beginning dominates due to the presence of the large transverse fields of the colliding nuclei. As the nuclei recede from the collision center, the Glasma fields have already reached the free streaming limit and therefore no negative longitudinal pressures are observed.
- •
In the results for thin nuclei in fig. 4.6 we recover negative longitudinal pressure. The colliding nuclei move away from the collision center fast enough, leaving behind longitudinal color flux tubes, which have not decayed yet. The still largely longitudinal fields generate negative pressure, which is characteristic for the early Glasma phase. We remark that the infrared regulator only weakly affects the pressure anisotropy.
4.3.2 Mismatch between electric and magnetic fields
In the following, we investigate the production of longitudinal chromo-magnetic fields and -electric fields characteristic for the Glasma at early times. In the boost invariant case the contributions to the energy density from magnetic and electric color flux tubes ( and respectively) should be roughly equal after averaging over initial conditions, see fig. 2.5. In the 3+1D setup with finite one observes that this is not the case and there is a dependency on the thickness parameter as well as the IR regulator . The results are presented in figs. 4.7 and 4.8.
Figure 4.7 shows the ratio of magnetic and electric longitudinal fields in the central region ( or ) for a range of values of the nucleus thickness and an IR regulator . Collisions of thick nuclei show a very small ratio of about after the collision. In the case of thin nuclei in fig. 4.7 (b) the ratio increases to roughly , which is still far away from the “canonical” value of in the boost invariant scenario. One observes that for early times the equal ratio is indeed realized. This is likely due to the fact that the first interactions of color fields (when the nuclei are about to overlap) are very similar to the boost-invariant case. Compared to the ultrarelativistic limit, the interaction of nuclei in our simulation takes place over a finite interaction time. The overlapping of the two colliding nuclei starts with the first few “layers” of the color fields. In this case, the elapsed interaction time is only very short, which could explain why we see similar results to the ultrarelativistic limit. The same effect was observed in section 4.2 where the correlation coefficient reached its maximum much earlier than the complete overlap of the two nuclei. Note that due to the small physical volumes used in the simulations of thin nuclei it is harder to achieve adequate statistics. As a result, the curves in fig. 4.7 (b) are not as smooth as in fig. 4.7 (a).
The results do not only depend on the thickness . In fig. 4.8 the results are shown for a fixed nuclear thickness () and a varying IR regulator. We observe that reducing to (which roughly corresponds to a correlation length of the color fields of the order of the confinement radius ) leads to better agreement with the boost-invariant case with a ratio of . This strong dependency of the ratio of magnetic and electric longitudinal fields on the IR regulator is not present in the boost invariant case.
The presented results seemingly suggest a suppression of chromo-magnetic flux tubes (or an overproduction of chromo-electric flux tubes) in the Glasma phase when introducing a finite nucleus thickness. However, the strong dependency of the magnetic to electric longitudinal field ratio on the IR regulator leads us to suspect that this discrepancy between our simulations and the boost invariant case is an artifact which can be attributed to the single color sheet approximation. In covariant gauge the three-dimensional color charge density is given by (see eq. 4.9)
[TABLE]
The charge density is “coherent” along over a length of . Switching to temporal gauge changes the charge density and modifies this longitudinal coherence slightly, but the nucleus largely remains smooth in the direction of . This stands in contrast to the random longitudinal structure of realistic nuclei and the continuum limit of the generalized MV model as described in section 2.1. On the other hand, the IR regulator enforces color neutrality on transverse length scales of . Consequently, the typical color structures in the initial conditions have a thickness proportional to and a transverse width of the order of . To be consistent with the picture of a highly Lorentz contracted nucleus modeled by classical Yang-Mills fields, one would demand that , such that nucleons within the nucleus are also contracted to flat “pancakes”. Therefore, if we move away from the limit , we can expect to see deviations from the boost invariant case, but these deviations may very well solely be due to the longitudinal coherence. This reasoning is consistent with our simulation results: In the case of a very thick nucleus with and the longitudinal magnetic fields are weakened as seen in fig. 4.7 (a). Here we have a value of , which corresponds to color structures which are prolonged in the longitudinal direction. We can compare this to the case of and as presented in fig. 4.8. The ratio of magnetic to electric fields is closer to and at the same time we have , which can be considered small. However, this is merely an observation and not a physical explanation for the suppression of the longitudinal chromo-magnetic field.
4.3.3 Non-conservation of total energy in the simulation
One of the fundamental assumptions made in the CGC framework is the separation of hard and soft degrees of freedom, which are modeled as external color charges and classical gauge fields respectively. As a result of the collision there is an energy exchange between the charges and the fields. However, since the nuclei are assumed to be recoilless, the hard sector acts as an inexhaustible energy reservoir for the gauge fields. The resulting field energy increase can be interpreted as the work done by the charges against the field. In the boost invariant case this effect is implicitly included in the initial conditions for the fields at . The total increase in energy is less apparent, because the 2+1 dimensional setup is restricted to the future light cone and the charges are not taken into account.
In our approach we are able to explicitly compute the energy increase during and after the collision. The change of the total field energy density as a function of time can be formulated in terms of an energy continuity equation
[TABLE]
which is the non-Abelian version of the Poynting theorem. The time dependence of is governed by two terms: the Poynting vector and . The integral over the divergence of the Poynting vector can be omitted in the continuum assuming that there is no net energy flux at asymptotic infinity. On the lattice this term vanishes due to the use of periodic boundary conditions. In the scalar product the only non-vanishing part of the current is the longitudinal component and therefore the expression reduces to . Consequently, the energy production is caused by longitudinal chromo-electric fields in the Glasma and must be centered around the collision event and the boundary of the forward light cone where the color currents are nonzero.
The energy increase as a function of time is shown in fig. 4.9. We observe that the total energy density is conserved before the onset of the collision when the external charges and the classical fields describing both nuclei are propagating through vacuum. There is a strong energy increase during, as well as after the collision. At later times there is an ongoing, but slowly decreasing energy production, which finally becomes almost constant. To check the stability of our results with respect to a change in the spatial and temporal resolution of the grid, we vary the spatial lattice spacing and the time step . Overall, there is a good agreement between results at different discretizations. The violation of eq. 4.32 is small and can be further reduced by using smaller time steps as can be seen in the lower plot of fig. 4.9 (b).
4.4 Rapidity profiles and broken boost invariance
In section 4.3 the focus was on observables evaluated directly in the laboratory frame and particularly at mid-rapidity. In the 3+1 dimensional setup the invariance under boosts is explicitly broken and therefore one expects that observables should depend on rapidity. Transforming quantities, such as the energy momentum tensor, from coordinates to the co-moving frame is conceptually simple, but there are some technical difficulties associated with this procedure that should be addressed properly.
First, in order to define the transformation from to , one has to choose a coordinate origin. While this is trivial in the case of infinitely thin nuclei, the introduction of finite thickness makes this choice more ambiguous. The simplest idea is to choose the time and position of complete overlap as we did in the previous sections. However, considering the results of fig. 4.9, it is clear that the production of the Glasma fields has not even reached its maximum by that time. One can therefore try a different approach and choose the origin dynamically as the space-time coordinates of maximum transverse pressure . The reasoning for this is as follows: the transverse pressure is related to the energy density contributions from longitudinal fields . Therefore, the produced Glasma is the only source of transverse pressure. As the Glasma evolves the initially longitudinal fields decay and decreases. therefore exhibits a maximum at close to the overlap event , but is slightly later as the generation of the Glasma is not instantaneous. Furthermore, in symmetric collisions (e.g. two gold nuclei in the center-of-mass frame) the spatial coordinate at overlap must be the same as on average. The origin defined by the maximum of is therefore naturally centered, but slightly later than the overlap which accounts for the slowed creation of the Glasma. Having defined the origin, the laboratory frame coordinates can be shifted such that . The future light cone is then parametrized using
[TABLE]
There is also another problem in the 3+1 dimensional setup: due to the finite thickness of the nuclei, there is no clear distinction between the fields of the Glasma and nuclei. This is particularly problematic when evaluating observables at high values of rapidity at early times . Given some directly after the collision, high absolute values of lead to coordinates that can end up within the receding nuclei. Since the color fields of the nuclei are generally much larger than the Glasma fields, getting too close to the nuclei strongly affects the observable one intends to compute. In practice, one is therefore restricted to much smaller, “safe” rapidity intervals, which depend on . We experimented with subtracting the fields of the nuclei in order to circumvent this problem and extend available rapidity intervals. Due to the large difference in size of the nuclear and Glasma fields, we found that this is a highly delicate matter and generally, subtracting did not lead to an overall improvement of results.
A further problem is related to lattice resolution in the laboratory frame: in the 3+1D setup we use a regular hypercuboid grid with lattice spacings . Observables, such as the energy momentum tensor, are evaluated at discrete positions and times and (the transverse coordinate is ignored for the moment). Transforming to the co-moving frame involves the smooth mapping of eqs. (4.33) and (4.34), which distorts the regular lattice. Therefore, before performing the coordinate transformation, we compute interpolation functions that can be readily evaluated at arbitrary and . However, one should be careful when evaluating observables at larger using this method. Finite resolution in the laboratory frame obviously translates to finite resolution in and . This implies that for some resolution is highest at mid-rapidity and decreases for higher values of . In order to still obtain reasonable numerical results at large , the longitudinal coordinate must be resolved by a very small lattice spacing with and .
4.4.1 An improved dependent MV model
Here, we revise the MV model used in sections 4.2 and 4.3. Until now the only parameter influenced by the collision energy is the longitudinal thickness via
[TABLE]
where is the nucleon mass and is the nuclear radius with nucleon number . In particular, we will be interested in Au-Au collisions as performed at RHIC with and . The radius of a gold nucleus is given by . The thickness at the relevant collision energies is for and for .
The saturation momentum should depend on as well. In order to quickly derive the relation between and we follow the reasoning presented in [50, 69, 111]: deep inelastic scattering (DIS) experiments at HERA (i.e. probing hadrons with leptons, for example using proton-electron collisions) have shown the phenomenon of “geometric scaling” [112]. Without going into too much detail, geometric scaling establishes a relationship between the saturation momentum of the proton and the longitudinal momentum fraction of a parton of the proton scattering off an electron:
[TABLE]
The values of the parameters , and can be obtained from fits to experimental data. In [113] it was found that , and . Assuming an average effective longitudinal momentum fraction given by (see e.g. [50])
[TABLE]
one can solve for to obtain
[TABLE]
where and is given in GeV. The saturation momentum in large, heavy nuclei can be related to the saturation momentum of the proton via
[TABLE]
where is the nucleon number and is a nuclear modification factor for which there exist a number of various models and fits (see the discussion section of [50] for a short summary and references). Assuming and for gold nuclei, one finds
[TABLE]
Depending on which fits are used and what one chooses for , the numerical factors can change slightly. Similar scaling behavior of can also be shown theoretically: within the CGC framework both the JIMWLK and the BFKL evolution equations (which are the weak field limit of the JIMWLK equations [114]) exhibit geometric scaling and predict energy dependence of [67, 115]. For simplicity, we use the fitted estimate of eq. 4.40 and apply it to the MV model. In our publication [2] a slightly different estimate was used:
[TABLE]
where is given in GeV and the numerical pre-factor was dropped. At RHIC we would therefore have and . Since the saturation momentum defines the relevant energy scale for the collision, it also determines the appropriate Yang-Mills coupling constant . Using an approximation of the one-loop running coupling (see [116] p. 141)
[TABLE]
where and active fermion flavors, we can estimate the coupling constant at . With this yields for and for . Given and we can determine the MV model parameter . To do this, we adopt the relationship from [24]. This yields for and for . All parameters are summarized in table 4.1. Regarding infrared and ultraviolet regulation we use for all simulations and . For the following simulations we use a numerical grid of size and a simulation volume , which implies . The time step is set to .
As before, we record the energy density contributions from chromo-electric and -magnetic fields as defined in section 3.4. In fig. 4.10 we show these quantities evaluated at mid-rapidity or as a function of time . The coordinate origin is determined by the maximum of transverse pressure . These results are for SU(2) and have not been rescaled using -dependent factors. In addition, we show the transverse pressure in fig. 4.11, which shows a pronounced maximum close to the collision event.
4.4.2 Energy density rapidity profiles
During the simulation we record the components of the energy-momentum tensor in the laboratory frame and average over a number of collision events to obtain the expectation value . In the MV model most of components of vanish after averaging over all initial conditions due to homogeneity and isotropy in the transverse plane. Therefore, the energy-momentum tensor reduces to (see section 3.4)
[TABLE]
where is the energy density, and are the longitudinal and transverse pressure components and is the longitudinal component of the Poynting vector. If we are interested in the energy density in the co-moving frame we perform the projection
[TABLE]
where is the timelike four-velocity of the co-moving observer. Plugging in everything, we find
[TABLE]
However, this is not the quantity that we will study for two reasons: first, this energy density is sensitive to the exact placement of the coordinate origin. As mentioned in the introduction of this section this choice is slightly ambiguous. Secondly, it is a priori not obvious that the frame is appropriate for the Glasma with broken boost invariance. A simpler and coordinate-free definition of the Glasma energy density is the local rest frame (LRF) energy density. The LRF corresponds to the frame in which the energy flux (or Poynting vector) vanishes. We can write a general Lorentz transformation of as a similarity transformation via
[TABLE]
or using matrix notation
[TABLE]
The matrix associated with an arbitrary boost in the direction is given by
[TABLE]
with and . Applying the transformation and requiring that the energy flux in vanishes in the new frame, we find the velocity of the LRF
[TABLE]
The energy density in the LRF is given by
[TABLE]
The above result can also be obtained by directly diagonalizing the energy-momentum tensor as and are related via a similarity transformation. The LRF energy density can be computed at any point without explicitly referring to the coordinates origin of the frame. However, since we are specifically interested in the rapidity dependence of observables, we compute everywhere in the future lightcone and then, using proper time and space-time rapidity , plot the profile of the LRF energy density as a function of for various values of .
In fig. 4.12 we see a calculation of the rapidity profile of the local rest frame energy density for a RHIC-like collision: choosing parameters and , we obtain an approximate Gaussian profile of the local rest frame energy density, i.e.
[TABLE]
in the range . Using a Gaussian fit we compute and extrapolate to higher . At higher values of the IR regulator the rapidity profiles become more narrow, which seems to suggest that the transverse size of the color-neutral regions has an effect on broken boost invariance. In the CGC literature it is well-known that quantities such as the initial energy density or the gluon multiplicity can be sensitive to the choice of the infrared regulator [39, 49, 50, 55]. Here we add another example of a strong dependency on the infrared regulation.
Even though the results of our model should be regarded as more qualitative than quantitative, we compare our findings to experimental results: it is an interesting observation that the rapidity profile of the energy density for agrees with the measured rapidity profile of pion multiplicities for the most central collisions at RHIC. Of course, a direct comparison of and profiles is not strictly valid: the gluon number distribution can be somewhat broader than the energy density [118]. The correct approach would be to use our results as initial conditions for hydrodynamic simulations in order to make a more direct connection with measured observables. The subsequent hydrodynamic expansion of the system likely increases the width of the profiles further as mentioned in [34], which would favor the more narrow curves (b) and (c) in fig. 4.12. Under these assumptions, the width of the rapidity profile of measured charged particle multiplicities can be seen as an upper limit for realistic rapidity profiles computed from our simulation. We also compare to the Gaussian rapidity profile of the hydrodynamic Landau model [119] with .
The profiles have been computed at , which roughly corresponds to the transition from the Glasma to the QGP. In general, we observe that the rapidity profiles quickly converge for and afterwards become independent of . We also observe free streaming behavior signaled by for a wide range of and longitudinal velocities . The longitudinal velocity can be computed directly from eq. 4.49 and is shown in fig. 4.13, where we plot in the case of and . We see that the velocity of the Glasma agrees with the longitudinal velocity of the boost-invariant Glasma . This implies that there is only negligible flow of energy between different rapidity directions and that the LRF of the 3+1 dimensional Glasma essentially corresponds to the co-moving frame. Our results suggest that the rapidity dependence of the Glasma is fixed early on in the collision and remains unchanged thereafter. To see the effects of increased longitudinal thickness, we repeat the calculation of the LRF energy density for in fig. 4.14. We fit to a Gaussian shape and find that, as expected, the profiles become narrower as compared to the case. This time, our results are also narrower than Landau’s prediction.
Investigating the reason behind the non-boost-invariant creation of the Glasma, we look at the space-time distribution of the transverse pressure in and along the forward light cone in fig. 4.11. The transverse pressure is solely due to the presence of longitudinal fields, i.e. is equivalent to the energy density generated by longitudinal fields . In the boost invariant case the longitudinal fields would be constant along the boundary of the light cone as determined by the boost invariant Glasma initial conditions at . In our simulations we see very different behavior: the longitudinal fields are mostly centered around the maximum in the extended collision region and decrease rather quickly along the boundaries. Since the initially longitudinal fields are the starting point of the evolution of the Glasma, this sharp decrease means that for some fixed proper time there is less Glasma at higher values of rapidity , leading to the observed Gaussian profiles of the energy density. Furthermore, we see a mostly constant spatial distribution of in fig. 4.11 for later times . A similar distribution has been found in holographic models of heavy-ion collisions [120, 121, 122].
4.4.3 Extension to larger rapidity ranges
In the previous subsections we computed the LRF energy density only within rather restricted rapidity intervals . The multiple reasons for this restriction were already outlined in the beginning of section 4.4, but here we present a method that at least circumvents the problem of reduced lattice resolution at higher values of rapidity.
This method is based on the fact that, at least in principle, the dynamics of Yang-Mills fields are invariant under the group of Lorentz transformations, since the underlying action is written in a Lorentz covariant manner. On the lattice Lorentz symmetry is broken, but is restored when taking the continuum limit. Therefore, if the parameters of the simulation are chosen well such that the numerical results are good approximations to the continuum system, the following is possible: first, we perform simulations of a symmetric collision as previously and record observables (such as the energy-momentum tensor). Next, we take the same initial conditions as before, but apply a longitudinal boost of some rapidity , effectively changing the widths of the nuclei and resulting in an asymmetric scenario. Again, we simulate the collision and record observables. Finally, due to Lorentz symmetry, it is possible to undo the Lorentz transformation of the boosted collision and find that the observables agree in both cases (symmetric and asymmetric collisions) up to small violations due to the lattice discretization.
We can exploit this property to extend the available “safe” rapidity intervals in the 3+1 dimensional setup. Consider a simulation of a symmetric collision, where the energy-momentum tensor can be computed within up to sufficient accuracy. Performing a boost by e.g. , this interval is shifted to (in terms of of the original frame). Observables can now be computed in this new interval with the same resolution as in the un-boosted case. Due to the overlap of intervals in , it is possible to check whether boosting leads to any violations of Lorentz symmetry of the discretized Yang-Mills system. In the extended interval observables at higher rapidity can be extracted. This process can be continued with further boosts , and so on.
We can perform longitudinal boosts on the initial conditions in a straightforward way using light cone coordinates. Recall from section 2.1 that the color field of a single nucleus in covariant gauge is given by the solution to the Poisson equation, see eq. 2.11
[TABLE]
In light cone coordinates a longitudinal boost along of rapidity is given by
[TABLE]
In the boost invariant scenario, the current remains unchanged under these boosts (see eqs. (2.63) and (2.64)). However, this is not the case anymore in the 3+1D setup. For instance, in case of the single color sheet approximation the color current for a right-moving nucleus is given by (see eq. 4.9)
[TABLE]
where is the Gaussian profile function of width . Applying the boost, the new current reads
[TABLE]
In the last two lines we have absorbed the factor into a redefined thickness parameter and defined a new (boosted) longitudinal profile function with thickness . Boosting in the longitudinal direction is therefore equivalent to simply rescaling the thickness parameter of the nuclei. Analogously, we find for the thickness of a left-moving, boosted nucleus. The associated color field is obtained by solving the Poisson equation (4.52) for the boosted color current.
As a proof-of-concept example we repeat the simulations of the previous subsections for and . The results for the two pressure components are shown in fig. 4.15. Here, the effects of the longitudinal boost on the single nucleus fields are clearly visible. For decreasing values of , the width of the left-moving nucleus decreases, while the profile of the right-moving nucleus becomes broader. The transverse pressure , which is generated by the Glasma fields, is “dragged along” by the left-moving nucleus towards negative . Finally, the rapidity profiles of the LRF energy density at are shown in fig. 4.16. In this plot it is evident that the boosted simulations agree in the overlapping regions, which demonstrates the (approximate) Lorentz symmetry of the numerical method. By combining the results of different boosted collisions, it is possible to extend the available range of observables by one full unit of rapidity. This method is however still limited: by boosting to higher (absolute) rapidities, the thickness of one of the nuclei can become rather small compared to the longitudinal lattice spacing . For this reason, the lattice resolution was increased for the highest boost in fig. 4.15.
Chapter 5 A semi-implicit variational integrator
In the previous chapters we presented two slightly different numerical methods to solve the classical Yang-Mills equations in the real-time lattice formalism. Both the Hamiltonian leapfrog and Lagrangian leapfrog method of section 3.2.1 are explicit schemes and accurate up to quadratic error terms in the time step and spatial lattice spacings . In terms of accuracy, the two schemes can be considered equivalent. As already mentioned in section 4.2, we find that using these schemes, our simulations exhibit a numerical instability that results in the production of unphysical spurious fields and an artificial increase of field energy over time. Generally, increasing lattice resolution dampens this instability up to a point where we can trust our numerical results to be accurate enough at the cost of increased computational resources. Unfortunately, using the standard schemes, no choice of lattice parameters is able to completely eliminate the instability. In this chapter, we identify it as a non-Abelian version of the numerical Cherenkov instability (NCI) [123] known from Abelian particle-in-cell simulations. Having identified the source of the instability, we are able to develop a new semi-implicit scheme for real-time lattice gauge theory, which almost completely avoids the problems associated with the standard schemes while keeping the same order of accuracy and lattice gauge invariance. This chapter is heavily based on the most recent publication [4] with some additional comments and plots added.
5.1 Numerical Cherenkov instability and numerical dispersion
We start our discussion by mentioning that the numerical instability already occurs when studying the propagation of single nuclei on the lattice. The problem is therefore not primarily related to the actual collision of nuclei. In principle, the propagation of a single nucleus along the beam axis should be static and stable. In the continuum, both the current and the color fields move at the speed of light while maintaining their shape along the longitudinal coordinate. For observers in the laboratory frame, the transverse color structure of both fields and charges appears static due to time dilation. One encounters this property when solving the Yang-Mills equations for single nuclei (see section 2.1): the propagating solutions to the field equations must be independent of one of the light cone coordinates or (at least in covariant and LC gauge).
Ideally, we would like to observe the same properties when solving the Yang-Mills equations numerically on the lattice: both fields and the ensemble of colored particles should move at the speed of light along the beam axis in a stable manner. Let us therefore re-examine the numerical methods of chapter 3 more closely with regards to propagation at the speed of light. The colored point-like particles which make up the color charge of the nucleus, always move at the speed of light: since we neglect any forces acting on particles, their continuous positions are updated at each time step such that they follow lightlike trajectories. Therefore, in terms of movement, the particle degrees of freedom behave ideally by construction. However, the color charges of particles can be affected by parallel transport. As colored particles move along the beam axis, their color charges are rotated via Wilson lines in the lightlike direction of their trajectory. In temporal gauge, the gauge fields of nuclei are purely transverse and therefore the color charges should ideally remain static. Color charges only rotate if they encounter non-zero longitudinal gauge fields.
In contrast, the propagation of color fields is governed by the discretized field equations of motion on the lattice. The propagating color fields of nuclei can be decomposed into a superposition of plane wave modes of varying momenta and associated frequencies . In the continuum, plane wave solutions of the Yang-Mills equations have a linear dispersion relation and as such they all exhibit the same phase velocity equal to the speed of light. This is generally not the case on the lattice. Typically, finite difference schemes such as the Hamiltonian or Lagrangian leapfrog method (see section 3.2.1) show modified dispersion behavior for high momenta . This phenomenon is known as numerical dispersion.
For concreteness, we now derive the wave dispersion relation of the Lagrangian leapfrog method presented in section 3.2.1. Starting from eq. 3.49 without external sources
[TABLE]
we first restrict the color field to an Abelian subspace, because plane wave solutions to the Yang-Mills equations are effectively Abelian. As a consequence, all commutator terms can be dropped. Consequently, the plaquette (see eq. 2.101) can be written as
[TABLE]
where the lattice field strength tensor is given by
[TABLE]
Here, is simply the first order forward finite difference. In this notation, the lattice gauge field is considered to be defined at the mid-point between and . The EOM (5.1) are still non-linear due to the use of the matrix exponential in eq. 5.2. Therefore, we go to the limit of small amplitudes and find from eq. 5.1 to lowest order in the field strength tensor
[TABLE]
which is simply a finite difference approximation of Maxwell’s equations. Similarly, we find the discrete Abelian Gauss constraint
[TABLE]
The wave dispersion relation can now be derived from these EOM by plugging in a plane wave ansatz
[TABLE]
which yields an equation that connects the frequency to the momentum . At this point we simply state the result and postpone the detailed derivation to section 5.4.1:
[TABLE]
Restricting the momentum to one particular direction and solving for yields the non-linear dispersion relation
[TABLE]
where is the spatial lattice spacing and . For small momenta we recover the linear result , but in general, is a non-linear function of . In fig. 5.1 we plot the dispersion relation and its associated phase velocity for a few parameter choices. In particular, we see that the phase velocity , i.e. the speed at which a plane wave mode with momentum propagates, is less than the speed of light. This effect is more pronounced at high momenta. Although choosing (or ) would render the dispersion relation linear, this particular parameter choice causes unstable wave propagation: requiring that arbitrary momenta are associated with real-valued frequencies when solving eq. 5.7, yields the stability criterion (see sections 5.3.1 and 5.4.1 for a derivation)
[TABLE]
which forbids the time step being equal to any of the lattice spacings for more than one spatial dimension. Violating the above stability constraint leads to modes with imaginary frequencies that can grow exponentially.
Deviations of phase velocity from the speed of light imply that in principle it is impossible to have stable propagation of color fields using the schemes derived in section 3.2.1. Over time, high momentum modes “fall behind” low momentum modes. This alone would not cause an instability, but merely change the longitudinal shape of the nuclear field as the simulation progresses. A problem arises when we take into account the coupling to the external color currents, which are not affected by numerical dispersion. Consider the Abelian scenario of a highly contracted charge density moving at the speed of light: the longitudinal part of Maxwell’s equations then reads
[TABLE]
where and are the electric field and current in the longitudinal direction. As in the non-Abelian case, the electric and magnetic field are purely transverse. This implies that there is a precise cancellation of and the curl of the magnetic field. In the continuum we therefore have . Any mismatch between the current and fields leads to the production of longitudinal electric fields. In the discretized scenario, the non-linear dispersion relation of the electromagnetic fields causes this mismatch and leads to a production of spurious field modes. Moving charges seem to be radiating and over time the total field energy of the system increases. In this case electric charges are traveling “faster than the speed of light”, in the sense that the speed of light is defined by the speed of propagation of the electromagnetic field. This is reminiscent of the physical phenomenon of Cherenkov radiation, where relativistic charges move through a medium in which electromagnetic radiation propagates at slower speeds. However, in the case of discretized systems subluminal phase velocities are numerical artifacts. This unphysical effect is therefore known as the numerical Cherenkov radiation or instability (NCI) [123]. In simulations of Abelian plasmas with dynamical particles, the instability arises because the artificial Cherenkov radiation further affects charged particles through the Lorentz force, which leads to unphysical acceleration of particles.
In contrast, the colored particles of our simulations are not affected by forces, but they interact with the color field via parallel transport. In temporal gauge a non-zero longitudinal color-electric field generally implies that the longitudinal gauge component becomes non-zero as well. The color charges of the nucleus start rotating and the initial mismatch is amplified further. The instability we are facing is therefore due to non-Abelian interactions.
Now that the main culprit is identified, one can think about possible solutions. In the last chapter the instability could be avoided (or rather postponed) by using very large grid sizes and finer lattice spacing. This helped because smaller increases the largest available momentum via . In fig. 5.1 the occupied field modes move towards smaller momenta (compared to ), where the dispersion relation is more linear. As a result, it takes longer for the color fields of the nucleus to disperse and start becoming unstable. However, this does not actually eliminate the instability and the computational cost of our simulations increases quickly. The only way to avoid the instability altogether is to eliminate numerical dispersion in the discretized system. A typical example of Abelian PIC simulations being heavily affected by the NCI is the simulation of laser wakefield acceleration [124]. A number of possible solutions has been developed: for example, instead of using a finite difference scheme in both spatial and temporal directions (finite difference time domain or FDTD), one can solve the equations of motion in Fourier space using so-called “spectral solvers”. In this method the dispersion relation can be modified more directly compared to finite difference schemes. In combination with filtering high momentum modes, the NCI is suppressed successfully [125]. Another solution is to carefully modify interpolation schemes and finite difference stencils in FDTD approaches, see [126, 127] or [128]. In the case of laser wakefield simulations it is in principle not necessary to eliminate the NCI for all propagation directions, but only along a particular beam axis. A numerical scheme based on this idea is developed in [129]. However, a particularly simple and elegant solution to the problem is the use of semi-implicit schemes to repair the dispersion relation (i.e. making it linear) for one direction of propagation [130], such that no mismatch between charges and fields can arise. For our 3+1 dimensional setup a corrected dispersion relation along the beam axis would be sufficient to mitigate the NCI. The development of a semi-implicit scheme is therefore the main goal of this chapter.
5.2 Lattice gauge invariance
One particularly important property of both the Hamiltonian and the Lagrangian leapfrog method for numerically solving the Yang-Mills equations is that there exists an exactly conserved Gauss constraint (GC). Due to this constraint being preserved throughout the time evolution, we are able to derive a covariant conservation equation for the color current and the color charge density on the lattice. Solving the continuity equation for single colored particles as in section 3.2.2 allows us to determine the correct color current necessary for solving the field equations of motion coupled to external sources. Without a preserved GC we have no clear notion of charge conservation on the lattice. It is therefore important that our sought-after semi-implicit scheme also possesses an associated, exactly conserved GC.
Keeping this in mind, it is probably not a good idea to develop new schemes by making arbitrary changes to the EOM, as it is highly unlikely that one ends up with a scheme that preserves the original GC. One then would have to “guess” a new GC without even knowing if such an equation exists. It is clear that a more systematic approach is necessary. Recall that the Lagrangian leapfrog method of section 3.2.1 was derived by first discretizing the Yang-Mills action (yielding the Wilson gauge action) and then varying this action. Analogously to the continuum system, variation with respect to spatial gauge links led to the discrete EOM and through variation with respect to temporal gauge links we found the associated GC. One can then explicitly prove that the GC is identically satisfied, see section C.4. It is no coincidence that the constraint holds exactly, because the action already implies charge conservation through lattice gauge invariance. In the continuum there is a direct connection between gauge invariance and charge conservation through the Noether theorem. Remarkably, this also holds for the discretized case. Consider a general action formulated in terms of gauge links such that approximates the Yang-Mills in the continuum limit . Further assume that is invariant under arbitrary lattice gauge transformations, i.e. with
[TABLE]
Now consider “small” gauge transformations with assumed to be infinitesimal. A single gauge link transforms according to
[TABLE]
where the gauge-covariant forward finite difference is given by (see eq. 2.110)
[TABLE]
The above relation can be interpreted as an infinitesimal perturbation of the gauge links, i.e.
[TABLE]
with . This corresponds to a perturbation of gauge fields by . Applying this transformation to the action and expanding in therefore yields
[TABLE]
where the derivative acting on gauge links is defined via (see eq. C.43)
[TABLE]
The above derivative follows from the variation of gauge links, see section C.2. Since we assume to be lattice gauge invariant, the term linear in must be identically zero (as well as terms of higher order in ). This yields
[TABLE]
Using the adjoint representation of the gauge-covariant finite difference
[TABLE]
where the gauge link in the adjoint representation is given by , we can use summation by parts to find
[TABLE]
This relation must hold for arbitrary and therefore we have
[TABLE]
which is similar to a gauge covariant continuity equation for the “current” . In temporal gauge we have and can identically rewrite the above equation as
[TABLE]
Recalling from section 3.2.1 (see also section C.3) that and yield the GC and the discrete EOM respectively, we see that if the EOM are satisfied, i.e.
[TABLE]
then the GC must be conserved
[TABLE]
or simply
[TABLE]
This shows that independent of the exact form of , the constraint will be conserved under discrete time evolution if the action is invariant under lattice gauge transformations. In particular, the Lagrangian leapfrog EOM, which are derived from the invariant Wilson gauge action, have this property. In the case of the Hamiltonian leapfrog EOM, which are derived by first taking the partial continuum limit and then re-discretizing the time coordinate, the conservation of the GC is less obvious but still holds. Note however, that the above proof does not apply to the case of the Hamiltonian leapfrog EOM.
A reasonable new numerical scheme for solving the Yang-Mills equations should therefore be derived from a gauge invariant action . Numerical schemes derived from a variational principle are known as variational integrators [131]. The main advantage of this type of numerical scheme is that any symmetries of the discrete action are inherited by the numerical scheme and any constraints associated with particular symmetries have discrete analogues. In the case of lattice gauge invariance, this leads to gauge covariant discrete EOM and exact constraint conservation. Moreover, variational integrators also tend to have improved energy conservation properties and it is possible to formulate a conserved symplectic structure, i.e. there is a notion of conserved phase space volume. For a pedagogical introduction see chapter 5 of [132].
The goal of the rest of this chapter is to formulate a new discrete action in terms of gauge links, which is invariant under lattice gauge transformations and approximates the Yang-Mills action with (at least) the same accuracy as the Wilson gauge action. Upon variation, this new action should yield a semi-implicit numerical scheme, which allows for stable, dispersion-free propagation along a particular axis on the lattice. In order to gain an understanding of how to formulate such an action, we begin with a simpler model, namely a free scalar field on the lattice. Then we continue with a semi-implicit scheme for Abelian gauge fields on the lattice, where we have a notion of lattice gauge invariance and charge conservation. Finally, we will be able to formulate the non-Abelian generalization of the Abelian semi-implicit scheme. Slightly breaking the conventions of earlier chapters we will choose as the longitudinal (or beam axis) direction.
5.3 Scalar field on the lattice
The basic ideas behind the numerical scheme we are after can be most easily explained using a simple toy model, namely the two-dimensional wave equation. We start by giving a few definitions and then derive three different numerical schemes by discretizing the action of the system in different ways and using a discrete variational principle. We will see how the exact discretization of the action affects the properties of the numerical scheme and in particular how numerical dispersion can be eliminated.
We consider a real-valued scalar field in 2+1D with mostly minuses metric signature and set the speed of light to . The action is given by
[TABLE]
which upon demanding that the variation of the action vanishes
[TABLE]
yields the equations of motion (EOM)
[TABLE]
We use Latin indices to denote the spatial components, is a shorthand for (no sum implied) and . Inserting a plane-wave ansatz
[TABLE]
with and into the EOM gives the dispersion relation
[TABLE]
Obviously, the phase velocity is constant, i.e. there is no dispersion.
5.3.1 Leapfrog scheme
One possible way of discretizing the action is
[TABLE]
where is the sum over all lattice sites and is the space-time volume of a unit cell. Introducing small variations of the discrete field at each point, the discrete variation of this action reads
[TABLE]
which upon setting it to zero yields the discretized EOM
[TABLE]
where and are second order finite differences
[TABLE]
In the above derivation, we made use of summation by parts, i.e.
[TABLE]
which is the discrete analogue of integration by parts. If the field is known in two consecutive time slices we can explicitly solve the EOM (5.32) for the field values in the next time slice:
[TABLE]
In fact, this scheme is identical to the explicit leapfrog scheme111 The connection to the leapfrog scheme becomes more apparent if we introduce an approximation of the conjugate momentum
(5.36)
which is defined naturally between time slices and (hence the index “” in our notation). The EOM can then be written as
(5.37)
(5.38)
, which is accurate up to second order in the time step and spatial lattice spacings . Using the plane-wave ansatz (5.28) we find the dispersion relation
[TABLE]
which is in general non-linear and only yields real-valued (stable) frequencies for all wave vectors if the Courant-Friedrichs-Lewy (CFL) condition holds
[TABLE]
The discretization errors of this finite difference scheme result in a non-linear dispersion relation, which is usually referred to as numerical dispersion, since this kind of artificial dispersive behavior of plane waves does not show up in the continuum. If it were possible to set (the so-called “magic time-step”) the leapfrog scheme would actually be non-dispersive along the lattice axes, but this choice of the parameters is forbidden by the CFL condition in higher dimensions than and would lead to unstable modes as discussed in section 5.1.
5.3.2 Implicit scheme
Let us consider a different discretization: we define a new action
[TABLE]
where is the temporally averaged field
[TABLE]
Note that only one of the spatial finite differences in the squared term is temporally averaged. Since this action differs from the leapfrog action (5.30) only up to an error term quadratic in , the numerical scheme derived from this action will have the same accuracy as the leapfrog scheme.
Repeating the steps as before we obtain the discretized EOM
[TABLE]
This is an implicit scheme, which is more complicated to solve compared to the explicit leapfrog scheme given by eq. 5.35. Here we have to find the solution to a system of linear equations, which can be accomplished using (for instance) iterative methods.
The dispersion relation for this scheme reads
[TABLE]
This relation can always be solved for real-valued frequencies and therefore the implicit scheme is unconditionally stable. Unfortunately, this does not solve the problem of numerical dispersion either, because there is no choice of lattice parameters that results in a linear dispersion relation.
We quickly summarize: the first action we considered given by eq. 5.30 gave us the explicit leapfrog scheme, which is rendered non-dispersive but unstable using the “magic time-step”. The second action, eq. 5.41, which we obtained by replacing one of the spatial finite differences with a temporally averaged expression, yields an implicit scheme. This scheme is unconditionally stable, but always dispersive. This suggests that a mixture of both discretizations might solve our problem.
5.3.3 Semi-implicit scheme
Finally, we consider the action
[TABLE]
where the derivatives with respect to (w.r.t.) are treated like in the leapfrog scheme and the derivatives w.r.t. involve a temporally averaged expression as in the implicit scheme. The EOM now read
[TABLE]
We call this numerical scheme semi-implicit, because the finite difference equation contains both explicitly and implicitly treated spatial derivatives. The dispersion relation associated with eq. 5.46 is given by
[TABLE]
which is stable if
[TABLE]
The CFL condition (5.48) now allows us to set . Looking at the dispersion relation (5.47) we notice that for , but the propagation becomes non-dispersive, i.e. . For and the propagation still exhibits numerical dispersion. The scheme defined by the action (5.45) therefore allows for non-dispersive, stable wave propagation along one particular direction on the lattice. This principle also extends to systems with more spatial dimensions, where one treats a preferred direction explicitly and all other spatial directions implicitly.
A comparison of dispersion relations for all three numerical schemes is shown in fig. 5.2. Although the semi-implicit scheme is ideal for wave propagation along the direction, it performs worse along the orthogonal direction . The semi-implicit scheme therefore offers a trade-off if modified dispersion along other directions is acceptable for the physical situation one is modeling numerically.
5.3.4 Solution method and numerical tests
To solve the EOM of the implicit or the semi-implicit scheme, one has to solve a linear system of equations. This can be accomplished for instance by inverting a band matrix. Alternatively the equations can also be solved in an iterative manner. Taking the latter approach will be readily applicable to lattice gauge theory. One example for an iterative method is damped (or relaxed) fixed point iteration: the idea is to first rewrite the EOM (5.46) as a fixed point equation
[TABLE]
and then, starting with an initial guess from e.g. the explicit leapfrog evolution, use the iteration
[TABLE]
to obtain a new guess . Here the real-valued parameter acts as a damping coefficient. Using fixed point iteration might induce other numerical instabilities not covered by the CFL condition (5.48). To analyze this we make the ansatz
[TABLE]
where is the true solution to eq. 5.49 and represents a time-independent error term. The growth of the error is determined by the modulus of . Employing a Fourier ansatz yields
[TABLE]
which is independent of the component and the corresponding lattice spacing . Requiring convergence for high- modes, i.e. for , we find
[TABLE]
where
[TABLE]
In dimensions, where we treat the direction explicitly and all others implicitly, the stability condition is given by eq. 5.53 with
[TABLE]
Note that for damping is not necessary. A similar stability condition can be derived also for the implicit scheme, which by itself (just from the plane wave analysis) is unconditionally stable. It is important to keep in mind that the use of fixed point iteration can introduce new instabilities depending on the lattice spacing and the time step.
Finally, we perform a crucial numerical test: we compare the propagation of a Gaussian pulse using the three different schemes to show the effects of numerical dispersion and in particular that the semi-implicit scheme is dispersion-free. For simulations using the implicit or semi-implicit method, we solve the equations using damped fixed point iteration. The results are shown in fig. 5.3.
The main insight of this section is that the specific discretization of the action completely fixes the numerical scheme of the discrete equations of motion (and the associated stability and dispersion properties) which one obtains from a discrete variational principle. The use of temporally averaged quantities in the action leads to implicit schemes. If we treat some derivatives explicitly and some implicitly we can end up with a semi-implicit scheme that can be non-dispersive and still stable for propagation along a single direction on the lattice. As it turns out, this is just what we need to suppress the numerical Cherenkov instability we encountered in our 3+1 dimensional setup for heavy-ion collision simulations. In the next section we will see how we can use the same “trick” for Abelian gauge fields on the lattice.
5.4 Abelian gauge fields on the lattice
Before tackling the problem of non-Abelian gauge fields on the lattice, it is instructive to see how we can derive a dispersion-free semi-implicit scheme for discretized Abelian gauge fields. We will approach the problem as before: starting with a discretization of the action, we apply a discrete variational principle to derive discrete equations of motion and constraints. Then we will see what modifications to the action are required to obtain implicit and semi-implicit numerical schemes. Since we are dealing with gauge theory we will take care to retain gauge invariance also for the discretized system.
In the continuum the action for Abelian gauge fields reads
[TABLE]
with the field strength tensor given by
[TABLE]
The field strength and the action are invariant under gauge transformations defined by
[TABLE]
where is an at least twice differentiable function which defines the gauge transformation.
Varying the action with respect to the gauge field leads to
[TABLE]
The term proportional to leads to the GC
[TABLE]
while the term proportional to the variation of spatial components gives the EOM
[TABLE]
It is trivial to see that the EOM imply the conservation of the GC. As we will see next, it is also possible to formulate a discretization of the system where this conservation of the constraint holds exactly.
5.4.1 Leapfrog scheme
We consider a discretized gauge field at the lattice sites . The field strength tensor at is defined using forward differences
[TABLE]
which is antisymmetric in the Lorentz index pair like its continuum analogue. Furthermore, the lattice field strength is invariant under lattice gauge transformations given by
[TABLE]
where defines the local transformation at each lattice site . A straightforward discretization of the gauge field action is given by
[TABLE]
where is the space-time volume of a unit cell. Due to the invariance of under lattice gauge transformations the action is also gauge invariant.
Performing the variation of eq. 5.64 with respect to spatial components yields the EOM. Using the variation of the magnetic part of the action
[TABLE]
and the variation of the electric part w.r.t. spatial components (denoted by )
[TABLE]
we find the discrete EOM
[TABLE]
which are of the explicit leapfrog type. We also obtain a discretized version of the GC by considering the variation w.r.t. temporal components (denoted by ). With
[TABLE]
we get the constraint
[TABLE]
Since both the discrete EOM and the constraint follow from the same discretized, gauge-invariant action (5.64), the GC is guaranteed to be automatically conserved under the EOM. We can show this explicitly via
[TABLE]
which is equivalent to
[TABLE]
This means that if the GC is satisfied in one time slice then the EOM will ensure that it remains satisfied in the next time slice. We can also give a more general proof analogous to the one given in section 5.2: consider an infinitesimal gauge transformation
[TABLE]
and expand the action for small . We then find
[TABLE]
Since is invariant for any it must hold that
[TABLE]
or written slightly differently
[TABLE]
If the equations of motion are satisfied in every time slice, i.e. , then the GC must be conserved from one slice to the next:
[TABLE]
This holds regardless of the exact form of the gauge invariant action . Consequently, it does not matter what kind of discretization of the action we use. As long as retains lattice gauge invariance in the sense of eq. 5.63, we are guaranteed to find that the discrete GC is conserved under the discrete equations of motion.
The EOM given by eq. 5.67 alone are not enough to uniquely determine the time evolution of the field : we must specify a gauge condition. As before, we use temporal gauge
[TABLE]
which we also use in the case of non-Abelian lattice gauge fields. The EOM then read
[TABLE]
Using a plane wave ansatz
[TABLE]
with amplitude , we find the same non-trivial dispersion relation and CFL stability condition as in the case of the leapfrog scheme for the wave equation in 2+1D (see eqs. (5.39) and (5.40)). To show this explicitly we first introduce some notation. Taking a forward (backward) finite difference of the plane wave ansatz yields
[TABLE]
which suggests the definition of the forward (backward) lattice momentum
[TABLE]
It holds that . We define the squared lattice momentum as
[TABLE]
Furthermore, we can render the lattice momenta dimensionless by multiplying with . We define the dimensionless lattice momentum as
[TABLE]
and
[TABLE]
The factor of in eq. 5.85 is introduced for convenience. For differences with respect to the time coordinate we find . Using these definitions the GC (5.69) for the plane wave ansatz can be reduced to
[TABLE]
The discrete EOM (5.78) can be written as
[TABLE]
One can eliminate the mixed terms using the GC and find the dispersion relation
[TABLE]
which is equivalent to the dispersion relation of the wave equation (5.39), i.e.
[TABLE]
5.4.2 Implicit scheme
An implicit scheme analogous to the one we derived for the wave equation (see eq. 5.41) can be found using the action
[TABLE]
where we introduce the temporally averaged field-strength
[TABLE]
The averaged field strength differs from only by an error term proportional to . Replacing one of the field strengths in the quadratic term in the action with its temporally averaged expression is analogous to replacing the wave amplitude with in the action (5.41). The averaged field strength is also invariant under lattice gauge transformations:
[TABLE]
Varying the action as we did for the leapfrog scheme we find the EOM
[TABLE]
and employing temporal gauge we have
[TABLE]
where the fields have been replaced with temporally averaged expressions on the right-hand side. The GC that arises from varying w.r.t. temporal components is simply (5.69), because we did not modify the term involving or introduce new dependencies on . The discrete EOM (5.94) still conserve the GC due to lattice gauge invariance. Performing the same steps as in the previous section, we can show that this implicit scheme is unconditionally stable and exhibits the same non-trivial dispersion relation as the implicit scheme for the wave equation, see eq. 5.44.
5.4.3 Semi-implicit scheme
In this section we want to develop a semi-implicit scheme for Abelian gauge fields. Specifically we need an action that allows for dispersion-free propagation of waves in the direction of the coordinate (which we refer to as the longitudinal direction). We call and the transverse coordinates. In this section, Latin indices refer to transverse indices and the longitudinal index will always be explicit. Our goal is to define the action in such a way that we obtain equations of motion that include explicit differences in the direction, but temporally averaged finite differences in the and direction. This means that we have to modify the term of the leapfrog action (5.64) such that it results in terms like . A first guess for a semi-averaged version of field strength that could accomplish this is
[TABLE]
with . However, it turns out that such a term is not invariant under lattice gauge transformations (5.63). The problem is that transforms differently than . We have
[TABLE]
which yields
[TABLE]
where the last term breaks gauge invariance. To fix this we introduce the “properly” averaged field strength given by
[TABLE]
where the last term transforms as
[TABLE]
Here we made use of the exact relation
[TABLE]
Therefore, the transformation property of the properly averaged gauge field is the same as :
[TABLE]
It still holds that in the continuum limit the properly averaged gauge field is the same as up to an error term quadratic in . While the definition (5.100) seems a bit arbitrary at first sight, this way of averaging is more natural using the language of lattice gauge theory: in section 5.5.3 we will find an intuitive picture in terms of Wilson lines that reduces to eq. 5.100 in the Abelian limit for small lattice spacings.
The properly semi-averaged gauge-invariant field strength is then given by
[TABLE]
In order to keep the field strength explicitly antisymmetric we define . Using these definitions we can propose the action
[TABLE]
where the indices denote transverse components. Since we have built the new action from gauge invariant expressions it is also invariant under lattice gauge transformations. The use of in introduces new terms in the action (5.105) dependent on the temporal component of the gauge field. Although these terms disappear in temporal gauge (our preferred choice), they still have an effect on the scheme since we have to perform the variation before choosing a gauge. Therefore, we will obtain a modified GC compatible with the equations of motion derived from the action (5.105) even after setting .
At this point the question might arise if the “proper” averaging procedure has any effect on the implicit scheme of the previous section as well. It turns out that if one defines the averaged field-strength of eq. 5.92 using the properly averaged gauge field , the action of the implicit scheme (and by extension the EOM and the constraint) remains unchanged. This is due to the fact that the terms proportional to in eq. 5.100 cancel:
[TABLE]
Therefore, no such modification is required in the implicit scheme.
Varying this action w.r.t. temporal components yields the modified GC
[TABLE]
Since explicitly depends on , we obtain a correction term to the standard leapfrog GC (5.69). The discrete EOM read
[TABLE]
We have separate EOM for the longitudinal and transverse components of the gauge field. By replacing the averaged expressions and with , the EOM reduce to the leapfrog equations as expected.
The propagation of waves in the semi-implicit scheme turns out to be more complicated compared to the leapfrog or implicit scheme: given a wave vector and a field amplitude (such that the GC (5.107) is satisfied) the dispersion relation becomes polarization dependent, i.e. the scheme exhibits birefringence. The amplitude of an arbitrary plane wave
[TABLE]
has to be split into a longitudinal and two momentum-dependent transverse components
[TABLE]
where are dimensionless lattice momenta given by eq. 5.85. In appendix D we find that the components and have the dispersion relation
[TABLE]
and the component has a second different dispersion relation
[TABLE]
Analyzing the stability of the scheme using the two dispersion relations, yields that it is stable if
[TABLE]
which, when requiring stability for all modes, reduces to
[TABLE]
If we consider the special case of a plane wave with a purely transverse amplitude and which propagates in the direction (i.e. setting the transverse momenta ), we find that both dispersion relations agree:
[TABLE]
This dispersion becomes linear if we set . This explicitly shows that the semi-implicit scheme for Abelian gauge fields allows for dispersion-free, stable propagation along the longitudinal direction if we use the “magic time-step”. The dispersion relations also agree if we set the longitudinal momentum to zero for arbitrary transverse momenta. In general however, wave propagation in this scheme is birefractive. It would be interesting to see if there are alternative discretizations of the action which allow for dispersion-free propagation without being birefractive.
The main result of this section is the action (5.105) which gives rise to the semi-implicit scheme. Here we used a combination of differently averaged field strengths, and , in the action. Our next goal is to generalize these expressions to non-Abelian gauge fields.
5.5 Non-Abelian gauge fields on the lattice
Recall that the continuum action for non-Abelian Yang-Mills fields is given by
[TABLE]
where the field strength is
[TABLE]
Through variation of the action we obtain the GC and the EOM:
[TABLE]
where the gauge-covariant derivative acting on an algebra element is given by
[TABLE]
In real-time lattice gauge theory, instead of gauge fields , we use gauge links (or link variables) as already discussed previously.
The standard discretization using plaquettes leads to the Wilson gauge action (see section 2.3.2)
[TABLE]
where denotes a sum over all spatial components. Using
[TABLE]
it is clear that the action (5.5) is a discretization of the continuum Yang-Mills action
[TABLE]
where we made the split into temporal and spatial components explicit. Since it is built from gauge-invariant expressions, the Wilson gauge action (5.5) is invariant under lattice gauge transformations.
At this point we remark that the continuum Yang-Mills action (5.124) and its discretization (5.5) look very different: while the Yang-Mills action is given in terms of squares of the field strength tensor, the Wilson gauge action is linear in plaquette variables. In terms of plaquettes it is not immediately clear how we might generalize our approach from section 5.4. Fortunately, the action (5.5) can be written differently so that its functional form is more similar to eq. 5.124. We define (see for instance p. 94 of [62])
[TABLE]
which transforms non-locally
[TABLE]
For comparison to the plaquette , the path traced by is shown in fig. 5.4 on the right. In the continuum limit can be identified (up to constant factors) with the field strength : expanding for small lattice spacing we find
[TABLE]
Most noteworthy is the exact relation
[TABLE]
with which we can identically rewrite the action as
[TABLE]
This functional form of the rewritten action (5.129) is now virtually the same as the continuum case (5.124). We will exploit this when generalizing the implicit (5.91) and semi-implicit schemes (5.105) to non-Abelian gauge fields.
5.5.1 Leapfrog scheme
For comprehensiveness we summarize the derivation of the Lagrangian leapfrog method. Following the same procedure as in the case of Abelian gauge fields, we vary eq. 5.129 w.r.t. spatial components to obtain the discrete EOM and w.r.t. temporal components to find the GC.
Starting with the constraint we get (see section C.3.1 for details)
[TABLE]
where
[TABLE]
Varying the spatial link variables we obtain the EOM (see section C.3.2 for more details)
[TABLE]
which upon using the definition of can be written in the more familiar form
[TABLE]
As before, time evolution under the EOM conserves the GC exactly (see section C.4). In order to actually solve the equations we specify the temporal gauge
[TABLE]
which enables us to compute the spatial link variables of the next time-slice using past links and the temporal plaquette:
[TABLE]
The temporal plaquette has to be determined from eq. 5.133. For SU(2) this can be done explicitly (see section 3.2.1).
5.5.2 Implicit scheme
Guided by what we learned from the Abelian case in section 5.4, we would like to replace one of the expressions in the Wilson action (5.129) with a time-averaged equivalent. At the same time we need to retain the gauge invariance of the action. Simply using the temporally averaged expression
[TABLE]
is not enough because and transform differently. A solution is to include temporal gauge links in order to “pull back” and to the lattice site . This leads us to the definition of the “properly” averaged field strength
[TABLE]
which transforms like , i.e.
[TABLE]
This gauge-covariant averaging procedure can be generalized: consider an object that transforms like
[TABLE]
As an example, could be a Wilson line connecting points and along some arbitrary path. A time-averaged version of is given by
[TABLE]
where is simply shifted up (or down) by one time step. It still transforms like eq. 5.139, i.e.
[TABLE]
Using this we can write
[TABLE]
Condensing our notation even further we write
[TABLE]
and
[TABLE]
It still holds that
[TABLE]
so the use of (instead of ) does not change the accuracy of the scheme. Note that temporal gauge renders all temporal link variables trivial and eq. 5.140 reduces to the simple time average.
The inclusion of temporal links in breaks a symmetry on the lattice: the diagonal from to is now a preferred direction, which can be seen in fig. 5.5. The loss of symmetry can be mitigated by also including terms like in the action. Therefore we propose the action
[TABLE]
where we explicitly include terms with negative spatial indices using the sum over positive and negative components to keep the action as symmetric as possible. The action is also invariant under time reversal, real-valued (see section C.5 for a proof) and gauge invariant. While we have not made any changes to the terms involving temporal plaquettes, the spatial plaquette terms now include temporal links and therefore we will also obtain a modified GC like in the semi-implicit scheme for Abelian fields.
Performing the variation the same way we did for the leapfrog scheme, we obtain the GC (see section C.5.1)
[TABLE]
where denotes the sum over positive and negative indices . The left hand side (LHS) of eq. 5.148 is the same as in the leapfrog scheme (5.130), but now there is also a new term on the right hand side (RHS) from varying the spatial part of the action. Performing the continuum limit for the GC (after multiplying both sides with ), the RHS term vanishes as . This shows that the RHS is not a physical contribution, but rather an artifact of the implicit scheme. The correct continuum limit of the constraint (and the EOM) is already guaranteed by the action.
In a similar fashion as before, we perform the variation w.r.t. spatial link variables to get the discrete EOM. We find (see section C.5.2)
[TABLE]
which is formally similar to the leapfrog scheme (5.132) with in place of , a sum over positive and negative components (instead of just positive indices) and an additional factor of to avoid overcounting. As a simple check one can replace with in eq. 5.149 (only introducing an irrelevant error term quadratic in ) and recover eq. 5.132.
Compared to the leapfrog scheme, solving eq. 5.149 is more complicated: it is not possible to explicitly solve for the temporal plaquette anymore because on the RHS involves contributions from both “past” and “future” link variables. Completely analogous to the case of Abelian gauge fields on a lattice, we obtain an implicit scheme by introducing time-averaged field strength terms in the action.
Similar to what we did in section 5.3.3, we propose to solve eq. 5.149 iteratively using damped fixed point iteration. Starting from an initial guess for the future link variable , for instance by performing a single evolution step using the leapfrog scheme, we iterate
[TABLE]
using in from the last iteration step to determine . In the above equations has been replaced with the unknown variable . We first solve eq. 5.150 for and then update using eq. 5.151. The parameter is used as a damping coefficient to mitigate numerical instabilities induced by fixed point iteration. For SU(2) we construct the temporal plaquette from using eqs. (3.60) and (3.61) and using temporal gauge we update the link variables via
[TABLE]
Then we repeat the iteration until convergence.
This iteration scheme can be used to solve the EOM (5.149) until the GC (5.148) is satisfied up to the desired numerical accuracy. Conversely, this means that unlike the leapfrog scheme, where the GC (5.130) is always satisfied up to machine precision in a single evolution step, the implicit scheme, solved via an iterative scheme, only approximately conserves the GC eq. 5.148. However, in section 5.6 we will show that using a high number of iterations the constraint can be indeed fulfilled to arbitrary accuracy. In practice however, we find that a lower number of iterations is sufficient for stable and acceptably accurate simulations at the cost of small violations of the constraint.
It is also immediately obvious that solving the implicit scheme requires higher computational effort compared to the leapfrog scheme. Considering that one has to use the leapfrog scheme for a single evolution step once (as an initial guess) and then use fixed point iteration, where every step is at least as computationally demanding as a single leapfrog step, it becomes clear that the use of an implicit scheme is only viable if increased stability allows the use of coarser lattices while maintaining accurate results.
5.5.3 Semi-implicit scheme
Using our knowledge from sections 5.4 and 5.5.2 we can now generalize the semi-implicit scheme to real-time lattice gauge theory. An appropriate generalization of the semi-averaged field strength (5.104) is given by
[TABLE]
where
[TABLE]
We also define . Using the time-averaging notation (see eq. 5.140) this can be written more compactly as
[TABLE]
where
[TABLE]
Note that
[TABLE]
which shows that the semi-averaged field strength only differs from by an irrelevant error term. Taking the Abelian limit (i.e. neglecting commutator terms) of eq. 5.157 and expanding for small lattice spacing yields
[TABLE]
We find that the linear term of the gauge-covariant average (5.157) agrees with the expression we constructed in the Abelian semi-implicit scheme (5.100). While we had to include an “arbitrary” correction term in the Abelian scheme to fix gauge invariance, the link formalism of lattice gauge theory forces us to only consider closed paths constructed from gauge links in the action and thus naturally leads us to the “proper” averaging procedure.
The Wilson line path traced by , as compared to , is shown in fig. 5.5 on the right. As before, we keep the action as symmetric as possible by also including terms with negative transverse directions, i.e. . Inspired by the Abelian semi-implicit case (5.105), we define the new action as
[TABLE]
where the sum over and only run over transverse coordinates and is the longitudinal coordinate. The purely transverse part of the action uses the same terms as the implicit scheme, see eq. 5.147. The longitudinal-transverse part is now given in terms of and analogous to eq. 5.91. We have to explicitly include the hermitian conjugate in order to keep the action real-valued.
Varying with respect to temporal components yields the GC (see section C.6.1)
[TABLE]
where we use the shorthand
[TABLE]
Varying w.r.t. spatial links the discrete semi-implicit EOM read (see section C.6.2)
[TABLE]
and for transverse components
[TABLE]
where simply means summing over the terms with positive and negative longitudinal directions. We now have two sets of equations: one for longitudinal components, and one for the two transverse components.
The equations of motion can be written more compactly by introducing the symbol
[TABLE]
where can be exchanged for corresponding expressions with or and can be exchanged for its temporally averaged version . The longitudinal component of the EOM then reads
[TABLE]
and the transverse components are given by
[TABLE]
These equations can be solved numerically using damped fixed point iteration completely analogously to section 5.5.2. First, one obtains an initial guess for “future” link variables from a single leapfrog evolution step using eqs. (5.133) and (5.135). Then one iterates from until convergence:
Compute the next iteration using damped fixed point iteration: in eqs. (5.5.3) and (5.5.3) replace and , solve for the unknown ’s and update the temporal plaquettes using
[TABLE]
and analogously for and , where is the damping coefficient. 2. 2.
For SU(2) we can reconstruct the full temporal plaquette from its components with the identity (see eq. C.30)
[TABLE]
for and . 3. 3.
Using and , compute the spatial links and via
[TABLE] 4. 4.
Repeat with .
As with the implicit scheme, our approach to solving the equations in the semi-implicit scheme is an iterative one. The GC (5.5.3) is only approximately satisfied, depending on the degree of convergence222This is also true for the Abelian semi-implicit scheme. If the equations of motion are solved only approximately using an iterative method, then the conservation of the GC is also only approximate depending on the degree of convergence..
5.5.4 Coupling to external color currents
Up until now we have only considered pure Yang-Mills fields. In the continuum we can include external color currents by adding a term to the action.
[TABLE]
The EOM then read
[TABLE]
and due to gauge-covariant conservation of charge we have
[TABLE]
which is the non-Abelian continuity equation. On the lattice we can simply add a discrete term to the action as well:
[TABLE]
where includes a factor of (“lattice units”). We also made the split into 3+1 dimensions explicit using and . The variation of simply reads
[TABLE]
which gives the appropriate contributions to the GC and the EOM. In the leapfrog scheme the GC now reads
[TABLE]
and the EOM read
[TABLE]
The constraint taken together with the EOM imply the local conservation of charge (see section C.4)
[TABLE]
which is the discrete version of the continuity equation. For the implicit and semi-implicit schemes the procedure is the same: including the term simply leads to the appearance of on the RHS of the GC (see eqs. (5.148) and (5.5.3)) and on the RHS of the EOM (see eqs. (5.149) and (5.5.3) – (5.5.3)). Due to conservation of the GC without external charges, the continuity equation for the implicit and semi-implicit scheme is simply eq. 5.178 as well. This implies that our treatment of the external currents in terms of parallel transport does not require any modifications when using the newly derived schemes.
5.6 Numerical tests
In this last section we test the semi-implicit scheme on the propagation of a single nucleus in the CGC framework. For an observer at rest in the laboratory frame, the nucleus moves at the speed of light and consequently exhibits large time dilation. As the nucleus propagates, the interactions inside appear to be frozen and the field configuration is essentially static. On the lattice we would like to reproduce this behavior as well, but depending on the lattice resolution we run into the numerical Cherenkov instability, which leads to an artificial increase of the total field energy of the system.
As described in section 5.1, the root cause of the instability is numerical dispersion: in the CGC framework a nucleus consists of both propagating field modes and a longitudinal current generating the field around it. It is essentially a non-Abelian generalization of the field of a highly relativistic electric charge. In our simulation the color current is modeled as an ensemble of colored point-like particles moving at the speed of light along the beam axis. The current is unaffected by dispersion, i.e. it retains its shape perfectly as it propagates. The field modes suffer from numerical dispersion, which over time leads to a deformation of the original longitudinal profile. The mismatch between the color current and the field leads to the creation of spurious field modes, which interact with the color current non-linearly through parallel transport (color rotation) of the current. This increases the mismatch further and more spurious fields are created. As the simulation progresses, this eventually leads to a large artificial increase of total field energy. The effects of the instability can be quite dramatic as seen in fig. 5.6. The main difference to the numerical Cherenkov instability in Abelian PIC simulations is that in electromagnetic simulations the spurious field modes interact with the particles through the Lorentz force [123]. In our simulations we do not consider any acceleration of the particles (i.e. their trajectories are fixed), but interaction is still possible due to non-Abelian charge conservation (5.178), which requires rotating the color charge of the color current. Therefore our type of numerical Cherenkov instability is due to non-Abelian effects.
We now demonstrate that the instability is cured (or at least highly suppressed for all practical purposes) by using the semi-implicit scheme. We test the scheme’s ability to improve energy conservation in the following way: we place a single gold nucleus described by the McLerran-Venugopalan model in a simulation box of volume and set the longitudinal extent of the nucleus to roughly correspond to a boosted nucleus with Lorentz factor . The values for , and are determined from as detailed in section 4.4. After setting up the initial condition we let the nucleus freely propagate along the longitudinal axis. As the simulation runs we record the total field energy
[TABLE]
where and are the color-electric and -magnetic fields at each time step. On the lattice the electric and magnetic fields are approximated using plaquettes (see eqs. (3.4), (3.148) and (3.149)):
[TABLE]
We compute the relative change , which we plot as a function of time . In the continuum we would have , but due to numerical artifacts and the Cherenkov instability this is not the case in our simulations.
The numerical results are shown in fig. 5.7. We see that the leapfrog scheme leads to an exponential increase of the total energy over time, which can be suppressed using finer lattices. On the other hand, the semi-implicit scheme leads to better energy conservation even on a rather coarse lattice. Therefore, the resolution that is usually required to obtain accurate, stable results is lowered by using the semi-implicit scheme. However, using the new scheme might not always be economical: finer lattices suppress the instability as well and since the leapfrog scheme is computationally cheaper than the semi-implicit scheme, the leapfrog can be favorable in practice. On our test system (a single GB node on the VSC 3 cluster) the simulation using the semi-implicit scheme (SI 256) takes hours to finish, while the same simulation using the leapfrog with (LF 1024) takes roughly hours and with (LF 2048) hours. Even though energy conservation is not as good as SI 256, the longitudinal resolution is much better in comparison, enabling us to extract observables with higher accuracy. It should be noted however that our implementation of the leapfrog scheme is already highly optimized, while the implementation of the semi-implicit scheme is very basic and should be considered as a proof of concept. Further optimizations and simplifications of the semi-implicit scheme might make it the better choice in many cases.
As a second test we look at the violation of the GC. The leapfrog scheme (5.133) conserves its associated GC (5.130) identically, even for finite time-steps . In numerical simulations this conservation is not exact due to floating point number errors, but the violation is zero up to machine precision. On the other hand, the semi-implicit scheme has to be solved iteratively and therefore the results depend on the number of iterations used in the fixed point iteration method. We define the relative violation of the GC as the ratio of the absolute (squared) GC violation to the total (squared) charge on the lattice. In the case of the leapfrog scheme this reads
[TABLE]
where the sum runs over the spatial lattice of a single time slice at . The numerator depends on the GC of the scheme and has to be adjusted according to the implicit or semi-implicit method (either eq. 5.148 or (5.5.3) including the charge density on the RHS as discussed in section 5.5.4). In fig. 5.8 we show how the GC violation converges systematically towards zero as we increase the number of iterations. Therefore, even though we can not use an arbitrarily high number of iterations due to limited computational resources, the semi-implicit scheme conserves the GC in principle. The same holds for the purely implicit scheme. In practice it is not necessary to satisfy the constraint up to high precision, as observables such as the energy density seem to converge much faster to satisfying accuracy.
Chapter 6 Closing remarks
In this last chapter I present a short summery of my thesis and outline some extensions and generalizations that could be attempted in the future.
6.1 Summary
Starting with chapter 2, I recapitulated the theory and numerical methods for the boost invariant case, where the color currents of the colliding nuclei are assumed to be infinitesimally thin, i.e.
[TABLE]
see eq. 2.2. The particular structure of the color currents allows one to unambiguously divide space-time along the light cone centered around the collision event into four separate regions. Along the boundary of the future light cone ( and ), a set of matching conditions can be derived, which turn out to be the initial values for the Glasma fields at proper time (see eqs. (2.72) and (2.73)). Due to invariance of eq. 6.1 under longitudinal boosts, the Glasma fields can be chosen to be independent of space-time rapidity , reducing the system to effectively 2+1 dimensions. Moreover, the color current does not contribute for , since its support is limited to the boundary of the light cone. Therefore, in order to study the evolving Glasma, one has to solve the 2+1 dimensional, non-linear, homogeneous Yang-Mills equations (eqs. (2.84) and (2.85)). Although the initial values for the fields in the future light cone are known exactly, the full analytical solution for is not known in general.
Fortunately, approximate solutions can be found numerically using the methods of real-time lattice gauge theory. One of the main features of the lattice formulation of gauge field theories is that the local fields are replaced by non-local gauge variables or links , which are Wilson lines connecting nearest neighbors on the underlying rectangular lattice. This reformulation of Yang-Mills theory allows us to derive discretized actions (the most simple one being the Wilson action, see section 2.3.2), which are exactly invariant under a subset of the full gauge symmetry, namely lattice gauge transformations. The most important symmetry characterizing Yang-Mills theories is therefore also realized when computing numerical approximations of its solutions.
As a model for the Glasma initial conditions we employed the McLerran-Venugopalan (MV) model. The MV model approximates relativistic nuclei as longitudinally infinitesimally thin, transversally infinitely extended walls of random, homogeneously distributed color charge. Neither a finite transverse radius nor nucleonic structure are described by the original MV model, although in principle such features can be added as well.
Equipped with these tools, numerical simulations were performed to study the boost invariant Glasma. In particular we studied the components of the energy-momentum tensor of the Glasma, where a number of observations can be made (see figs. 2.4, 2.5 and 2.7 in section 2.4):
- •
Initially, the Glasma consists of purely longitudinal color-electric and color-magnetic fields in the form of color flux tubes. These fields are correlated over lengths inversely proportional to the saturation momentum , which also roughly fixes the radius of the flux tubes. This initial state is highly anisotropic. The longitudinal pressure is negative with and the transverse pressure positive with , where is the energy density of the Glasma. The pressure anisotropy is maximal in the sense of and .
- •
As the system evolves the flux tubes start to expand in the transverse plane and transverse (color)-electric and -magnetic fields are generated. This leads to an increase of and a decrease of . Nevertheless, the system never becomes isotropic () as the change in pressures stops when the longitudinal pressure reaches zero. The system ends up in an anisotropic state where and .
- •
This anisotropic late time behavior is reached after and is referred to as free streaming. In this case the energy density of the system exhibits characteristic behavior given by .
- •
Obviously, by construction the system and all of its observables are independent of space-time rapidity .
Although the 2+1D framework can be extended to include the rapidity coordinate in a straightforward way, it is less clear how to formulate rapidity dependent initial conditions that satisfy the Gauss constraint (GC). Previous 3+1 dimensional approaches therefore focus on studying non-Abelian plasma instabilities due to small fluctuations on top of the boost invariant Glasma background field [35, 36, 39, 41, 42, 84, 85, 86]. If one is interested in a generalization of eq. 6.1, i.e. color currents of the form
[TABLE]
then the above outlined approach is too constraining. Due to a small, but finite longitudinal support of the color charge density of the nuclei, one cannot rely on a well-separated future light cone anymore. Therefore, in chapter 3 a numerical method is developed to simulate this generalized case. The most obvious difference to the 2+1 dimensional method is that in order to accommodate finite width currents (6.2), the whole Minkowski space-time must be considered and color currents have to be explicitly taken into account when solving the Yang-Mills equations on the lattice. This generalization implies a number of differences:
- •
The Yang-Mills equations are written in terms of flat coordinates using the laboratory frame time and the longitudinal coordinate instead of formulating them in the curved coordinate patch spanned by and . As a result, the field equations simplify due to lack of explicit -dependent terms. The time evolution in this frame is an evolution in instead of . Therefore, the natural choice of gauge for this case is the laboratory frame temporal gauge . The standard leapfrog equations of motion in the lattice formalism are derived in section 3.2.1 for this purpose.
- •
Instead of specifying the Glasma fields at , one has to set the fields of single nuclei at some initial time before the collision has happened. Analytically one can find the solution to the Yang-Mills equations, if at time both nuclei are spatially separated such that their color fields do not interact and each single nucleus can be treated independently. Furthermore, the initial conditions have to be compatible with the choice of gauge. The temporal gauge initial conditions are derived in section 3.1 for the continuum case where it is shown that the color field of a single nucleus in temporal gauge is identical to its light cone gauge version. The lattice discretization of these initial conditions is derived in section 3.3.2.
- •
The 3+1D method takes finite width color currents explicitly into account. In order to describe these currents numerically, one can make use of the colored particle-in-cell method originating from non-Abelian plasma simulations. The main idea is to replace the smooth color charge density of a nucleus with a large number of densely placed colored point particles. The movement of this ensemble of particles approximates the smooth color current of the nucleus. Although each particle’s trajectory is known in advance as all forces acting on them are neglected (the eikonal approximation), their color charge precesses due to being parallel transported in the non-Abelian color field. This color charge precession or rotation must be accounted for in order to satisfy the Gauss constraint. All the necessary equations to perform the numerical time evolution of fields and currents in this manner have been derived in section 3.2.
In chapter 4 the previously developed method is applied to study collisions of nuclei with color currents of the form
[TABLE]
where is the longitudinal profile function of the nucleus and determines the distribution of color charge in the transverse plane. Equation 6.3 can be viewed as a regularization of eq. 6.1 and a special case of the more general form eq. 6.2. The color charges are correlated along the longitudinal extent of the nucleus, i.e. the nucleus only consists of a single sheet of color charge with some finite width. This specific case was therefore referred to as the single color sheet approximation. The width of the profile function is set by the Lorentz-contracted longitudinal extent of the nucleus, i.e. , where is the nuclear radius and is the Lorentz factor.
First, the 3+1D method was put to the test by directly comparing to the boost invariant case in section 4.2. Immediately after the collision, the longitudinal color-electric field in the transverse plane was compared to the Glasma initial conditions. It was demonstrated that when the profile function approaches the delta distribution , the results from the 3+1D method and the boost invariant initial conditions agree, see fig. 4.3. Next, the pressure components were studied more closely. The 3+1D Glasma behaves very similar to the boost invariant case: the system exhibits free streaming and the pressure anisotropy remains large, see figs. 4.5 and 4.6. However, a subtle difference which strongly depends on the longitudinal width of the colliding nuclei, can be observed: in the boost invariant case the magnetic and electric components roughly contribute equally to the energy density of the Glasma. This does not appear to be the case for three-dimensional collisions. Immediately after the collision event the longitudinal magnetic field appears to be reduced compared to the electric component. This mismatch grows with larger widths and larger values of the infrared regulator and is only slowly cured over time, see figs. 4.7 and 4.8. As a further test of the validity of the three-dimensional simulations, local energy conservation was explicitly checked using the non-Abelian version of Poynting’s theorem, see fig. 4.9.
Having established that the 3+1D method is able to produce valid results, initial conditions were chosen that resemble the collisions at RHIC more closely. This was done to demonstrate the main feature of the method, namely the explicit breaking of boost invariance. In section 4.4 we extracted the rapidity profile of the energy density of the Glasma in its rest frame. For certain choices of parameters, the obtained profiles are roughly in agreement with experimental results at RHIC, see figs. 4.12 and 4.14. Surprisingly, infrared regulation seems to have a large effect on the widths of these profiles. It can also be shown that the longitudinal velocity of the Glasma is very close to the free streaming result, see fig. 4.13. These findings indicate that the rapidity profile of the energy density develops very shortly after the collision. Afterwards, the Glasma expands in a free streaming manner. Although the range of rapidity that can be probed with sufficient numerical accuracy is rather limited, a workaround method has been developed to extend this interval slightly (see section 4.4.3).
Unfortunately, the original numerical approach presented in chapters 3 and 4, which is based on the standard Wilson action, suffers from a non-Abelian numerical Cherenkov instability. This instability leads to an artificial increase of field energy over time even when only considering the propagation of a single nucleus. In chapter 5 the main culprit was identified as the different numerical dispersion behaviors of the two types of degrees of freedom in the simulation: while the color current, composed of many colored point particles, keeps its shape as it moves, the discretized gauge field suffers from numerical dispersion. As a result, the growing mismatch between the color fields of nuclei and their color currents, drives an instability that is only cured when taking the continuum limit. To resolve this issue, an alternative lattice discretization of the Yang-Mills action was proposed, which eliminates the problematic dispersion behavior along the trajectories of the nuclei. It was shown that this action is invariant under lattice gauge transformations and exhibits the same order of accuracy as the original Wilson action. Varying this action led to semi-implicit equations of motion for the gauge field, which can still be solved numerically, see section 5.5.3. Compared to the leapfrog method, the semi-implicit method requires higher computational effort, but was shown to successfully suppress the numerical Cherenkov instability (see section 5.6).
6.2 Outlook
A straightforward continuation of the work presented in this thesis is the study of random fluctuations within the longitudinal extent of nuclei and their effects on the Glasma. Most of the focus in the CGC approach to the initial state in heavy-ion collisions has been on realistically modeling the transverse structure of nuclei. For instance, the IP-Glasma model includes impact parameter dependence for non-central collisions and nucleonic and sub-nucleonic fluctuations [23]. On the other hand, the random longitudinal structure of nuclei has been largely neglected, except for the approximation of the path ordered Wilson line as an infinitely thin stack of uncorrelated sheets of color charge [49] (see also section 2.1). The 3+1D setup using the numerical improvements of chapter 5 should make the study of initial conditions including these fluctuations feasible.
In the spirit of [49], the starting point of this study could be a generalization of the MV model including random fluctuations in the longitudinal extent. Extending eq. 2.36, the general charge density two-point function of this model (in covariant gauge, for a right-moving nucleus) can be expressed as
[TABLE]
where and refer to transverse and light cone coordinates, is the coupling constant, and is a generalization of the phenomenological MV model parameter. In this model, is considered to be a function of and . The first argument defines the overall longitudinal shape of the nucleus, while the second argument controls the correlation of random color charges along the longitudinal extent. The one-point function is assumed to be identically zero to implement color neutrality on average. A specific realization of eq. 6.4 is given by (see [49])
[TABLE]
where the functions and represent Gaussian profile functions with widths and respectively and is a constant. The parameter defines the longitudinal width (as in the single color sheet approximation, see chapter 4) and is a new parameter to control the amount of random longitudinal fluctuations in the nucleus. The two-point function can be interpreted as a regularization of eq. 2.36.
The Wilson line associated with the transformation to LC gauge is given by (see eq. 2.20)
[TABLE]
with as the infrared regulator to control the divergence in the model. For we obtain the MV model with finite width in the single color sheet approximation which has been thoroughly investigated in chapter 4. The boost invariant limit is achieved by taking the limit and with the order of limits being . The color field of the nucleus is then treated as a shock wave and the dependence of becomes trivial as detailed in section 2.1. Random longitudinal structure is accounted for using the color sheet regularization, see eqs. (2.38), (2.39) and (2.40).
For the present purpose it would be interesting to study a different limit: in the case of while keeping , we have and the generalized model reduces to the three-dimensional MV model with finite longitudinal width , see eq. 2.36. In this case the nucleus exhibits random longitudinal structure and it is this limit which has not been studied previously.
In principle, the model for initial conditions described above is simple to implement and simulate with the 3+1D method presented in this thesis. The procedure for setting up initial conditions discussed in section 3.3 is general enough to work with an arbitrary charge density regardless of the exact dependence on . It would then be interesting to study the Glasma created from collisions of such nuclei, in particular with regards to the breaking of boost invariance and the problem of pressure isotropization and (physical) instabilities. The effects of random variation in the longitudinal extent should be compared to the setting that has been thoroughly investigated in the past decade, namely rapidity-dependent fluctuations on top of a boost invariant background field [36, 41, 42, 54, 86]. The main difference to previous works is that the rapidity-dependent fluctuations which seed the Glasma instability are of different origin than usually assumed. In the proposed approach, deviations from boost invariance are due to the color structure of the colliding nuclei and not of quantum origin as e.g. in [86]. Furthermore, instead of having to model the Glasma by superimposing fluctuations and background fields, the Glasma field in the 3+1D method arises directly from solving the Yang-Mills equations. By construction, these fields fulfill the Gauss constraint. The 3+1D setting also differs from previous approaches in more subtle ways: in the “background plus fluctuations” setup, which is generally formulated in the future light cone using coordinates, the axis along the rapidity coordinate is assumed to have periodic boundary conditions. In the 3+1D setup with coordinates, the boundary condition in is more similar to an open boundary.
One can even conceive a slightly artificial, but more direct way of showing how much the variation in the longitudinal extent affects the Glasma: starting with nuclei in the single color sheet approximation, one could make arbitrary changes to the Wilson line inside the nucleus without affecting the asymptotic value
[TABLE]
In the boost invariant scenario, these different initial conditions would be considered identical, as only the asymptotic Wilson line is used as input for the initial Glasma fields. In contrast, using the 3+1D setup it should be possible to differentiate between different initial conditions which share the same asymptotic Wilson line, if random variations within the longitudinal extent strongly affect the produced Glasma. Therefore, a more thorough study of these effects and more direct comparisons to the boost invariant case could prove interesting and provide new insights into the creation and the dynamics of the Glasma with broken boost invariance.
Another phenomenon that calls for further investigation using 3+1D simulations is the production of axial charge in the earliest stages of the collision, which is related to the chiral magnetic effect [133, 134]. The chiral magnetic effect predicts an anomalous separation of charges in non-central collisions where large magnetic fields are present. This anomalous transport phenomenon is due to a chirality imbalance or axial charge, i.e. a difference in the numbers of left-handed and right-handed fermions in the system. Both the theoretical aspects and the experimental signatures of the chiral magnetic effect are so far inconclusive and active topics of current research. The importance of the chiral magnetic effect can be seen by the fact that in 2018 a new kind of experiment within Run-18 at RHIC has been devoted exclusively to the exploration of this effect: by performing collisions of nuclei pairs and comparing to isobars, which are nuclei with the same mass but different number of protons, one can perform controlled experiments where the only changing parameter is the electric charge and thus a 10% change in the magnetic field induced by its current [134].
A theoretical uncertainty in the chiral magnetic effect is how large the produced axial charge is and how it is created. Locally the generation of axial charge is due to the non-zero covariant divergence of the axial charge current which exhibits the axial anomaly
[TABLE]
Here refers to the number of fermion flavors and is the dual to the non-Abelian field strength tensor . Because of this anomaly the axial charge density is not conserved, but arises due to parallel (color-)electric and magnetic fields. The right hand side is readily computable in Yang-Mills simulations and has been investigated in boost invariant, Glasma-like scenarios using numerical and analytical methods [135, 136, 137]. While the boost invariant scenario allows for strong local axial charge production due to the nature of longitudinal (and therefore parallel) color-electric and magnetic flux tubes early in the evolution of the Glasma, there are topological arguments that disallow a net change in the total axial charge [20, 135]. However, this no-go theorem only strictly applies to the boost invariant case where the system is effectively 2+1 dimensional, so it might be possible to find different behavior in our 3+1 dimensional setup. Even if no net axial charge will be found, one could still study local variations of axial charge density, make direct comparisons to boost invariant simulations and investigate how the mismatch of color-electric and color-magnetic fields (see section 4.3.2) affects axial charge production.
An important quantity that has not been computed yet in the 3+1D setup is the gluon occupation number. In the boost invariant case this quantity can be readily computed as mentioned in section 2.4. It has been used to study the gluon spectrum in the Glasma [138], to fix phenomenological parameters of models (such as in the MV model) and to determine the average transverse momentum of gluons [32, 68]. Furthermore it can be used to directly compute flow coefficients from the initial state from a Fourier decomposition with respect to the azimuthal angle [81, 139].
There are two technical problems that complicate the extraction of gluon numbers in 3+1 dimensions. The first problem is related to gauge fixing: in the boost invariant case the common procedure is to use the transverse Coulomb gauge condition
[TABLE]
at some proper time . In the 3+1D method, which is formulated in the laboratory frame of the collision, the analogue to this gauge condition would be the three-dimensional Coulomb gauge
[TABLE]
fixed at some time . Equation 6.10 does not necessarily imply eq. 6.9 and therefore it is possible that the gluon numbers extracted in 3+1 dimensions differ from the definition used in the boost invariant scenario. It might be necessary to choose a different gauge condition than the obvious choice of eq. 6.10.
A second problem is related to the choice of frame: as the frame only covers the future light cone, one only obtains the gluon numbers of the Glasma. In the 3+1D method, both the Glasma and the color fields of the nuclei are part of the simulation. Therefore, unwanted contributions from the nuclei would be present in the gluon occupation numbers. Direct comparisons between the boost invariant and the 3+1D method are only possible if these contributions can be somehow canceled or filtered out.
Once these issues are solved, it would be interesting to revisit the study of rapidity profiles from section 4.4. Instead of using the rest-frame energy density, one could use the gluon occupation number as a function of momentum rapidity, which would allow for better comparison to experimental data.
Finally, a much harder problem to tackle would be to study the production of quarks in the background of fully non-perturbative Glasma fields in 3+1D. The main idea is to solve the Dirac equation coupled to an external non-Abelian gauge field, which is treated as a background field (i.e. no backreaction from the fermion fields to the gauge fields) and then compute the number of quark-antiquark pairs that arises from the interaction using the solution of the Dirac equation. This problem can be solved analytically for Abelian fields [140] and in the special case of proton-nucleus collisions [141, 142], where one can consider the color field of the proton to be weak compared to the color field of the nucleus.
The fully non-perturbative problem for nucleus-nucleus collisions can only be solved numerically and turns out to be quite complicated: early numerical investigations have been performed first in a 1+1 dimensional toy model [143] and later in boost invariant Glasma background fields [144]. However, there are multiple numerical difficulties arising in these earlier approaches: it turns out that even though the color fields are boost invariant (effectively 2+1 dimensional) one has to solve the Dirac equation in 3+1 dimensions. A further complication is the choice of coordinate systems. The boost invariant Glasma fields are naturally solved in co-moving coordinates in the future light cone, but the solution of the Dirac equation also depends on information about the field in other parts of the light cone. Consequently due to mixing of different reference frames, one has to transform the Glasma fields into a more suitable coordinate system and special care has to be taken with coordinate singularities associated with the co-moving frame. Although there has been progress in improving the mathematical formalism [145], the fully numerical calculation with Glasma fields in 3+1D could still provide new results.
In the light of these technical issues, the 3+1D approach using laboratory frame (Cartesian) coordinates seems exceptionally well-suited to revisit the problem of quark pair production in heavy-ion collisions: the simulations presented in this thesis already provide information about the color fields in all parts of the light cone and by construction there are no ambiguities with coordinate singularities. Thus, even though it would be a difficult task to undertake due to the large computational effort, the full numerical solution could be attainable.
6.3 Conclusions
In conclusion, the 3+1D CPIC method for simulating collisions within the CGC framework developed in this thesis was shown to be viable and in agreement with the boost invariant scenario as a limiting case. It genuinely extends the boost invariant method by including the longitudinal coordinate and is able to describe the fully rapidity-dependent dynamics of the Glasma. In particular, it allows for studying the effects of finite longitudinal width of the colliding nuclei and their longitudinal structure. Other extensions, such as axial charge production, the rapidity dependence of gluon occupation numbers, and quark pair production in the Glasma background field, could provide interesting new insights into the pre-equilibrium stage of heavy-ion collisions in future works.
Appendix A Notation and conventions
In this chapter of the appendix we summarize special notation and conventions that are used throughout the thesis.
A.1 Natural units
In this thesis we use natural units as defined by and . This allows us to express particle masses and momenta in units of energy. In particle physics the preferred unit of energy is the electron volt denoted by eV. For example, the mass of the proton is . Length and time scales are given in inverse units of energy, i.e. . It is convenient to also use femtometers (or Fermi) for lengths, e.g. the radius of a nucleus with mass number is given by . A gold nucleus with mass number has a radius of .
Using the relation
[TABLE]
one can translate between inverse energy units and Fermi. For a gold nucleus one finds . It is common to express time scales in Fermi divided by the speed of light , e.g. . In some cases, such as energy densities , we use a mix of Fermi and electron volts:
[TABLE]
A.2 Yang-Mills theory
Throughout the thesis we use the mostly minus metric convention
[TABLE]
where Greek indices run from [math] to . The symbol is used as a shorthand for the partial derivative with respect to the coordinate .
The Lagrangian density of Yang-Mills theory with colors is given by
[TABLE]
are the components of the non-Abelian field strength tensor given by
[TABLE]
where is a non-Abelian gauge field and is the Yang-Mills coupling constant. are the real-valued, totally antisymmetric structure constants of the Lie algebra of the Lie group . The color component index runs from to .
Often it is convenient to work with algebra elements instead of color components . The gauge field can be written as
[TABLE]
where are the traceless, hermitian generators of the gauge group . In the fundamental representation of , we use the normalization
[TABLE]
The commutator of two generators is related to the antisymmetric structure constants via
[TABLE]
The anti-commutator defines the real-valued, symmetric structure constants
[TABLE]
Using generators we can write the field strength tensor as an algebra element
[TABLE]
The Lagrangian density can be expressed as
[TABLE]
The Lagrangian density of Yang-Mills theory is invariant under local gauge transformations. The gauge field transforms according to
[TABLE]
where is a space-time dependent gauge transformation. It holds that
[TABLE]
The gauge transformation of the field strength tensor is given by
[TABLE]
Compared to eq. A.12 the transformation of does not involve any partial derivatives acting on . The invariance of is easy to show using the above transformation law and the cyclicity of the trace
[TABLE]
The gauge-covariant derivative , acting on an algebra element , is given by
[TABLE]
If transforms like
[TABLE]
then transforms in a similar fashion
[TABLE]
where we have used .
The anti-hermitian, traceless part of a matrix in the fundamental representation is given by
[TABLE]
It holds that
[TABLE]
and
[TABLE]
If is hermitian, i.e. , then . The anti-hermitian traceless part is related to the components . They are defined by
[TABLE]
where “” denotes the imaginary part. For the fundamental representation of , it holds that
[TABLE]
which can be shown using the Fierz identity for the generators of
[TABLE]
where are fundamental representation matrix indices.
A.3 Lattice notation
In order to make calculations in discretized systems less cumbersome, we use a compact notation for discretized fields on lattices. We define the four-dimensional, regular hypercubic lattice as the set of lattice sites:
[TABLE]
where (no sum implied) with unit vectors . The time step is (sometimes ) and the spatial lattice spacings are . We use a shorthand notation to denote the neighbors of a lattice site . For example, corresponds to the lattice site shifted by . A set of points which share the same time coordinate are referred to as a time slice.
A discretized field is the set of field values defined at every lattice site . The neighbors of are denoted as . For example, is the field shifted by one time step. is shifted by one lattice spacing along the axis.
Forward and backward finite differences are defined as
[TABLE]
We use the “” symbol to denote equivalence under the sum over all lattice sites. For example, using explicit sums it holds that
[TABLE]
where we have simply used the shift . We write “” when we omit the explicit sum, i.e.
[TABLE]
If is matrix valued (e.g. when working with real-time lattice gauge theory) and an expression involves a trace, “” also implies equivalence under the trace. For example the expression
[TABLE]
is shorthand for
[TABLE]
where we have used the cyclic property of the trace.
As an example of a derivation using this shorthand notation, we prove a formula for summation by parts, which is the discrete analogue of integration by parts. It holds that
[TABLE]
To see this, we write first write the LHS of the above equation as
[TABLE]
and then perform simple manipulations using the “” symbol:
[TABLE]
or simply
[TABLE]
Appendix B Boost invariant Yang-Mills theory
B.1 Yang-Mills theory in , coordinates
The Yang-Mills action in flat coordinates (e.g. laboratory frame coordinates ) is given by
[TABLE]
The goal is to rewrite this action in terms of the coordinate set, which is defined via
[TABLE]
In terms of light cone coordinates
[TABLE]
one finds
[TABLE]
The gauge field transforms as a four-vector:
[TABLE]
The temporal component expressed in terms of and is therefore
[TABLE]
Similarly, one finds
[TABLE]
In , coordinates we have the infinitesimal line element
[TABLE]
which implies the metric
[TABLE]
The infinitesimal space-time volume is given by . Therefore we can write
[TABLE]
Lowering all indices we find
[TABLE]
One could define the electric fields in the frame as
[TABLE]
and thus we can rewrite the action as
[TABLE]
However, interpreting as the actual longitudinal electric field is problematic, since it only has units of energy instead of units of energy squared. Relating the field strength tensor in coordinates to the laboratory frame field strengths one finds [74]
[TABLE]
The additional factor in the relation for restores the correct units. We also define the canonical momenta using functional derivatives
[TABLE]
We start varying the action (B.13) with respect to without specifying any gauge condition. This yields
[TABLE]
Integrating by parts and making use of gives
[TABLE]
which by requiring yields the Gauss constraint
[TABLE]
or in terms of canonical momenta
[TABLE]
The only remaining variations to be done are with respect to and . We perform both already employing Fock-Schwinger gauge
[TABLE]
which is equivalent to temporal gauge . Furthermore, we require that the gauge fields are boost invariant, i.e. they do not depend on . The action then simply reads
[TABLE]
In order to always guarantee that the fields are independent of , gauge transformations must be independent of as well. This restriction affects the transformation behavior of , which now transforms as a scalar
[TABLE]
In fact, the action (B.28) is equivalent to the action of a boost invariant, color charged scalar:
[TABLE]
with .
Varying the action (B.28) with respect to yields
[TABLE]
while varying with respect to gives
[TABLE]
The equations of motion therefore read
[TABLE]
or in terms of canonical momenta
[TABLE]
with
[TABLE]
Under the assumption of boost invariance the Gauss constraint (B.26) simplifies to
[TABLE]
B.2 Boost invariant real-time lattice gauge theory
In order to solve the boost invariant Yang-Mills equations numerically, we employ real-time lattice gauge theory techniques. In this section we present a non-rigorous, but quick derivation of the discretized equations. A more thorough discussion of real-time lattice gauge theory techniques follows in appendix C.
Starting from the action (B.28), we replace the continuous transverse plane with a finite, rectangular grid while keeping proper time continuous. Instead of the transverse gauge field , we use the transverse gauge links
[TABLE]
where the index refers to the coordinate in the transverse plane and is the lattice spacing in the direction. The gauge links are elements of the gauge group in the fundamental representation, i.e. unitary matrices. As a conjugated momentum we use (instead of the continuous canonical momentum ) which we define to be acting on the left via
[TABLE]
which reduces to eq. B.38 in the continuum limit . Note that is an algebra element of the gauge group in the fundamental representation, i.e. a hermitian, traceless matrix. Consequently, it is guaranteed that the matrices remain unitary under time evolution, since
[TABLE]
and similarly .
The gauge links are to be interpreted as Wilson lines connecting nearest neighbors on the rectangular lattice. Under discrete gauge transformations (i.e. gauge transformations defined at the discrete lattice sites ), the gauge links transform according to
[TABLE]
where refers to the lattice site shifted by one cell in the direction of . Since we are employing temporal gauge , the discrete gauge transformations must be independent of . Furthermore, since we implement boost invariance at the level of gauge fields, gauge transformations are also independent of .
The smallest possible Wilson loops that can be formed on the lattice are the transverse plaquettes
[TABLE]
where we define . The plaquette at transforms locally, i.e.
[TABLE]
In the continuum limit the plaquette can be related to the field strength tensor (see eq. 2.102)
[TABLE]
which yields the approximate, gauge invariant relation (see eq. 2.103)
[TABLE]
On the other hand, the gauge field component does not need to be replaced by a link, since the direction is kept continuous. It is simply replaced by defined at the transverse lattice sites . Under gauge transformations, this field transforms locally
[TABLE]
since the transformations are by definition independent of . The discretized canonical momentum of is simply defined by
[TABLE]
Gauge-covariant derivatives acting on in the transverse plane have to be replaced by gauge-covariant finite difference expressions. We define the first-order forward and backward derivatives
[TABLE]
The (semi-)discretized, boost invariant action in temporal gauge can now be defined as
[TABLE]
and the discretized Gauss constraint reads
[TABLE]
In order to derive the discretized equations of motion we have to vary the above action with respect to and . The variation with respect to the scalar component is performed in a straightforward way with . The variation then reads
[TABLE]
The terms involving the derivative can be integrated by parts in the usual way. For the discrete gauge-covariant differences we can perform a summation by parts, i.e.
[TABLE]
where the forward difference turns into a backward difference. Defining (no sum implied), the variation then reads
[TABLE]
and with we can read off the discrete equation of motion
[TABLE]
which is completely analogous to the continuous equation (B.36).
The variation of section B.2 with respect to the transverse gauge links is slightly more involved. It is insufficient to simply vary the matrix elements of independently, because gauge links are elements of , i.e. they satisfy the constraints and . However, we can perform a variation that leaves these constraints fulfilled (see section C.2):
[TABLE]
Starting with the kinetic term of section B.2, we find after some algebra and integration by parts
[TABLE]
Varying the plaquette term of section B.2 we find
[TABLE]
With the shorthand
[TABLE]
we can write this as
[TABLE]
Varying the term involving of section B.2 with respect to we find
[TABLE]
where we used the shorthand notation .
Combining all terms, the variation of the action reads
[TABLE]
With the equations of motion for read
[TABLE]
Multiplying everything with and using (see eq. A.24)
[TABLE]
and
[TABLE]
we can write this equation directly in terms of algebra and group elements
[TABLE]
While the continuum limit of the commutator term in the above equation on the right is obvious, the term involving plaquettes is more subtle. Using and , we can expand for small lattice spacing to arrive at
[TABLE]
Therefore the equations of motion for reduce to eq. B.35 in the continuum limit .
Finally, in order to solve the time evolution numerically, the time coordinate also has to be discretized in eqs. (B.41), (B.49), (B.53), (B.57) and (B.68). A simple approach is to use a leapfrog time integrator. We introduce the finite time step and discrete times . We then define that the scalar component and the transverse gauge links are evaluated at whole-numbered time steps , while the canonical momenta and are evaluated in-between these steps, i.e. at . The discretized equations (i.e. update equations) then read
[TABLE]
and the discrete Gauss constraint consistent with these update equations is given by
[TABLE]
The Gauss constraint contains momenta evaluated at , while gauge fields and links are evaluated at . It is straightforward to prove that the Gauss constraint is conserved under the update equations.
Appendix C Real-time lattice gauge theory
C.1 Parametrization for
In this section of the appendix we list a few formulae and identities for the parametrization of the SU(2) group that is being used in the simulation. This parametrization can also be used for elements of and anti-hermitian, traceless matrices in .
We can represent any SU(2) group element using the exponential map applied to an algebra element. The dimension of the Lie algebra (2) of SU(2) is . Given the three-component vector and the generators with the Pauli matrices , we can write any element in the fundamental representation as
[TABLE]
The matrix is unitary and has determinant one. Exploiting the property of the Pauli matrices (no sum over implied)
[TABLE]
we write
[TABLE]
where we used and . The above exact relation provides a simple way of computing the exponential map in SU(2). Moreover, it allows us to represent any group element in SU(2) with only three real-valued parameters . Since computing trigonometric functions can be numerically costly, one can switch to a slightly different representation. By introducing the four real-valued parameters
[TABLE]
the group element can be written as
[TABLE]
The above is only a group element if the constraint
[TABLE]
holds. This parametrization is particularly useful, since it also allows us to represent other complex matrices as well, such as sums of unitary matrices. In the following sections we use more compact notation and drop the sum symbol.
C.1.1 Multiplication
Scalar multiplication is very simple in terms of parameters and . Let be a real number. We then have
[TABLE]
We just have to multiply the parameter vector by the scalar . Note that the norm of the parameter vector changes to
[TABLE]
Matrix multiplication in terms of parameter vectors is more complicated. Let and be two group elements with parameter vectors and respectively. Using the property
[TABLE]
we find
[TABLE]
Let be the parameter vector for . We then find
[TABLE]
The matrix is again a matrix due to
[TABLE]
and consequently
[TABLE]
If and are elements of (i.e. and ), then will also be in (i.e. ). This multiplication rule also works for general complex matrices using this parametrization.
C.1.2 Inverse of the exponential map
If we want to extract the original algebra vector from the parameters (see eqs. (C.5) and (C.6)), i.e. perform the inverse exponential map, then we have to do the following: First, we find the normalized vector using
[TABLE]
We can compute the norm from eq. (C.6)
[TABLE]
as long as . The original components are then given by
[TABLE]
If , one can also linearize:
[TABLE]
We then get
[TABLE]
C.1.3 Anti-hermitian traceless part of group elements
The anti-hermitian traceless part of is given by (see eq. A.20)
[TABLE]
Using , we find
[TABLE]
and therefore
[TABLE]
Similarly, we find
[TABLE]
It holds that
[TABLE]
This is an exact relation, which generally holds for the fundamental representations of (see section C.4).
In the case of SU(2), the relation between and also allows us to write in terms of . Using and assuming we find from eq. C.8
[TABLE]
C.1.4 Inverse of matrices
If SU(2), then is unitary and therefore . In terms of our parametrization we then have
[TABLE]
However, if is not an element of SU(2) (i.e. ), then we can still use the parametrization to compute the inverse. Writing
[TABLE]
and requiring we find
[TABLE]
with .
This is mainly used for the Glasma initial conditions on the lattice, see sections 3.3 and 2.135.
C.2 Variation of gauge links
When performing the variation of e.g. the Wilson action (see section 2.3.2) with respect to gauge links , we have to keep in mind that one considers a variation w.r.t. special unitary matrices, i.e. it holds that
[TABLE]
It is therefore incorrect to simply vary as a function of the matrix elements of . Instead, we must make sure that the constraints in eqs. (C.36) and (C.37) remain fulfilled.
For instance, imagine writing the Wilson action for in terms of the parametrization of SU(2) introduced in section C.1, such that every gauge link is given in terms of the four real-valued parameters and () with . One could then try to find the variation of by taking the independent derivatives w.r.t. the parameters and for each link. This would yield the wrong result, because the parametrization using four parameters is in fact overdetermined. In order for to be an element of SU(2), the parameters satisfy eq. C.8
[TABLE]
which describes a three-dimensional sphere embedded in . Taking independent derivatives w.r.t. and generally violates this constraint. Geometrically speaking, one would “leave” the sphere instead of “moving along” the sphere while performing the variation.
In order to solve this problem, we could implement the constraints using the method of Lagrange multipliers. However, this can be quite cumbersome: for each link we would need to introduce independent multipliers for eqs. (C.36) and (C.37) in the Wilson action and then additionally solve for the unknown multipliers.
In order to avoid the use of Lagrange multipliers, we use the method of constrained variation, which turns out to be much easier to handle. We introduce the infinitesimal variation of a gauge link variable
[TABLE]
or written shorter
[TABLE]
where the variation of the gauge field is traceless and hermitian and the factor is absorbed in a re-definition of . In the continuum limit becomes the infinitesimal variation of the gauge field . The infinitesimal variation preserves the unitarity of gauge links. Let . We then find
[TABLE]
The determinant is also unaffected for infinitesimal variations. Using Taylor expansion of the determinant we find
[TABLE]
where “” refers to the adjugate which is given by . In the second line of section C.2 we have made use of Jacobi’s formula for the derivative of the determinant. The variation therefore preserves the constraints and allows us to vary the action without the use of Lagrange multipliers, which dramatically simplifies the derivation of equations of motion.
The constrained variation of the gauge link further allows us to define the derivative
[TABLE]
This follows from
[TABLE]
C.3 Variation of the leapfrog action
In this section of the appendix we derive the discrete equations of motion (EOM) (as used in sections 3.2.1 and 5.5.1) obtained from the standard Wilson action (2.3.2) (or eq. 5.129) using the constraint preserving variation of link variables of the previous section. To make the calculation more organized, we first split the action into two parts: a part containing temporal plaquettes (“E” for electric) and a part containing spatial plaquettes (“B” for magnetic). We write
[TABLE]
where
[TABLE]
Furthermore, we use the identity introduced in section 5.5
[TABLE]
with
[TABLE]
This allows us to rewrite the electric and magnetic part as
[TABLE]
C.3.1 Gauss constraint
For the Gauss constraint (GC) in the leapfrog scheme we only have to consider the variation , since does not contain any temporal links. In the following sections we make use of the “” symbol, denoting equality under the sum over lattice sites and under the trace (see section A.3). We then have
[TABLE]
To go from the first to the second line, we applied a shift in the right term of the first line and made use of the cyclicity of the trace. In the third line we simply used the definition of the variation of gauge links. The variation of the action therefore reads
[TABLE]
where we used
[TABLE]
Replacing the terms with link variables we find
[TABLE]
and subsequently
[TABLE]
We require that the variation vanishes, i.e. . Since all gauge links can be varied independently we find
[TABLE]
C.3.2 Equations of motion
First, we consider the variation of w.r.t. spatial links. We find a result that is similar to the expression for the GC
[TABLE]
and
[TABLE]
For the variation of we use
[TABLE]
The variation then reads
[TABLE]
We set and, after canceling some constants, we find the discrete EOM
[TABLE]
which can also be written as
[TABLE]
C.4 Conservation of the Gauss constraint
We now explicitly show that the leapfrog EOM (C.63) preserve the associated GC (C.57). We use the identity for the fundamental representation of (see eq. A.24)
[TABLE]
where “ah” denotes the anti-hermitian traceless part of , see eq. A.20. The constraint and the equations of motion (including external charges) then read
[TABLE]
We take eq. C.66 and sum over . Due to being antisymmetric in the index pair we find
[TABLE]
Doing the same at and parallel transporting from to yields
[TABLE]
Subtracting the above two equations gives
[TABLE]
Using antisymmetry we have
[TABLE]
and
[TABLE]
Moreover, in temporal gauge we have and , which leads to
[TABLE]
where we used the GC in the last line to replace the temporal plaquette terms with charge densities. This yields the gauge-covariant continuity equation
[TABLE]
If there are no external charges, we simply have the conservation of the GC: assuming that the constraint in the previous time slice holds, i.e.
[TABLE]
then the EOM guarantee that it will also hold in the next one, i.e.
[TABLE]
C.5 Variation of the implicit action
We now consider the action introduced in section 5.5.2
[TABLE]
where
[TABLE]
The electric part is the same as in the leapfrog scheme. However, the magnetic part now contains temporal gauge links and gives a new contribution to the GC.
Before we vary this action, we must verify that is indeed real-valued. While it is easy to see that the original leapfrog action is real-valued because of the obvious hermicity of , the term is not hermitian in general. Still, we can show that the action is real: we have
[TABLE]
where we can rewrite
[TABLE]
Here we used the cyclicity of the trace. Then using a shift we have
[TABLE]
Likewise we have
[TABLE]
which leads to
[TABLE]
Incidentally, the RHS term is exactly the hermitian conjugate of the LHS term. In other words
[TABLE]
This means that under the sum over and the trace, the expression is indeed real-valued and by extension is real-valued as well. We have also shown that the time-average in can be “shifted” to the other term under the sum and trace, i.e.
[TABLE]
This is a useful property that we will need in the following derivation.
C.5.1 Gauss constraint
We perform the variation of and w.r.t. temporal links to derive the GC in the implicit scheme. Since is the same for all schemes, we do not have to repeat it. On the other hand, the variation of the magnetic part involves terms like , which we now discuss explicitly. First, we make use of
[TABLE]
Then, after some algebra we find
[TABLE]
and
[TABLE]
which yields
[TABLE]
The variation of the magnetic part therefore reads
[TABLE]
In the last line we consolidated terms with index and into a single term using the sum . Taking the result for from the leapfrog scheme, we find the GC
[TABLE]
C.5.2 Equations of motion
For the EOM we vary w.r.t. spatial links. Again, we already have the result for from the leapfrog scheme and only need to calculate . In particular we consider the term
[TABLE]
Since the variation only acts on spatial links we can shift the time-average of the right term from to . This gives
[TABLE]
The variation then proceeds analogously to the derivation of the leapfrog scheme in section C.3.2. We find
[TABLE]
With the result for we obtain the discrete EOM
[TABLE]
Introducing the shorthand
[TABLE]
allows us to write the EOM rather compactly as
[TABLE]
C.6 Variation of the semi-implicit action
In the semi-implicit scheme (see section 5.5.3) the action reads
[TABLE]
where is given by eq. C.50. The magnetic part comprises of , where is the same as from the implicit scheme except that the indices and in only run through transverse components:
[TABLE]
Therefore we can take the results from the previous section for , eqs. (C.5.1) and (C.94). The new part is given by
[TABLE]
C.6.1 Gauss constraint
We already know and from sections C.3.1 and C.5.1, so we only have to compute . The relevant terms are
[TABLE]
where we defined
[TABLE]
Using the same techniques as before, we find
[TABLE]
and the variation of reads
[TABLE]
which (with the previous results taken into account) yields the GC in the semi-implicit scheme, see section 5.5.3.
C.6.2 Equations of motion
For the variation w.r.t. spatial links we have to distinguish two cases: the longitudinal and the transverse links. Starting with the variation of longitudinal links we find
[TABLE]
Again we used the fact that the time-average can be shifted to other terms by exploiting the sum over and the cyclicity of the trace. The variation of then reads
[TABLE]
Combining the above with yields the longitudinal EOM (5.5.3). No contributions from are necessary because it does not include any longitudinal links.
For the transverse components of the EOM we vary w.r.t. , where is a transverse index. The relevant terms for are
[TABLE]
and
[TABLE]
For terms with negative component indices we can show that they are identical (under the sum over and the trace) to the last two terms except for the substitution :
[TABLE]
This allows us to write
[TABLE]
where stands for summing over terms with component indices and . The variation of then reads
[TABLE]
With the expressions for and we find the transverse components of the EOM (5.5.3).
Appendix D Stability of the semi-implicit Abelian scheme
In this section of the appendix we prove stability for the semi-implicit scheme for Abelian fields derived in section 5.4.3. Using temporal gauge, , the EOM (5.108) and (5.109) read
[TABLE]
where reduces to
[TABLE]
The Gauss constraint (5.107) reads
[TABLE]
Splitting the equations of motion into Laplacian terms and mixed derivative terms we find
[TABLE]
Using a plane wave ansatz
[TABLE]
we will use the Gauss constraint to first reduce the number of degrees of freedom and then compute the dispersion relation . Inserting the ansatz into the Gauss constraint we find
[TABLE]
where we use the dimensionless lattice momenta (5.85). The constraint equation can alternatively be written as
[TABLE]
where is a momentum-dependent factor given by
[TABLE]
The temporal average reduces to a multiplication with a frequency dependent factor:
[TABLE]
where we used the shorthand . Inserting the plane wave ansatz into the EOM yields
[TABLE]
Note that both and depend on . After making use of the Gauss constraint the longitudinal EOM reads
[TABLE]
and the two transverse equations read
[TABLE]
This system of equations can be written in matrix notation as an eigenvalue problem
[TABLE]
where the coefficient matrix is given by
[TABLE]
and the vector is simply
[TABLE]
The eigenvectors of are
[TABLE]
where we find two transverse, momentum dependent eigenvectors , and the longitudinal unit vector . The two transverse vectors and can be interpreted as transverse polarization modes and are orthogonal in the sense of . The three eigenvectors yield three different equations for the eigenvalue problem, namely
[TABLE]
where and are expressions which depend on the momenta and the frequency via . Solving the first equation for yields the dispersion relation of the longitudinal component:
[TABLE]
Solving the two transverse equations yields
[TABLE]
and
[TABLE]
It turns out that the expressions for associated with the different eigenvectors are not the same, except that and share the same dispersion relation, that is to say in general. This implies that for a given momentum , the amplitude of an arbitrary wave has to be split into two components: a part which is projected into the plane spanned by and which oscillates with and a part parallel to which oscillates with frequency . We interpret this as numerical (or artificial) birefringence.
We require the propagation of a wave to be stable, i.e. we require the frequencies to be real-valued. This is guaranteed if the arguments of the expressions in eqs. (D.23) and (D.24) are restricted to . Both dispersion relations remain stable if the CFL condition
[TABLE]
holds. Using
[TABLE]
and requiring stability for all values of yields the constraint
[TABLE]
This concludes the proof that the semi-implicit scheme, even though exhibiting peculiar wave propagation phenomena, is stable.
Appendix E Abbreviations
2+1D
2+1 dimensions / dimensional
3+1D
3+1 dimensions / dimensional
BFKL
Balitsky, Fadin, Kuraev, Lipatov
CFL
Courant, Friedrichs, Levy
CGC
Color Glass Condensate
CIC
cloud-in-cell
CPIC
colored particle-in-cell
EOM
equation(s) of motion
FDTD
finite difference time domain
FFT
fast Fourier transformation
GC
Gauss constraint
IR
infrared
JIMWLK
Jalilian-Marian, Iancu, McLerran, Weigert, Leonidov, Kovner
LC
light cone
LF
leapfrog
LHC
Large Hadron Collider
LHS
left hand side
LRF
local rest frame
MV
McLerran, Venugopalan
NCI
numerical Cherenkov instability
NGP
nearest grid point
PIC
particle-in-cell
QCD
quantum chromodynamics
QGP
quark-gluon plasma
RHIC
Relativistic Heavy Ion Collider
RHS
right hand side
SI
semi-implicit
SU
special unitary
UV
ultraviolet
w.r.t.
with respect to
YM
Yang-Mills
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1[1] D. Gelfand, A. Ipp and D. Müller, Simulating collisions of thick nuclei in the color glass condensate framework , Phys. Rev. D 94 (2016) 014020 , [ 1605.07184 ]. · doi ↗
- 2[2] A. Ipp and D. Müller, Broken boost invariance in the Glasma via finite nuclei thickness , Phys. Lett. B 771 (2017) 74–79 , [ 1703.00017 ]. · doi ↗
- 3[3] A. Ipp and D. Müller, Rapidity profiles from 3+1D Glasma simulations with finite longitudinal thickness , Po S EPS-HEP 2017 (2017) 176 , [ 1710.01732 ]. · doi ↗
- 4[4] A. Ipp and D. Müller, Implicit schemes for real-time lattice gauge theory , Eur. Phys. J. C 78 (2018) 884 , [ 1804.01995 ]. · doi ↗
- 5[5] U. Heinz and R. Snellings, Collective flow and viscosity in relativistic heavy-ion collisions , Ann. Rev. Nucl. Part. Sci. 63 (2013) 123–151 , [ 1301.2826 ]. · doi ↗
- 6[6] N. Brambilla et al., QCD and Strongly Coupled Gauge Theories: Challenges and Perspectives , Eur. Phys. J. C 74 (2014) 2981 , [ 1404.3723 ]. · doi ↗
- 7[7] PHENIX collaboration, K. Adcox et al., Formation of dense partonic matter in relativistic nucleus-nucleus collisions at RHIC: Experimental evaluation by the PHENIX collaboration , Nucl. Phys. A 757 (2005) 184–283 , [ nucl-ex/0410003 ]. · doi ↗
- 8[8] BRAHMS collaboration, I. Arsene et al., Quark gluon plasma and color glass condensate at RHIC? The Perspective from the BRAHMS experiment , Nucl. Phys. A 757 (2005) 1–27 , [ nucl-ex/0410020 ]. · doi ↗
