Coloured Particle-in-Cell simulations of
non-Abelian plasma dynamics in glasma
Graduate School of Science, Hiroshima University
Department of Physical Science

M162632 Tomohiro Honda

Supervisor : Prof. Toru Sugitate
Chief examiner : Prof. Toru Sugitate
Second examiner : Associate Prof. Takuya Morozumi
April 2, 2018

Abstract
Ultrarelativistic heavy ion collisions create quark gluon plasma(QGP), a
deconfined state of quarks and gluons. There have been a lot of experiments to understand the nature of QGP. From these attempts, elliptic flowof QGP was found and it is an evidence for a collective motion of partons.
The results are consistent with hydrodynamical simulations, therefore, it isconsidered that QGP is in a local equilibrium state. Its thermalisation process, however, is still one of the big questions in heavy ion collisions. One
possible explanation is that plasma instabilities take significant role in thermalising preliminary matter, glasma. Hence, it is important to understandthe plasma dynamics in glasma for exploring the very beginning of heavy ion
collisions. In this thesis, I investigate fundamental physics of plasma dynamics in glasma. Assuming high energy limit, high momentum gluons can betreated as plasma particles, low momentum ones as background fields. For
basics equations of calculation, I adopted one dimensional Wong-Yang-Mills
equations, based on Vlasov equation and Yang-Mills equations. Analysing
kinetic evolution of plasma numerically, I implemented coloured Particle-inCell(PIC) method which is the expansion of PIC method commonly used
for plasma simulations.
I set initial condition such as counter flow of particles and a few types ofelectro-magnetic fields. I discuss the time developments of particle distribution, electro-magnetic field, and the eﬀect of nonlinear term in Yang-Mills
equations. In the case of zero initial fields, hole structure in phase space
was constructed as is the case with classical electro-magnetic(EM) plasma.
This results shows the possibility that phenomena similar to classical electromagnetic plasma can be observed with small eﬀects of nonlinear terms. In
order to verify the applicability of various concepts in classical plasma, it isnecessary to increase the accuracy and the spatial dimension of simulations.

1Contents
Abstract

1

1 Introduction
1.1 Heavy ion collisions . . . . . . . . . . . . . . . . . . . . . . . .
1.2 Quark gluon plasma . . . . . . . . . . . . . . . . . . . . . . .
1.3 Colour glass condensate . . . . . . . . . . . . . . . . . . . . .
1.4 Thermalisation process . . . . . . . . . . . . . . . . . . . . . .
1.5 objective . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3
3
4
5
5
6

2 Basic concepts
2.1 Yang-Mills equations . . . . . . . . . . . . . . . . . . . . . .
2.2 Wong’s equations . . . . . . . . . . . . . . . . . . . . . . . . .

7793 Simulation methods
11
3.1 Particle-in-Cell method . . . . . . . . . . . . . . . . . . . . . 11
3.2 Discretisation of basic equations . . . . . . . . . . . . . . . . 12
3.3 Colour-electric current on a grid . . . . . . . . . . . . . . . . 15
4 Results
16
4.1 Results without nonlinear eﬀects . . . . . . . . . . . . . . . . 16
4.2 Results for Yang-Mills plasma . . . . . . . . . . . . . . . . . . 17
5 Conclusion

20

Acknowledgements

212

Chapter 1

Introduction
In this chapter, I introduce background knowledge of my research.

1.1

Heavy ion collisions

Ultra-relativistic heavy ion collisions are extreme high energy physics on theearth. Experiments create quark gluon matter with high energy and density,
called quark gluon plasma(QGP). It is considered that QGP was existed in
the very beginning of the universe. Therefore, research of QGP dynamics
has significant role for going way back to the early time physics. There are

Figure 1.1: Time evolution of heavy ion collision.

several steps in heavy ion collisions as following: (1)First, accelerated nuclei
at almost speed of light collide. They are pancake-like shape due to Lorentz
3contraction.
(2)The energy density in the centre of collision get extremely high.
(3)Through thermalisation process, quarks and gluons go through phasetransition to be in deconfined state.
(4)Haronisation occurs as the energy density decreases due to its expansion.
Finally chemical freezeout takes place, then, kinds and number of hadron isdetermined.

1.2

Quark gluon plasma

QGP is deconfined state of quarks and gluons at extremely high temperatureand energy density. Quarks and gluons are in confined state in hadronic sector, called hadron, however, phase transition to QGP takes place in extreme
condition.
QGP is considered to be in local equilibrium state from experimental results. Lattice simulation[1] shows saturation of entropy density of quark
gluon matter. Hydrodynamical simulation provide a good explanation for a

Figure 1.2: Entropy density saturate as time advances. It is one of theevidence of thermalisation of colliding nuclei[1]lot of experimental data, such as radial and elliptic flow.

4

1.3

Colour glass condensate

QCD theory shows the possibility of asymptotic freedom in which the coupling constant decreases at the limit of high energy. In such weak coupling
multi-body system, it is considered that a state called colour glass condensate exists. For such state, L. McLerran and R. Venugopalan gave a
explanation, MV model[5]. This description tell that a lot of gluons are
created in heavy ion collision due to its high energy density.
Fig1.3 show the parton distribution in proton while deep inelastic scatter-

Figure 1.3: Parton distribution in proton, measured in deep inelastic scattering at HERA.[2]ing at HERA. In fact, lattice calculation for parton distributions at highenergy shows that gluons are outnumbering all the other particles. Therefore, valence quarks are negligible compared with gluons, and sea quarks aresuppressed by one power of the coupling constant since they are produced
from gluons g → q q̄.

1.4

Thermalisation process

One of the big question of heavy ion collision and creation of QGP is its
thermalisation process of quark gluon matter. From Hydrodynamics, ther-5malisation time scale is about 0.6fm/c. However, calculation based on perturbative scattering show its time scale a few times larger than experimental data. One possible reason why scattering process estimate larger results
is that this calculation ignores collective behavior of plasma, for example,
plasma instabilities.

1.5

objective

The purpose of this thesis is to analyse the fundamental physics of nonAbelian plasma dynamics. The dynamics is important for investigating themechanism of heavy ion collisions and the very beginning of the universe.
Yang-Mills plasma dynamics for counter flow initial conditions are shown inthis thesis.

6Chapter 2

Basic concepts
2.1

Yang-Mills equations

I consider the system with only gluons because gluons occupy the centre
part of collision and approximate the fields as classical one in the following.
Because of low coupling constant, I divide into gluons with larger momentum and smaller momentum. In other words, I approximate gluons withlarge and small momentum to classical plasma particles and classical gauge
fields, respectively. Metric of the system is (1,-1,-1,-1). Classical Yang-Mills
equations in CGS unit is√1g√ ∂ν Faµν − √ fabc Abν Fcµν = − 4πgjaµ
4π
4π
a
a
a
Fµν = ∂ν Aν − ∂ν − gfabc Abµ Acν ,

(2.1)
(2.2)

where Faµν is electro-magnetic tensor of gluon fields and Aaν is vector potential of gluon fields, jaµ is colour electric current. In addition, a is inertial
degrees of freedom of gluons with 8 kinds for SU(3) theory. Here, lettingf (x, p, Q) be the distribution function in phase space for one particle, jaµ ,Faµν and Aaν are∫µja (x) = dpdQpµ Qa f (x, p, Q)
(2.3)
(2.4)


0

 −ExaFaµν = 
 Eya
Eza


Exa EyaEza
0 −Bza
Bya 

a
Bz
0
−Bxa 
−Bya Bxa0Aν = (ϕ, A).

(2.5)

7

I use these equations for colour particle simulations of gluon system. First,
I discretise continuous distribution function into particles. Then f is
f (x, p, Q) =(2π)3 ∑
δ(x − xi (t))δ(p − pi (t))δ(Q − Qi (t)),
Ntest

(2.6)

i

where xi (t), pi (t), Qi (t) are position, momentum and colour charge of theith particle. Therefore, colour electric current is∑
ja = gQai vi δ(x − xi ).
(2.7)
ivi is velocity of ith particle , and the absolute value is the same as light.
I use temporal gauge for gauge condition
Aa0 = 0.

(2.8)

From classical Yang-Mills equation and gauge condition,∂Aa= −E a∂t(2.9), and magnetic fields are
Bxa = −∂z Aay + ∂y Aaz − gfabc Abz Acy(2.10)

Bya = −∂x Aaz + ∂z Aax − gfabc Abx Acz
Bza = −∂y Aax + ∂x Aay − gfabc Aby Acx(2.11)(2.12)

Then, using fabc Abz Acy = − 12 fabc [Ab × Ac ]x from perfect antisymmetry of
structure constant fabc , colour magnetic field is1
B a = ∇ × Aa + gfabc Ab × Ac .
2

(2.13)

Substituting µ = 0 into (2.2),∂ν Fa0ν − gfabc Abν Fc0ν = −4πgja0 ,(2.14)

therefore, using Fa00 = 0,
∂1 Fa01 + ∂2 Fa02 + ∂3 Fa03 − gfabc (Ab1 Fc01 + Ab2 Fc02 + Ab3 Fc03 ) = −4πρa .
(2.15)
This is Gauss’s law for Yang-Mills equation.
∇ · E a + gfabc Ab · E c = 4πρa(2.16)
8

Then, for the case of µ = 1 (2.2) is∂ν Fa1ν − gfabc Abν Fc1ν = −4πjxa(2.17)

Substituting gauge condition and the values of fields, the x component is∂Exa= ∂y Bza − ∂z Bya + gfabc (Aby Bzc − Abz Byc ) − 4πjxa .
∂t(2.18)

Adopting the same operation to µ = 2, 3, we obtain∂t E a = ∇ × B a + gfabc Ab × B c − 4πj.

(2.19)

This is Ampere-Maxwell law of Yang-Mills equation. Then, diﬀerentiating
eq.(2.13), we obtain
1∂t B a = ∇ × ∂t Aa + gfabc (∂t Ab × B c + Ab × ∂t Ac ).
2(2.20)∂t B a = −∇ × E a − gfabc E b × Ac

(2.21)

then,This is Faraday’s law of Yang-Mills equation. Finally, from the divergenceof eq. (2.13), we obtain Gauss’s law of magnetic field.
∇ · B a = −gfabc Ab · ∇ × Ac

(2.22)

These equations and Maxwell equation are the same when the nonlinear
term is zero. Therefore the nonlinear term is important for the dynamics ofheavy ion collision.

2.2

Wong’s equations

For large momentum gluons, I deal with them as plasma particle. In thissection, I introduce the equations of motion for plasma particles.
Considering zero mass particles flying at the speed of light, the equation of
motion, called Wong’s equations are as follows.
dxi
= vi
dt
dpi
= gQai (E a + vi × B a )dt
dQai= gfabc vi · Aa Qci
dt
p
v=
ϵ

(2.23)(2.24)(2.25)(2.26)

9Here, ϵ is energy of particles.
One diﬀerence between classical Maxwell particles and colour particles is
the time evolution of colour charge Qa . For SU(N) gauge fields, Qa satisfies
the equation below.
Qa Qa = N.

(2.27)The sum of colour charges is conservative for time evolution.

10Chapter 3

Simulation methods
3.1

Particle-in-Cell method

To solve the motion of plasma particles, we need information of electromagnetic field at positions of particles. Particle positions are continuous,
however, field is defined only on grids. In addition, we need to know electric
current to solve equations showed in previous chapter.
Here, field Fp at particle position xp is∑∑
Fp (xp ) =F (xi )S(xi − xp ) =Fi Si (xp ),
(3.1)
i

i

using F (xi ) defined on the grid. S is shape factor of particles. I use nearest
grid point method in my research defined as below.
{
1
S=
0

|x −x |

if p∆x i ≤ 12
else

(3.2)

Figure 3.1: Nearest grid point method.

11

3.2

Discretisation of basic equations

In this section, I discretise basic equations in terms of Particle-in-Cell method.
Basic equations are classical Yang-Mills equation for colour electro-magnetic
fields and Wong’s equations for motion of plasma particles.

Figure 3.2: electro-magnetic fields on the grid.[4]∂Aa= −E a∂t(3.3)

1
B a = ∇ × Aa + ϵabc Ab × Ac
2
∇ · E a + ϵabc Ab · E c = ρa

(3.4)
(3.5)

∂t E a = ∇ × B a + ϵabc Ab × B c − j.

dx
=v
dt
dp
= gQa (E a + v × B a )dt
dQa= gfabc v · Ab Qc
dt
p
p
v= =
ϵ
|p|
dϵ
= Qa v · E a
dt

(3.6)

(3.7)
(3.8)
(3.9)
(3.10)(3.11)

12

Calculation step is as follows:(1)update of electric fields(2)update of vector potential(3)update of magnetic fields
We can adopt central diﬀerence calculus:-Assigning E, A, j to integer grids, B to half integer grids with respect tothe spatial direction.
-Assigning E to integer grids, BjA to half integer grids with respect to thetime direction.
(i is spatial grid, n is time step)
a,n
a,n−1
Ex,i
− Ex,i

∆t
a,n
a,n−1
Ey,i
− Ey,i∆tb,n+ 12c,n+ 12= ϵabc (Ay,iBz,ia,n+ 1=−b,n+ 21− Az,ic,n+ 12By,ia,n+ 12) − jx,i(3.12)

a,n+ 1Bz,i+ 12 − Bz,i−]f2rac122∆xb,n+ 12+ ϵabc (Az,ic,n+ 21Bx,ib,n+ 12− Ax,ic,n+ 12) − jy,ic,n+ 12) − jz,iBz,ia,n+ 12(3.13)
a,n
a,n−1
Ez,i
− Ez,i∆t=−a,n+ 12y,i+ 12B

a,n+ 21− By,i−]f rac12∆xb,n+ 12+ ϵabc (Ax,ic,n+ 21By,ib,n+ 21− Ay,i(3.14)
a,n+ 21Ax,ia,n− 12− Ax,ia,n= −Ex,i(3.15)

a,n= −Ey,i(3.16)

a,n= −Ez,i.

(3.17)

1b,n+ 1 c,n+ 1
a,n+ 1b,n+ 1 c,n+ 1
Bx,i+ 12 = ϵabc (Ay,i+ 12 Az,i+ 12 − Az,i+ 12 Ay,i+ 12 )
2
2
2
2
2
2

(3.18)∆ta,n+ 21Ay,ia,n− 12− Ay,i∆ta,n+ 21Az,ia,n− 12− Az,i∆t

Ba,n+ 12y,i+ 12B

a,n+ 21z,i+ 12=−=−a,n+ 1a,n+ 1

2

2

Az,i+ 12 − Az,i− 12∆xa,n+ 1a,n+ 1

2

2

Ay,i+ 12 − Ay,i− 12∆x1b,n+ 1 c,n+ 1
b,n+ 1 c,n+ 1− ϵabc (Az,i+ 12 Ax,i+ 12 − Ax,i+ 12 Az,i+ 12 )
2
2
2
2
2
(3.19)
1b,n+ 1 c,n+ 1
b,n+ 1 c,n+ 1− ϵabc (Ax,i+ 12 Ay,i+ 12 − Ay,i+ 12 Ay,i+ 12 ).
2
2
2
2
2
(3.20)

13

Bx,ia,n+ 12.

Figure 3.3: Updates of physical values.[4]Then, I discretise the equation of motion of particle. We can adopt central
calculus:-Assigning Qa , x, ϵ to integer grids, p, v to half integer grids with respect tothe time direction.
11

pn+ 2 − pn− 2= Qna (Ean + v n × Ban ) + O(∆t2 )∆t(3.21)

We can obtain the values at integer time steps using backward diﬀerence.
11

pn+ 2 − pn− 2
p =+ O(∆t2 )2pn
vn = n .
|p |n(3.22)(3.23)

By backward diﬀerence, particle energy isϵn − ϵn−1∆2= Qna v · Ean + O(∆t).

(3.24)

1

Substituting v n = v n− 2 + O(∆t), we obtain111ϵn = ϵn− 2 + Qna ∆tv n− 2 · Ean + O(∆t2 ).
2(3.25)

Particle position is also obtained.
1

xn+1 = xn + v n+ 2 ∆t + O(∆t2 )14(3.26)

The diﬀerence between classical electro-magnetic plasma and Yang-Mills
plasma is time evolution of colour charge.
1
Qn+1
− Qna
n+ 1 n+ 1a= ϵabc v n+ 2 · Ab 2 Qc 2
∆t

3.3

(3.27)

Colour-electric current on a grid

We need to define electric current and charge density on grids because thefields are defined only on grids. I define electric current in this section.
Electric current for ith particle with colour a is defined as∑ qja = gQi vi δ(x − xi ).
(3.28)i

Defining shape factor S∫
Si,k =
Vk S0 (x − xi )dx,

(3.29)∆

I define electric current on grid k as the mean value of the cell. Then, we
obtaing ∑ a
jka =Qi vi Si,k .
(3.30)∆Vki

15

Chapter 4

Results
4.1

Results without nonlinear eﬀects

In this section, I show some results for classical Maxwell plasma to justify
the program. For this purpose I use pcans[2] code.
I employ electron-electron two steam instability. Two components of electron fluid have relative velocity, and there are finite mass ion in the background in initial condition. To, analyse the instability, wave number in theboundary of two components is important. When instability grows, large
wave number components increase in Fourier expansion. I obtained the results of growth curve from simulation and theory. I adopt Fourier expansion
to mode four and five with respect to space components.

Figure 4.1: Growth curve for mode 4,5. [4]16This results shows that calculation agrees the theoretical results. It justifies
results using my program.

4.2

Results for Yang-Mills plasma

In this section, I provide results for classical Yang-Mills plasma. Here, I
employ simple initial condition 1,2. initial condition 1 has zero initial fields,
and 2 has finite electric field for x direction.
Table 4.1: initial condition 1
directionx y z
electric field E
0 0 0magnetic field B
0 0 0
vector potential A 0 0 0

Table 4.2: initial condition 2
directionxy z
electric field E
0.5 0 0magnetic field B00 0
vector potential A0
0 0

Initial position distribution of particles are given uniformly and velocity
distribution is given by Maxwell distribution by random numbers in bothconditions.

Figure 4.2: left: initial distribution of particles in phase space. right: distribution at t=200.

This result shows that there is a hole structure in phase space as is the caseof classical Maxwell plasma. It indicates the existence of hole in heavy ion
collision.

17

Hole structures of classical electro-magnetic plasma are shown in recent research[8]. To, confirm whether the structure exist, multi-dimentional simulations are necessary. Plasma hole is one the significant topic for this research.

Figure 4.3: Time evolution of energy of electro-magnetic fields vs. totalenergy.

In Figure4.3, the time evolution of energy ratio of electro-magnetic fields andtotal energy is shown. At first, the energy rise quickly, and then, decrease torelaxation. Further work need to analyse the linear growth of its evolution
to make sure which instability takes place.
Figure4.4 shows the time evolution of energy ration of electro-magnetic fields
and total energy. Blue is the result with initial condition 2. Green andblack show the results with initial condition 1 and without nonlinear termsof Yang-Mills equations and with initial condition 1 respectively. In theseresults, no significant diﬀerences cannot be confirmed.

18

Figure 4.4: Time evolution of energy of electro-magnetic fields vs. totalenergy. Blue: with finite initial electromagnetic fields. Green: with zero
initial fields. Black: without nonlinear terms of Yang-Mills equations andwith zero initial fields.

19Chapter 5

Conclusion
I study the fundamental physics of classical Yang-Mills plasma dynamics.
This is the key to analyse the thermalisation process of heavy ion collisions.
I employ initial condition which has counter flow of particles. The resultsshows the existence of hole structure which is obtained in simulations ofclassical Maxwell plasma with weak initial fields.
It is reported that plasma hole can be shown not only in simulation but also
in realistic phenomena. In this thesis, the possibility that there are hole
structure in thermalising plasma in heavy ion collisions. In addition, shockwaves in collisionless plasma[9] are valuable topic to analyse. Thermalisinggluon plasma can be regarded as collisionless plasma. Therefore, there ispossibility that shock wave exist in such plasma and it can aﬀect the thermalisation process of heavy ion collisions.
For future prospectives, it is necessary to increase the accuracy and spatial
dimension of simulations for more various research including realistic initial conditions. To analyse specific phenomena such as collisionless shock
and plasma hole, we need to expand spatial region and spatial dimension ofsimulations.

20

Acknowledgements
I would like to thank my supervisor Dr.Takahiro Miyoshi for discussion andadvice of research. I acknowledge all staﬀs in Quark physics laboratory,
Prof. Toru Sugitate, Prof. Kenta Shigaki and Dr.Kensuke Homma. All
staﬀs and members of laboratory gave me a lot of important suggestions at
the meeting.

21

Bibliography[1] A. K. Chaudhuri, ”Centrality dependence of elliptic flow and QGP
viscosity”, arXiv:0910.0979v3 [nucl-th] (2010).
[2] F. Gelis, ”Color Glass Condensate and Glasma” , arXiv:1211.3327v2[hep-ph].
[3] Wong, S. K., ”Field and particle equations for the classical Yang-Mills field and particles with isotopic spin”,
https://doi.org/10.1007/BF02892134.
[4] http://www.astro.phys.s.chiba-u.ac.jp/pcans/.
[5] Dumitru, Adrian and Nara, Yasushi and Strickland, Michael, ”Ultraviolet avalanche in anisotropic non-Abelian plasmas”, 10.1103/PhysRevD.75.025016(2007).
[6] Adrian Dumitru,Yasushi Nara,Phys.Lett.B,621,(2005)89-95.
[7] Takashi Okazaki,”Plasma Heating by a Relativistic Electron Beam. II.
Two-Stream and Return Current Driven Instabilities”,Journal of the
Physical Society of Japan, 49, 1532-1541 (1980).
[8] I. H. Hutchinson, ”Electron holes in phase space: What they are andwhy they matter”,Physics of Plasmas 24, 055601 (2017).
[9] G. K. Parks, E. Lee, S. Y. Fu, N. Lin, Y. Liu, Z. W. Yang,”Shocks in collisionless plasmas”, Reviews of Modern Plasma Physics 10.1007/s41614017-0003-4(2017).

22