USING LATTICE QUANTUM CHROMODYNAMICS FOR NON-PERTURBATIVE
RENORMALIZATION OF THE GLUON MOMENTUM FRACTION OF THE
NUCLEON

By

Matthew David Zeilbeck

A THESIS

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

Physics—Master of Science

2023

ABSTRACT

Lattice quantum chromodynamics, often abbreviated as lattice QCD or LQCD, presents a

non-perturbative method for calculating phenomenological quantities in quantum field theory

(QFT). However, scale-dependent quantities such as the gluon momentum fraction

must

x
⟨

⟩g

be renormalized in order to compare lattice calculations to other theory and even other lattice

calculations. The gluon momentum fraction is useful for understanding hadron structure and

calculating the parton distribution function (PDF). I obtain the energy spectrum and ground

state matrix elements of the gluon energy-momentum tensor (EMT) from the 2pt and 3pt

correlators calculated from LQCD. Then, I calculate of the NPR renormalization constants

using the regularization-independent momentum subtraction scheme (RI/MOM) for various

lattice spacings with a pion mass Mπ

≈

310MeV. This calculation required various steps in

Chroma and Mathematica, and the workflow for this process is comprehensively detailed.

I dedicate this thesis to anyone who has ever supported me or believed in me and my
abilities and to anyone who receives any kind of benefit or positive outcome from reading it.

iii

ACKNOWLEDGEMENTS

This work would not have been possible without the assistance and support of various people.

First, I would like to thank my advisor, Dr. Huey-Wen Lin, for her support through my

struggles and shortcomings. I am especially grateful for her selflessness as she allowed me

to work on projects in her group despite knowing that I did not intend to complete my

doctorate and that my research interests lie elsewhere. I also thank Dr. Remco Zegers, the

physics graduate program director for most of my graduate education, for his support and

guidance about what I could do to achieve my goals. Both of these people supported my

desire to transfer in order to pursue research in quantum gravity and string theory, and even

though this has yet to materialize, they nonetheless did everything in their power to set me

up for success when others would have told me that this goal is too lofty and unrealistic.

Next, I thank Zhouyou Fan for allowing me to take part in the NPR research that was

ultimately his project. He (and Dr. Lin) gave me the opportunity to become part of this

project, which allowed be to get published as a co-author. Of course, this work was also the

basis for this thesis, so this thesis would not have been possible without his research acumen.

Lastly, I thank my parents for their unending support and for always being there for

me when I was down. Many times throughout my graduate education I needed to vent my

frustrations to someone, and I could always rely on them to listen. They may not have

always been able to provide the most directly useful advice, but their first and foremost

concern was my happiness. Their countless sacrifices and support is why I was able to attain

my educational success in the first place. I certainly wouldn’t be where I am today without

them.

iv

TABLE OF CONTENTS

CHAPTER 1

1.1 Particle Physics and Baryon Structure
1.2 The Path Integral Formulation of Quantum Mechanics

INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . .
. . . . . . . . . .
1.2.1 The Schrödinger Equation, Operator Formalism . . . . . . . . . . . .
1.2.2 Feynman’s Path Integral Approach . . . . . . . . . . . . . . . . . . .
1.2.3 Euclidean Time . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . .

1.2.3.1 Quantum Mechanics in Euclidean Time

1.3 Quantum Field Theory Path Integral

CHAPTER 2

QUANTUM CHROMODYNAMICS ON THE LATTICE . . . . .
2.1 Continuum QCD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.1.1 The QCD Lagrangian, Gauge Fields, and Color Components . . . . .
2.1.2 Local Gauge Invariance . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 QCD Actions on the Lattice . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.1 Lattice Gauge Invariance and Link Variables . . . . . . . . . . . . . .
2.2.1.1 Fermionic Part . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.1.2 Gauge Part . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . .
2.3.1 Calculating Lattice Quantities . . . . . . . . . . . . . . . . . . . . . .
Spectral Decomposition of Lattice Quantities . . . . . . . . . . . . . .
2.3.2

2.3 Evaluating Lattice Quantities

CHAPTER 3

ANALYZING LATTICE CORRELATOR DATA . . . . . . . . . .
3.1 The 2pt and 3pt Correlators . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . .
3.2 Fitting the 2pt Correlator
. . . . . . . . . . . . . . . . . . . . . . . . . .
3.3 Fitting the 3pt Correlator

CHAPTER 4

NONPERTURBATIVE RENORMALIZATION (NPR)

. . . . . .

4.1 Regularization-independent Momentum Subtraction (RI/MOM)

Scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2 Cluster-decomposition Error Reduction (CDER) . . . . . . . . . . . . . .
4.3 Renormalized Gluon Momentum Fraction . . . . . . . . . . . . . . . . . .
4.4 NPR Calculation Workflow . . . . . . . . . . . . . . . . . . . . . . . . . .
4.4.1 Chroma . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Implementation of the NPR Calculation . . . . . . . . . . .
4.4.1.1
4.4.1.2 Outstanding Issues . . . . . . . . . . . . . . . . . . . . . . .
4.4.2 Mathematica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.4.2.1 My Contribution to the Notebook . . . . . . . . . . . . . . .

CHAPTER 5

CONCLUSION . . . . . . . . . . . . . . . . . . . . . . . . . . . .

REFERENCES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

v

1
1
2
2
3
4
5
7

8
8
8
9
10
10
10
12
13
13
15

17
17
18
21

26

27
28
30
32
33
38
41
41
45

46

47

CHAPTER 1

INTRODUCTION

1.1 Particle Physics and Baryon Structure

All matter in the universe is comprised of particles called leptons and hadrons. Leptons

have no further internal composite structure — they are fundamental — while hadrons are

made of quarks and are thus composite particles. Hadrons are further classified into mesons

and baryons. Mesons are lighter and are made of a quark-antiquark pair while baryons are

heavier and made of three quarks. In order for quarks to bind together to make a hadron,

they must interact with each other and form a bound state. The nature of these interactions

is quite complicated. For example, a proton has only three valence quarks (up, up, down

or uud) that are always detectable, but inside the proton, there is a complicated web of

interactions involving gluons and virtual quark-antiquark pairs that affect the nature of the

proton, e.g. see. Fig. 1.1. It used to be believed that the properties of hadrons like the

proton were almost entirely determined by their valence quarks . However, we now know

that this is not the case; it is as a result of both the properties of quarks and the nature of

the interactions between them that hadrons acquire various physical properties such as spin,

magnetic moment, observed mass1, effective size, etc. In addition, the gluons can carry a

nontrivial fraction of the total momentum of a hadron, and calculating this fraction is part

of understanding hadron structure.

In addition to furthering our understanding of hadron structure, the gluon momentum

fraction is also useful for understanding the gluon parton distribution function (PDF). The

PDF g(x) is a probability density for finding a parton with momentum fraction x in a hadron.

The gluon momentum fraction

x
⟨

⟩g

is the first moment of the PDF,

(cid:90) 1

x
⟨

⟩g =

0

xg(x) dx .

(1.1)

The gluon momentum fraction serves another purpose besides the study of the PDF mo-

1It is primarily the energy of the interactions with gluons that we observe as the a hadron’s mass, not

the sum of the bare quark masses.

1

Figure 1.1: Comparison of our understanding of the composition of the proton over time.
The left image shows the previous model of the proton as consisting only of its valence
quarks since gluon effects were not thought to be important. The right image shows the
current understanding of the dynamics of the proton with valence quarks, virtual quarks,
and gluons. The larger spheres that represent the valence quarks are made larger for
emphasis and do not reflect the relative sizes of particles. The coil-like lines represent
gluons that interact with the quarks to which they are attached. Image source [1].

ments. PDFs are of interest in other parts of high-energy physics, so it is worthwhile to

calculate them in LQCD. One method for calculating the PDF is the “pseudo-PDF" ap-

proach [2, 3, 4]. This method does not directly obtain g(x) but instead yields the quantity

xg(x)/

x
⟨

⟩g

. Consequently,

x
⟨

⟩g

must be known in order to eliminate it from this quantity.

For these reasons, one goal of this work is to calculate the gluon momentum fraction of the

nucleon.

1.2 The Path Integral Formulation of Quantum Mechanics

1.2.1 The Schrödinger Equation, Operator Formalism

Lattice QCD developed from quantum field theory, which itself was developed from quan-

tum mechanics. In non-relativistic quantum mechanics, one can solve physical problems by

2

obtaining the wavefunction of a system which obeys the Schrödinger equation

(cid:32)

N
(cid:88)

ℏ2

−

i=1

2mi ∇

(cid:33)

2
i + V (xi, t)

Ψ(xi, t) = iℏ ∂
∂t

Ψ(xi, t),

(1.2)

where mi is the mass of the ith particle, V (xi, t) is the potential, Ψ(xi, t) is the wavefunc-

tion, and the sum is over N particles. xi is shorthand for the positions of each of the N

particles

xi
{

}

. The wavefunction is overlap of the time-dependent state

with the po-

Ψ(t)
⟩
|

sition eigenstate

xi

|

⟨

such that Ψ(xi, t) =

xi
⟨

. Then we evaluate expectation values of
Ψ(t)
⟩
|

physical quantities with

(t) =

O

⟨

⟩

ˆO(ˆx, ˆp)
Ψ(t)
Ψ(t)
|
|
⟨

⟩

=

(cid:90)

Ψ∗(x, t)O(x,

iℏ

∇

−

)Ψ(x, t) d3x

(1.3)

where ˆO(ˆx, ˆp) is the operator associated with the observable

O
⟨

, ˆp =
⟩

iℏ

−

∇

in the position

basis, and the hat is used to denote an operator. All operators are functions of the canonical

variables (ˆx, ˆp) in operator form.

1.2.2 Feynman’s Path Integral Approach

When Richard Feynman was in graduate school, he noticed a remark in Paul Dirac’s

book that states that exp

(cid:20)

(cid:90) t2

(cid:21)

i

S

“corresponds to"

x2, t2|
⟨

x1, t1⟩

where

t1

(cid:90) t2

t1

S =

L(x, ˙x)
ℏ

, L(x, ˙x) = T ( ˙x)

V (x).

−

In this expression, T is the kinetic energy, V is the potential energy, L is the classi-

cal Lagrangian, and S is the classical action. By inserting completeness relations ˆI =
(cid:90)

x, t
|

⟩⟨

x, t
|

dx and splitting up the time interval t2 −

t1 into infinitesimal intervals, it can

be shown that the probability amplitude of a particle located at initial position xi at time

ti being found at some other final location xf at time tf is

xf , tf
⟨

xi, ti
|

⟩

= lim
N →∞

(cid:19)N/2 (cid:90)

(cid:18) m
2πi

tf

N

−

ti

dx1 dx2 · · ·

dxN −1 eiS.

(1.4)

This expression is for 1 dimension, but it can easily be generalized to 3 dimensions to obtain

xf , tf
⟨

xi, ti
|

⟩

= lim
N →∞

(cid:19)3N/2 (cid:90)

(cid:18) m
2πi

tf

N

−

ti

3

d3x1 d3x2 · · ·

d3xN −1 eiS.

(1.5)

This expression becomes exact when N goes to infinity, so it may seem like the path integral

is impossible and impractical to evaluate. However, this result generalizes very nicely to

lattice formalism where spacetime is discrete and finite. Thus, path integrals of the form of

Eq. 1.5 can be evaluated on a lattice where each integral over d3xi corresponds to a time slice

ti of the Euclidean lattice. Since there is now a finite number of integrals, lattice formulation

of Eq. 1.5 will have discretization errors but becomes more accurate as the lattice resolution

becomes sharper and more continuum-like.

1.2.3 Euclidean Time

Lattice theory is often formulated using an imaginary time t =

iτ , τ > 0. This is

−

called Wick rotation. This means that we use a Euclidean metric gµν = δµν where δµν is

the Kronecker delta, and the scalar product of two 4-vectors xµ = (

−

i˜x0, x), yµ = (

−

i˜y0, y)

(upper vs.

lower indices no longer matter with a Euclidean metric), with ˜x0, ˜y0 being the

analog to τ , are now represented as xµgµνyν. When t is replaced with

iτ and time derivates

−

are understood to be derivatives with respect to τ , the action S becomes

(cid:90)

(cid:90)

S =

(T + V ) dt = i

HE dt = iSE, SE

(cid:90)

≡

HE dt

(1.6)

where HE =

1
2m

2

(cid:12)
(cid:12)
(cid:12)
(cid:12)

dx
dτ

(cid:12)
(cid:12)
(cid:12)
(cid:12)

+ V (x) is the Euclidean Hamiltonian. Now the path integral becomes

xf , τf
⟨

xi, τi
|

⟩

= lim
N →∞

(cid:19)3N/2 (cid:90)

(cid:18) m
2π

τf

N

−

τi

d3x1 d3x2 · · ·

d3xN −1 e−SE .

(1.7)

In classical mechanics, allowed trajectories x(t) satisfy the least action principle, δS = 0.

That is, they minimize the action S. In contrast, in quantum mechanics, all paths contribute

to the path integral. However, eiS in Eq. 1.5 rapidly oscillates for paths with large S such that

paths far from the classically allowed trajectory do not contribute much to the integral. That

is, the rapidly oscillating term causes a lot of cancellation. In the Euclidean path integral

Eq. 1.7, the oscillating eiS is now an exponentially decaying e−SE , so that rapid oscillations

are now replaced with exponential suppression for non-classical paths in the integral where

the amount of suppression increases with deviation from the classical trajectory.

4

Although I used the substitution t =

iτ , in practice, t is used in LQCD equations, and

−

it should be understood to be the Euclidean time.

1.2.3.1 Quantum Mechanics in Euclidean Time

Changing to Euclidean time affects other equations and tensorial objects in quantum

mechanics in order to preserve various expressions involving 4-vectors. For instance, the

Dirac equation’s gamma matrices in Minkowski spacetime γM
µ

satisfy

(cid:8)γM

µ , γM
ν

(cid:9) = 2ηµν

(1.8)

where ηµν = diag(1,

1,

1,

−

−

1) is the Minkowski metric. When we use Euclidean time, the

−

metric should be replaced with the Euclidean metric so that

γµ, γν

{

}

= 2δµν.

(1.9)

Thus, we must modify the γ matrices from their typical Minkowski definitions and define

the Euclidean γµ. The anticommutation relation will hold if we define the Euclidean gamma

matrices as

−
where i denotes a spatial index. Note that while 0 is used as the temporal index in Minkowski

γ4 = γM
0 ,

γi =

iγM
i

, i = 1, 2, 3

(1.10)

spacetime, 4 is customarily used in Euclidean spacetime. Any γµ should be assumed to be

defined with the Euclidean definition unless otherwise specified.

This definition of the gamma matrices modifies the appearance of the Dirac equation as

well. Normally the Dirac equation reads

(cid:0)pµγM

µ −

m(cid:1) ψ(xµ) = 0

where pµ represents (p0, p) and γM
µ

3 ). Summations over repeated
indices (Einstein summation convention) are implied throughout. In Euclidean spacetime,

represents (γM

1 , γM

0 , γM

2 , γM

4-vectors are denoted as

vµ = vµ = (v1, v2, v3, v4) = (v1, v2, v3, iv0)

5

where the distinction between contravariant (upper) and covariant (lower) indices does not

matter, and the components of each 4-vector are equal to their contravariant Minkowski

counterparts (except that the time component is imaginary). Therefore, we can rewrite the

operator in the Dirac equation as

pµγM

µ −

m =

−

i (cid:0)p1γ1 + p2γ2 + p3γ3(cid:1)

ip4γ4

−

m =

ipµγµ

m.

−

−

−

Since the Euclidean metric makes no distinction between upper and lower indices, both

indices in the implied summation can be written lower as they are in the above equation.

However, Greek indices are still used to indicate that the sum is over 1–4 rather than just the

spatial indices. However, Euclidean time treats time symmetrically with respect to space.

As a consequence, the Dirac equation now reads

(

ipµγµ

−

−

m) ψ(xµ) = 0.

(1.11)

Similarly, the completeness relations for the Dirac spinors are also modified, now reading

(cid:88)

s=1,2

u(s)u(s) = N (

ipµγµ + m) ,

−

(cid:88)

s=1,2

v(s)v(s) = N (

−

ipµγµ

m) ,

−

(1.12)

where N is a normalization constant that depends on the normalization used for the bispinors.

u(s) are the regular solutions to the Dirac equation (

ipµγµ

m) u(s)(xµ) = 0, v(s) are the

−
ipµγµ + m) v(s)(xµ) = 0, and the adjoint bispinor u(s) is

−

negative energy solutions to (
u(s) = u(s)†

γ4.

−

The above definition for the Euclidean γ matrices is not unique. For instance, one can

define

but this also requires vi =

−

Then, the Dirac equation is (ipµγµ

−

γ4 = γM
4 ,

γi = iγM
i

,

vM,i such that the Euclidean 4-vector is vµ = (

vM ,

−
m) ψ(xµ) = 0, and the completeness relations

v0).

−

(Eq. 1.12) are similarly modified by eliminating the negative sign in

ipµγµ.

−

6

1.3 Quantum Field Theory Path Integral

Quantum field theory (QFT) is the theory that comprehensively describes relatistic quan-

tum mechanics and the fundamental interactions, excluding gravity. In classical and quantum

mechanics, the primary degrees of freedom are the coordinates x, y, z. In a field theory, the

degrees of freedom are the value of the field at each point in spacetime. Since spacetime is

continuous, this results in infinite degrees of freedom. The Lagrangian formalism of classical

mechanics can be generalized to QFT with the goal of finding a Lagrangian density

that the Euler-Lagrange equations,

∂
L
∂ϕ −

∂µ

∂
L
∂(∂µϕ)

= 0,

such

L

(1.13)

yield the differential equation that the field satisfies, e.g. the Klein-Gordon equation, Dirac

equation, etc. The path integral can be generalized to field theory by integrating over every

value of the field at every value of space in addition to every value of time:

ϕ′(xµ)
ϕ(yµ)
⟨
⟩
|
(cid:90)

=

(cid:90)

(cid:90)

[ϕ]eiS,

D

[ϕ]

D

≡

(cid:90)

lim
N →∞
(cid:90) (cid:89)

ϕ1 · · · D

D

ϕN −1,

dϕ(x, ti) ,

(1.14)

ϕi

D
(cid:90)

≡

x

S =

L

(ϕ, ∂µϕ) d4x .

7

CHAPTER 2

QUANTUM CHROMODYNAMICS ON THE LATTICE

2.1 Continuum QCD

2.1.1 The QCD Lagrangian, Gauge Fields, and Color Components

The theory of the strong interaction is quantum chromodynamics (QCD). It is mediated

by gluon fields which are described by the potential Aµ(x)cd, which is a 4-vector (with Lorentz

index µ) of 3

×

3 hermitian, traceless matrices (color indices c and d). These matrices are

elements of the Lie algebra su(3). Together with the bispinor fields ψ(x)α
c

(here α is a Dirac

index), the fermionic part of the QCD Lagrangian density can be written as

(cid:2)ψ, ψ, A(cid:3) =

F

L

Nf
(cid:88)

f =1

(f )

ψ

(x)α
c

(cid:0)(γµ)αβDµ(x)cd + m(f )δαβδcd

(cid:1) ψ(f )(x)β

,

d

where

Dµ(x)cd = δcd∂µ + igAµ(x)cd

(2.1)

(2.2)

is the covariant derivative, g is the coupling strength between the fermionic and gauge fields,

f is the flavor index, and Nf is the total number of flavors in consideration. The gauge field

part of the Lagrangian density is

G[A] =

L

1
2

Tr[Fµν(x)Fµν(x)].

(2.3)

The field strength tensor Fµν(x) is defined as

Fµν(x)

i[Dµ(x), Dν(x)] = ∂µAν(x)

≡ −

−

∂νAµ(x) + i[Aµ(x), Aν(x)]

(2.4)

where the color indices have been suppressed. Since Aµ(x)

su(3), we can decompose Aµ(x)

∈

into a superposition of the generators of SU(3) (which are the elements of the algebra su(3)):

Aµ(x) =

8
(cid:88)

i=1

A(i)

µ (x)Ti, Ti =

1
2

λi,

(2.5)

where λi are the Gell-Mann matrices. This first summation is explicitly written despite the

repeated index to emphasize that this sum runs from 1 to 8 (over color components). The

8

4-vectors A(i)

µ (x) are the color components of the gauge field. Thus, the field strength tensor

can be written as

Fµν(x) =

8
(cid:88)

i=1

F (i)

µν (x)Ti

µν (x) = ∂µA(i)
F (i)

ν (x)

∂νA(i)

µ (x)

−

−

fijkA(j)

µ (x)A(k)

ν (x),

(2.6)

where the structure constants fijk are defined by the relation [Ti, Tj] = ifijkTk. Thus, the

gauge field part of the Lagrangian density can be written as

G[A] =

L

1
4

µν (x)F (i)
F (i)

µν (x),

where the trace of products of Ti is evaluated with the relation Tr[TiTj] =

QCD action is thus

S (cid:2)ψ, ψ, A(cid:3) =

(cid:90)

d4x (cid:0)
L

F [ψ, ψ, A] +

G[A](cid:1) .

L

(2.7)

δij. The full

1
2

(2.8)

2.1.2 Local Gauge Invariance

The Lagrangian is defined in the form of Eqs. 2.1 and 2.3 so that it is invariant under

local SU(3) gauge transformations. These transformations modify the fields as

ψ′(x) = Ω(x)ψ(x), ψ

′

(x) = ψ(x)Ω(x)†, A′

µ(x) = Ω(x)Aµ(x)Ω(x)† + i (∂µΩ(x)) Ω(x)†,

(2.9)

for Ω(x)

∈

SU(3). The transformation properties of Aµ(x) are defined like this so that the

covariant derivative transforms as

D′

µ(x) = Ω(x)Dµ(x)Ω(x)†.

Because Fµν(x) is the commutator of covariant derivatives, it transforms as

µν(x) = Ω(x)Fµν(x)Ω(x)†.
F ′

This ensures that S

(cid:104)
ψ′, ψ

′

, A′(cid:105)

= S (cid:2)ψ, ψ, A(cid:3). That is, the action is gauge-invariant.

(2.10)

(2.11)

9

2.2 QCD Actions on the Lattice

The following subsection is largely based on [5], and the notation used there largely

matches the notation used here.

2.2.1 Lattice Gauge Invariance and Link Variables

2.2.1.1 Fermionic Part

The Lagrangian density and action introduced in the previous section apply to the con-

tinuum, but the expressions in Eqs. 2.1, 2.3, and 2.8 are not suitable for the lattice in their

current form. This is because of the derivative term which must be discretized on the lat-

tice. Let x = an where a is the lattice spacing and n = (n, n4) is a vector of integers which

denotes which lattice site is referenced, e.g. ni

1
}
NT lattice. From hereon, ψ(n) is understood to mean the value of the field ψ at

0, 1, . . . , NT

0, 1, . . . N

1
}

∈ {

−

−

, n4 ∈ {

on an N 3

×

position an. Thus, the gradient on the lattice is given by

∂µψ(x) =

ψ(n + ˆµ)

ψ(n

ˆµ)

−

−
2a

(a),

+

O

(2.12)

where ˆµ is a vector of length a in the µ-direction. Then the discretized version of the

fermionic part of the action for a free fermion becomes

(cid:2)ψ, ψ(cid:3) =

SF

(cid:90)

(cid:2)ψ, ψ(cid:3) d4x

F

L

a4 (cid:88)

n∈Λ

→

(cid:18)

ψ(n)

γµ

ψ(n + ˆµ)

ψ(n

ˆµ)

−

−
2a

(cid:19)

+ mψ(n)

,

where Λ is the set of all spacetime points on the lattice, and the integral is replaced with a

discrete sum. This presents a problem for gauge invariance because when ψ is replaced with

ψ′, we are left with terms like

′

ψ

(n)ψ′(n

±

ˆµ) = ψ(n)Ω(n)†Ω(n

ˆµ)ψ(n

±

±

ˆµ)

= ψ(n)ψ(n

ˆµ)

±

that spoil the gauge invariance that we had in the continuum theory. The solution is to

introduce a new field Uµ(n) that transforms under a gauge transformation as

U ′

µ(n) = Ω(n)Uµ(n)Ω(n + ˆµ)†.

(2.13)

10

̸
These are called link variables because they link site n with the neighboring site in the ˆµ-

direction (site n + ˆµ) through this gauge transformation property. A negative direction may

also be specified, and Uµ(n) is defined such that

U−µ(n)

Uµ(n

ˆµ)†,

−

≡

(2.14)

as shown in Fig. 2.1.

Figure 2.1: The link variable Uµ(n) (right) links lattice site n with site n + ˆµ and is
directed in the ˆµ-direction. The link variable U−µ(n) (left) links n with n
ˆµ and is
directed in the
−
also displayed. Image source [5].

ˆµ-direction. The definition the hermitian conjuage of the link variables is

−

Now we note that

′

ψ

(n)U ′

±µ(n)ψ′(n

±

ˆµ) = ψ(n)U±µ(n)ψ(n

ˆµ)

±

is a gauge-invariant construct and try as a candidate for the fermionic action

(cid:2)ψ, ψ, U (cid:3) = a4 (cid:88)

SF

(cid:18)

ψ(n)

γµ

Uµ(n)ψ(n + ˆµ)

n∈Λ

U−µ(n)ψ(n

−
2a

ˆµ)

−

(cid:19)

+ mψ(n)

.

(2.15)

Although this is gauge invariant, we must still determine what U could be. It turns out that

the path-ordered exponential

(cid:18)

(cid:90)

(cid:19)

P exp

i

Aµ dsµ

,

C

along some contour

C

, has the gauge transformation property required by Eq. 2.13. Here

the coupling g has been absorbed into the field Aµ. Since Uµ connects adjacent lattice sites

that are a small distance apart, the integral can be approximated such that

Uµ(n) = exp (iaAµ(n)) .

(2.16)

11

To see that Eq. 2.15 with U instead of A is correct, we see that when Eq. 2.16 is expanded

through order

O

(a), the action becomes

(cid:2)ψ, ψ, A(cid:3) = a4 (cid:88)

SF

(cid:20)

ψ(n)

γµ

(cid:18) ψ(n + ˆµ)

n∈Λ

ψ(n

ˆµ)

−

−
2a

+ iAµ(n)ψ(n) +

O

(cid:19)

(cid:21)

(a)

+ mψ(n)

,

(2.17)

which matches Eq. 2.1 if one attempts to naively discretize that expression to use on the

lattice. Thus, it is clear that the correct gauge degrees of freedom on the lattice are Uµ(n)

instead of Aµ(n).

2.2.1.2 Gauge Part

Since the correct gauge field to use on the lattice is U , the expression for the gauge part

of the action SG[A] must also be modified. The trace of a product of link variables that

trace out a closed contour is gauge-invariant, so this is the type of object that will be used

to construct the action. The simplest contour of this type is the plaquette as illustrated in

Fig. 2.2. It is a closed contour made of 4 link variables that make a planar square.

Figure 2.2: A counter-clockwise plaquette in the ˆµ- and ˆν-directions. The links that appear
to go against the orientation of the loop (Uµ(n + ˆν) at the top and Uν(n) to the left) are
actually written to go in the positive ˆµ- and ˆν-directions because the diagram is meant to
show the relevant link variables written in standard form such that they point in the
positive directions. Thus, the link variables that appear in the plaquette are Uµ(n + ˆν)†
and Uν(n)†. Image from [5].

The plaquette specifies two directions µ and ν for the directions of the links. If we wish

12

to traverse in the µ-direction first, then we write the plaquette as

Uµν(n) = Uµ(n)Uν(n+ ˆµ)U−µ(n+ ˆµ+ˆν)U−ν(n+ˆν) = Uµ(n)Uν(n+ ˆµ)Uµ(n+ˆν)†Uν(n)†. (2.18)

The Wilson gauge action is constructed by summing over all plaquettes in a single orientation

on the lattice:

SG[U ] =

2
g2

(cid:88)

(cid:88)

n∈Λ

µ<ν

Re Tr [1

−

Uµν(n)] .

(2.19)

The second summation requires µ < ν because when µ and ν are swapped, the plaquette

reverses orientation, and if µ = ν, then the plaquette has nonsensical meaning (actually it

does not contribute to the summation because Uµµ = 1). If Uµν(n) is Taylor expanded and

rewritten with the Baker-Campbell-Hausdorff formula up to the first term with a commuta-

tor, then to lowest non-vanishing order in Eq. 2.19, we find that

SG[U ] =

2
g2

(cid:88)

(cid:88)

n∈Λ

µ<ν

Re Tr [1

−

Uµν(n)] =

a4
2g2

(cid:88)

(cid:88)

(cid:0)Tr (cid:2)Fµν(n)2(cid:3) +

n∈Λ

µ,ν

(cid:0)a2(cid:1)(cid:1) .

O

This matches Eq. 2.3 to lowest order, except once again, A has been rescaled to absorb a

factor of g such that 1/g2 appears in the expression. The sum over µ < ν has been replaced

with an unconstrained sum over µ and ν with a factor of 1/2 to account for the fact that

this double-counts terms (also the sum over µ, ν is explicitly written out instead of implied

like in Eq. 2.3). Note that all of these traces are over color indices only.

2.3 Evaluating Lattice Quantities

2.3.1 Calculating Lattice Quantities

Now we can write the entire action as

S (cid:2)ψ, ψ, U (cid:3) = SF

(cid:2)ψ, ψ, U (cid:3) + SG[U ],

and we can use this to evaluate the path integral on the lattice.

In Eq. 1.7, there are

similarities to the types of integrals that one performs in statistical physics calculations

where the Boltzmann factor is e−SE . In addition, the amplitude

xf , T
⟨

xi, 0
⟩
|

can be written

as

xf
⟨

e−T ˆH
|

xi
|

⟩

where e−T ˆH looks like the canonical density operator ˆρ = e−β ˆH where β

T .

→

13

That is, the size of the time evolution/the temporal extent of the lattice acts as an inverse

temperature in the Euclidean LQCD theory. This suggests defining the partition function

as

ZT = Tr

(cid:104)

e−T ˆH(cid:105)

=

(cid:90)

ψ0D

ψ0D

D

U0

(cid:10)ψ0, ψ0, U0

(cid:12)e−T ˆH(cid:12)
(cid:12)

(cid:12)ψ0, ψ0, U0

(cid:11) .

(2.20)

However, the matrix element in the integrand can be evaluated with the path integral for-

malism if we generalize Eq. 1.7 to fields. Then the partition function on the lattice becomes

(cid:90)

ZT =

(cid:2)ψ, ψ(cid:3)

D

D

[U ]e−S[ψ,ψ,U],

(cid:2)ψ, ψ(cid:3) =

(cid:89)

(cid:89)

n∈Λ

f,α,c

D

dψ(f )(n)α
c

dψ

(f )

(n)α
c

,

[U ] =

D

(cid:89)

4
(cid:89)

n∈Λ

µ=1

dUµ(x) .

(2.21)

Overall multiplicative factors that appear outside of the integral in equations like Eq. 1.7 are

omitted because they cancel out in relevant calculations. This motivates defining expectation

values as

O2(t)O1(0)

⟩ ≡

⟨

1
ZT

(cid:104)

(cid:105)
ˆρ ˆO2(t) ˆO1(0)

Tr

=

1
ZT

(cid:104)
e− ˆH(T −t) ˆO2e− ˆHt ˆO1

(cid:105)

,

Tr

(2.22)

where in the second equality, the Heisenberg representation time-evolution of operators is

used, i.e.

where ˆO

≡

ˆO(t) = e ˆHt ˆO(0)e− ˆHt,

(2.23)

ˆO(0) and the exponential is real because of Euclidean time. This particular

example is a 2pt correlation function, but it can be generalized to simpler or more complex

observables as well. It will be assumed that the operators ˆO are functions of the relevant

fields (and not any conjugate fields, e.g. conjugate momentum), i.e. ˆO = ˆO

(cid:104) ˆψ, ˆψ, ˆU

(cid:105)
. The

same methodology used to derive Eqs. 1.5, 1.7, and 2.21 shows that Eq. 2.22 can also be

evaluated in the path integral approach and becomes

O2(t)O1(0)
⟩
⟨

=

(cid:90)

(cid:8)

1
ZT

(cid:2)ψ, ψ(cid:3)

D

D

[U ]O2

(cid:2)ψ(n, nt), ψ(n, nt), U (n, nt)(cid:3)

(cid:2)ψ(n, 0), ψ(n, 0), U (n, 0)(cid:3) e−S[ψ,ψ,U](cid:111)

O1

×

(2.24)

14

Because these integrals are identical in form to the integrals used to evaluate quantities

in statistical physics, the same numerical techniques can be used here to calculate lattice

quantities as well. However, the generation of lattice configurations and the calculation of

lattice quantities from actual lattice configurations are outside the scope of this thesis.

2.3.2 Spectral Decomposition of Lattice Quantities

Once an observable has been calculated from Eq. 2.24, it is analyzed through a spectral

decomposition of the observable and the partition function. Now it is better to evaluate the

trace using a complete basis of energy eigenstates in addition to the identity operator as

the sum over projections to energy eigenstates I =

becomes

(cid:88)
n
n |

.1 Thus, the partition function
n
|

⟩⟨

(cid:88)

ZT =

e−EnT

1 for sufficiently large T,

(2.25)

≈
where En+1 > En. Then Eq. 2.24 becomes

n

O2(t)O1(0)
⟨

⟩ ≈

(cid:88)
m,n ⟨

n
O2|
m
|

⟩ ⟨

n

O1|
|

m
⟩

e−Em(T −t)e−Ent

(cid:88)
n
O2|
0
|
n ⟨

⟩ ⟨

≈

O1|
n
0
⟩
|

e−Ent,

(2.26)

where in the last approximation it is once again assumed that T is sufficiently large.

Now that the groundwork for lattice theory has been laid down, the rest of this thesis

is organized as follows. Chapter 3 explains how lattice correlator data is analyzed in order

to extract physical observables.

I begin with the definitions of the 2pt and 3pt lattice

correlators. Then I explain how the correlator data is fitted and how useful quantities such

as the matrix element associated with the gluon momentum fraction are extracted from the

fit (Sec. 3.1). Next, in Chapter 4, I explain how to calculate the renormalization constants

to perform non-perturbative renormalization (NPR) of the gluon energy-momentum tensor

(EMT) in order to renormalize the gluon momentum fraction calculated in Chapter 3.

I

explain what quantities are needed for the calculation and how to calculate different parts

1This definition can be modified depending on normalization conventions, e.g. I =

(cid:88)

Gordon fields with Lorentz-invariant normalization, i.e.

p′
p
|

⟨

⟩

= 2E(
|

)(2π)3δ3 (p
p
|

−

n
p′).

|

n
n
⟩⟨
2En

|

for Klein-

15

of the renormalization constant using the RI/MOM scheme (Secs. 4.1–4.3), including the

CDER technique to improve the signal (Sec. 4.2). I also explain the workflow for the different

steps required to turn the lattices into results (Sec. 4.4). Finally, I present a conclusion and

summary of my work in Chapter 5.

16

CHAPTER 3

ANALYZING LATTICE CORRELATOR DATA

3.1 The 2pt and 3pt Correlators

The particular lattice quantities required for the calculation of the gluon momentum frac-

tion of the nucleon

x
⟨

⟩g

are the 2pt and 3pt correlation functions (the correlation functions

are also simply referred to as correlators). For this particular project, I analyzed spin-1/2

particles that obey the Dirac equation, and the fermion fields were normalized such that

(cid:88)

I =

ns

MN
EN |

ns

ns

, where the sum is over spin s as well as energy eigenstates n. The 2pt
|

⟩⟨

correlators are defined very similarly to Eq. 2.24 as

C 2pt(p, t) =

χB(p, t)χB(p, 0)
⟨

⟩

=

(cid:88)

n

p
Zn(
|

)
|

χB(p)
n
0
|
|
⟨

⟩ ⟨

n

χB(p)
0
⟩
|
|

e−Ent,

(3.1)

where χB(p, 0) is an operator that creates baryon B with momentum p at time 0 and χB(p, t)
annihilates it at time t. Usually the baryon is instead created at some position on the lattice

(typically the origin) and annihilated after undergoing some displacement x at time t, then

Fourier transformed to the momentum of interest.

The 3pt correlators are defined as

C 3pt (pf , pi, tsep, tins) =

=

χB(pf )V (q)χB(pi)
⟨
(cid:88)

pi
Zm(
|

pf
)Zn(
|
|

)
|

⟩

m,n

χB(pf )
n
0
|
|
⟨

⟩ ⟨

n

V (q)
m
|
|

⟩ ⟨

χB(pi)
m
0
|
⟩
|

(3.2)

e−En(tsep−tins)e−Emtins.

×

This corresponds to a particle created at time 0 with 3-momentum pi at time 0 that interacts

with operator V (q = pf

−

pi) at time tins (insertion time) that imparts momentum transfer

q that is then annihilated at time tsep (separation time) with 3-momentum pf . tsep is also

known as the source-sink separation, and tsep is also sometimes called tsnk. The V operator

corresponds to some property of the particle that we want to study. As a result, the 3pt

correlator is only defined for 0

tins ≤

≤

tsep since the particle cannot be inspected after it

is annihilated. The Z terms in Eqs. 3.1 and 3.2 are kinematic factors that depend on the

17

normalization, e.g. Zn(
p
|

p
) = M/En(
|
|

) for spin-1/2 particles. For spin-1/2 particles, the
|

creation and annihilation operators have Dirac indices (not shown in Eqs. 3.1 and 3.2), so

the expressions in Eqs. 3.1 and 3.2 must be contracted by a projection operator to select out

the relevant spin projections. The details of the projector depend on the quantity that we

are trying to calculate by using the correlators.

In order to calculate the gluon momentum fraction of the nucleon, the necessary operator

V is the gluon energy-momentum tensor, and in particular, we need

V
0
|

⟨

, the ground
0
⟩
|

state diagonal matrix element. This can be extracted from the 3pt correlator by fitting

sampled 3pt correlator data and using fitted sampled 2pt correlator data to divide out the

unwanted terms and kinematic factors in Eq. 3.2.

3.2 Fitting the 2pt Correlator

The following sections involve analyzing data that comes from lattice simulations of QCD,

i.e. the 2pt and 3pt correlators. For my analysis, these were provided to me to analyze; I

was not involved in taking the data.

In addition, any values of and variables related to

time (t, tins, or tsep) that appear in the following sections should be interpreted to be integer

values, i.e. they correspond to lattice grid indices. This means that the exponential terms

in equations like Eqs. 3.1 and 3.2 contain the unitless energy aEn and the unitless integer

time t as the lattice spacing a is attached to the energy rather than the time. However, the

energy will still be written as En even though it is implied to be the unitless aEn. This is

reflected by the lack of units in tables containing energies.

The first step to analyzing correlator data is performing 2pt correlator fits. Given Ncfg

samples of the 2pt correlator as a function of time from lattice calculations, the first step is

jackknife resampling of the Ncfg configurations. Jackknife resampling is important because

of the nonlinearity of the quantities calculated with the correlators. By using the raw data

without resampling, the sample mean of such quantities has a bias that goes as O(1), i.e.

increasing the number of samples does not improve the convergence to the true mean. With

jackknife resampling, the bias goes as O (1/Ncfg), so the sample mean will converge to the

18

true mean with more data [6].

For the purposes of curve fitting, Eq. 3.1 is written as

C 2pt(t) =

(cid:88)
n |

An

2e−Ent =
|

2e−E0t +

A0|
|

2e−E1t + . . .

A1|
|

(3.3)

such that the specific matrix elements are ignored for fitting purposes and normalization

or kinematic factors are absorbed into the An. The fit parameters are then

An
|

2 and En.
|

Although it is not written for brevity, the 2pt correlator and its fit parameters are functions

of momentum. Thus, different sets of data for different momenta are required. The number

of states/terms included in the fit depends on how many can be resolved in the data, which

can be discerned by the quality of attempted multistate fits.

I present an example of 2pt correlator (and subsequently 3pt correlator analysis) using the

data from the lattice ensembles that are used for non-perturbative renormalization (NPR)

in Chapter 4, where the ensembles are explained in more detail (see Table 4.1). The analysis

in this chapter uses data from the a12m310 ensemble, and it serves to double-check the

calculations in [7]. For consistency with [7], I used a 2-state fit for this data. However,

not all of the data is used for fitting. The lowest values of t have too much excited state

contamination to be used, and the highest values of t end up having too much noise compared

to signal. Therefore, a proper subset of the range [0, Lt] is used, and selecting it is non-trivial.

In order to determine a good t range to use, I start by making a so-called effective mass 1

plot. The effective mass is given by ln [C 2pt(t)/C 2pt(t + 1)]. As t increases, Eq. 3.3 becomes

dominated by just the ground state term, the effective mass asymptotically approaches E0,

and the plot should approach a horizontal line for this range of t

0. This is visible in

≫

Fig. 3.1 where the effective mass is plotted with the reconstructed fit band obtained by using

the fitted parameters in the fitted equation. This was performed for 2pt correlators projected

to pz

[0, 5]

×

∈

2π/(aL).

1Effective mass is a bit of a misnomer because it only corresponds to the mass if the 2pt correlator is
projected to zero momentum; otherwise it corresponds to the ground state energy, which will have a kinetic
contribution. However, I’ve never heard anyone call it a ground state energy plot.

19

∈

[0, 5] in units of 2π/(aL). Fits were performed with t

Figure 3.1: Effective mass plots and reconstructed fit bands for the a12m310 ensemble for
[2, 11]. No upper bound on the
pz
fit bands was imposed (they exceed t = 11) to demonstrate that the effective mass
approaches E0. The y-axis range was chosen so that the most important features of the
data can be better visualized, but this caused some of the data points for the highest t to
not be visible for pz = 4, 5.

∈

The reconstructed fit bands in Fig. 3.1 are created by fitting the 2pt correlator to Eq. 3.3

and then calculating the effective mass with ln [C 2pt(t)/C 2pt(t + dt)] /dt with the fitted pa-

rameters. Note that dt replaces 1 (and there is division by dt) because the fitted band uses

many interpolated values at t that have smaller spacing than that of the data points. Since

I used a 2 state fit, the fitting function was the first two terms of Eq. 3.4 (up to n = 1):

C 2pt(t) =

2e−E0t +

A0|
|

2e−E1t.

A1|
|

(3.4)

In order to determine a suitable fit range, I examined the effective mass plots (without

the fit bands because they are calculated from the fit which had yet to be performed) and

picked t bounds such the fit range includes some of the data points before they approach a

horizontal line and most of the data after this occurs. However, not all of the available data

in the asymptotic region is used because eventually the data starts to become noisier and also

deviates from the asymptotic value of the other nearby data points. For example, for pz = 0,

the effective mass at t

∈

[12, 14] is noticeably lower than the asymptotic value approached

20

02468101214t0.50.60.70.80.91.0Meffpz=002468101214t0.50.60.70.80.91.0Meffpz=102468101214t0.70.80.91.01.11.2Meffpz=202468101214t0.80.91.01.11.21.3Meffpz=302468101214t1.01.11.21.31.41.5Meffpz=402468101214t1.21.31.41.51.61.7Meffpz=5at t

∈
range t

[6, 11]. Even though these phenomena appear at different t for different pz, the fit

[2, 11] was used for all pz for consistency. All of the fits have χ2/dof < 1 except

∈

for pz = 1, 2 where they are still close to 1. This suggests that the fits are good quality,

and with the lowest χ2/dof being around 0.4, the data does not appear to be overfitted.

The fitted bands decrease as t increases and become asymptotically horizontal as the excited

state terms decay. Eventually, the fit band approaches a horizontal line as the 2pt correlator

becomes dominated by only the ground state term. The results of the fits are summarized

in Table 3.1.

pz
0
1
2
3
4
5

2

A0|
|
7.705(97)e
2.684(28)e
8.469(120)e
2.444(57)e
6.057(205)e
1.225(56)e

10
−
9
−
9
−
8
−
8
−
7
−

E0
0.681(2)
0.730(1)
0.862(2)
1.044(3)
1.245(5)
1.445(7)

2

A1|
|
2.088(30)e
6.491(59)e
1.839(15)e
4.542(41)e
9.529(95)e
1.728(16)e

9
−
9
−
8
−
8
−
8
−
7
−

E1
1.455(14)
1.448(10)
1.531(12)
1.696(20)
1.894(32)
2.092(44)

χ2/dof
0.808(254)
1.234(322)
1.100(360)
0.692(248)
0.397(356)
0.397(587)

Table 3.1: Fit parameters from fitting the 2pt correlator to Eq. 3.4 and the corresponding
χ2/dof for pz
t
∈
units.

2π/(aL) for the a12m310 ensemble. Fits were performed using
[2, 11]. The energy parameters are actually the unitless aEn, i.e. they are in lattice

[0, 5]

×

∈

3.3 Fitting the 3pt Correlator

After fitting the 2pt correlator, the 3pt correlator fit is performed using parameters from

the 2pt correlator fit. Just as Eq. 3.1 was rewritten as Eq. 3.3, in practice, we rewrite Eq.3.2

as

C 3pt (pf , pi, tsep, tins) =

(cid:88)
m,n |

Am(pf )

An(pi)

||

V
m
|

| ⟨

e−Em(tsep−tins)e−Entins,

(3.5)

n
⟩
|

where normalization or kinematic factors are once again absorbed into the Am and An. It is

important to fit the 2pt correlator first because in practice, we do not fit the energy parame-

ters in Eq. 3.2 but rather insert the values obtained from the 2pt correlator fits. This makes

finding the optimal parameters much easier as the time-dependence in the 3pt correlator

is more complex than that of the 2pt correlator. Thus, it is preferable to avoid having fit

21

parameters that are directly multiplied by time. For this particular work, the operator V

was the gluon energy-momentum tensor (EMT) which I will label as Og, consistent with

Ref. [7], given by

−
In particular, we used Og,tt (the time-diagonal component, alternatively written as Og,44), for

Og,µν = FµρFνρ

gµνFρσFρσ.

(3.6)

1
4

which

Og,44|
0
0
⟩
|
⟨

is proportional to

, where

x
⟨

⟩g

0
⟩
|

is the nucleon ground state and

is the

x
⟨

⟩g

fraction of momentum carried by the gluons. The EMT was projected to zero momentum

transfer so that pi = pf .

Just as with the 2pt correlator, I fitted to the first excited state diagonal term, i.e. up

to m = n = 1 in Eq. 3.7. Thus, the fit function was

C 3pt (tsep, tins) =

+

2

A0|
|
A1||
|

0
⟩
|

Og
0
⟨
|
Og
A0| ⟨
0
|

e−E0tsep +

Og
A1| ⟨
0
|

A0||
1
⟩
|
|
e−E1(tsep−tins)e−E0tins +
A1|
|

1
⟩
|

e−E0(tsep−tins)e−E1tins

2

Og
1
|
⟨

1
⟩
|

e−E1tsep,

(3.7)

where momentum labels have been omitted for brevity. This equation simplifies from the

form of Eq. 3.2 because of the identical initial and final momenta/lack of momentum transfer.

Because the energies from the 2pt correlator fit are inserted into this equation, the only fitted

parameters are the products

||
momenta are the same, the states

Am
|

An

m
|

| ⟨

Og
m
|
and

n
|
n
|

. Additionally, because the initial and final
⟩
are at the same momentum, so when m = n,

⟩
they describe identical states. Thus, the m = n terms do not depend on tins, and the off-

⟩

diagonal matrix elements are related as

Og
1
|
⟨

matrix elements of the EMT are real, we have

=

0
⟩
|
Og
1
|
⟨

Og
0
|
=

⟨
0
⟩
|

∗. Furthermore, because the
1
⟩
|
Og
0
|

. This means that there

1
⟩
|

⟨

are only 3 parameters to fit rather than 4 (8 if the real and imaginary parts of complex

parameters are considered separate) despite the presence of 4 matrix elements in Eq. 3.7.

After fitting, the

are divided out so that we are left with only the matrix element

|
(and kinematic factors that must also be removed).

||

Am
|

An

In order to choose the fit range, we observe that for 0

tins ≪

≪

tsep, Eq. 3.5 is dominated

by the ground state term (m = n = 0). This corresponds to the middle-most tins between

0 (creation) and tsep (annihilation). If suitable tins and tsep are used such that both the 3pt

22

and 2pt correlators are dominated by the ground state, then the special ratio

=

R

C 3pt(pf , pi, tsep, tins)
C 2pt(pi, tsep)

(cid:115)

C 2pt(pi, tins)C 2pt(pi, tsep)C 2pt(pf , tsep −
C 2pt(pf , tins)C 2pt(pf , tsep)C 2pt(pi, tsep −

tins)
tins)

(3.8)

approaches the ground state matrix element. Since pi = pf = p, the square root term

is 1 and Eq. 3.8 becomes a simple ratio of C 3pt(p, tsep, tins)/C 2pt(p, tsep).

In the limit

0

tins ≪

≪

tsep, a plot of Eq. 3.8 as a function of tins approaches a plateau at the ground state

matrix element (for pi

= pf , it can look more like a stationary inflection rather than a local

maximum). To determine a suitable fit range, we plot this ratio as a function of tins for mul-

tiple choices of tsep and examine the plateaus. Because agreement with this limit is strongest

when tins is exactly in the middle of 0 and tsep, we symmetrically exclude endpoints such that
tcut points are excluded from each side, so the fit range becomes tins ∈
each tsep. Much like the effective mass plot analysis, when excited states are desired in the

[tcut, tsep −

tcut] for

fit, we make sure to include some tins outside of the plateau region since the plateau arises

when only the ground state is dominant. When 3pt correlator data with multiple tsep are

fitted as part of a single fit, this is called a simultaneous fit or “sim fit," and the fit function

is essentially a function of multiple variables (tins and tsep). Since we included terms in the

fit function up to the first excited (two states altogether), we refer to these 3pt correlator fits

as “two-sim fits." Using the same lattice configurations that were used for the 2pt correlator

analysis in the previous section, I fitted the 3pt correlator for the a12m310 ensemble using

[7, 11] and tcut = 1. The fitted parameters and fit functions for both the 3pt and 2pt

tsep ∈
correlators are then used to make reconstructed fit bands of Eq. 3.8, and these bands are

plotted with the original data points and a horizontal band to denote the location of the

ground state matrix element

Og,tt
0
|
⟨

0
⟩
|

in Fig. 3.2.

Only the real part of the correlators and ratio produce a signal for the gluon EMT.

Because different values of tsep result in different ranges of tins such that the largest tins is
tsep, the data is plotted vs. tins −
centered around the middle of the plot at tins −
on the same plot, they would not align very well as the data for different tsep would end at

tsep/2 (rather than tins) so that the data for each tsep is

tsep/2 = 0. Otherwise, if all tsep were plotted

23

̸
∈

[0, 5]

Figure 3.2: Plots of the 3pt/2pt correlator ratio defined in Eq. 3.8 vs. tins −
the data for different tsep align and can be compared with each other) for
2π/(aL) for the a12m310 ensemble. Fits were performed with tsep ∈
pz
tcut = 1. Data for different tsep are offset horizontally from each other so that data points
that otherwise overlap can be distinguished. The gray band denotes
pz = 4, 5
increases with tsep and with momentum.

2π/(aL), some of the data goes outside the bounds of the plot because the error

tsep/2 (so that

. For
0
⟩
|

[7, 11] and

Og,tt
0
|
⟨

×

×

24

−4−2024tins−tsep/2−0.2−0.10.00.10.20.30.40.50.6Re[R](pz=0)a12m310tsep=7tsep=8tsep=9tsep=10tsep=11−4−2024tins−tsep/2−0.2−0.10.00.10.20.30.40.50.6Re[R](pz=1)a12m310tsep=7tsep=8tsep=9tsep=10tsep=11−4−2024tins−tsep/2−0.20.00.20.40.6Re[R](pz=2)a12m310tsep=7tsep=8tsep=9tsep=10tsep=11−4−2024tins−tsep/2−0.4−0.20.00.20.40.60.8Re[R](pz=3)a12m310tsep=7tsep=8tsep=9tsep=10tsep=11−4−2024tins−tsep/2−0.50−0.250.000.250.500.751.00Re[R](pz=4)a12m310tsep=7tsep=8tsep=9tsep=10tsep=11−4−2024tins−tsep/2−1.5−1.0−0.50.00.51.01.5Re[R](pz=5)a12m310tsep=7tsep=8tsep=9tsep=10tsep=11different tins. With this choice of the horizontal axis, the ratio data for different tsep should
overlap at tins −
qualitatively if the data behaves as expected. In order to convert

, so it is easier to determine
0
⟩
|

tsep/2 = 0 and approach the value of

Og,tt
0
|
⟨

into

, a

Og,tt
0
|

0
⟩
|

⟨

x
⟨

bare
g
⟩

kinematic factor must be divided out.

In Fig. 3.2, the gray band denoting

Og,tt
0
|

0
⟩
|

⟨

and plateau region (near tins −

tsep/2 = 0)

of the reconstructed fit bands intersect as expected. For pz

2π/(aL), the error on the

3

×

≥

data points becomes noticeably larger such that the data for the largest tsep no longer form

a clear plateau. The same occurs for increasing tsep, although the increase in the amount of

error is greater as pz increases. This causes the data for the largest pz and tsep to no longer

form a clear plateau. As a result, the vertical axis range for pz = 4, 5

2π/(aL) is reduced

×

so that the gray band and data with smaller error can be seen more clearly. This results in

error bars for the largest tsep extending past the plots’ vertical range. All of the fits have

χ2/dof < 1. The fitted parameters and matrix element are displayed in Table 3.2.

2

2

pz
0
1
2
3
4
5

A0|
M00
|
3.488(460)e
10
−
1.320(129)e
9
−
4.146(494)e
9
−
1.369(218)e
8
−
3.767(917)e
8
−
9.538(3569)e
8
−

A1|M01
A0||
|
5.950(933)e
10
−
2.127(246)e
9
−
6.639(781)e
9
−
2.014(275)e
8
−
4.928(963)e
8
−
1.203(307)e
7
−
Table 3.2: Fit parameters for the 3pt correlator two-sim fits to Eq. 3.7 and the
corresponding χ2/dof for pz

A1|
M11
|
1.158(741)e
8
−
2.629(1409)e
8
−
1.255(395)e
7
−
3.196(1627)e
7
−
5.387(7204)e
7
−
3.281(28207)e
7
−

Og,tt
0
0
⟨
⟩
|
|
0.453(59)
0.492(48)
0.490(58)
0.560(89)
0.622(151)
0.778(289)

2π/(aL) for the a12m310 ensemble.

−
−
−
−
−
−

χ2/dof
0.541(507)
0.307(386)
0.454(401)
0.259(252)
0.354(360)
0.221(194)

ij =

i
⟨

Og,tt
|

j
|

[0, 5]
. Fits were performed using tsep ∈

×

∈

M
⟩
pz), the error on some of the parameters is on the same order as the parameters
themselves. For pz = 5, the error on
2
A1|
|
than the parameter.

must be divided by a kinetic factor to convert it to

M11 is actually an order of magnitude larger
bare
x
g
⟩
⟨

Og,tt
0
|

0
⟩
|

⟨

.

[7, 11]. For the noisiest data (largest

25

CHAPTER 4

NONPERTURBATIVE RENORMALIZATION (NPR)

Some quantities calculated via LQCD must be renormalized in order to be compared to

the same quantities calculated by other means, including other lattice calculations that use

different definitions of quantities. For example, the gluon momentum fraction of the nucleon

is independent of the definition of the gauge energy-momentum tensor (EMT) (e.g. clover

vs. plaquette definitions) after nonperturbative renormalization (NPR) is applied to the

EMT [8].

This chapter references a paper of which I am a co-author [7], and the use of first person

singulars like “I" indicates that I am referring to the work that I specifically completed

whereas plurals like “we" refer to the work in that paper completed mostly by the primary

author (though my contributions involved double-checking some of the work). We used lattice

configurations of Nf = 2 + 1 + 1 highly-improved staggered quarks (HISQ) [9] generated by

the MILC collaboration [10]. The a12m310 ensemble analyzed in Chapter 3 comes from

these configurations, and their parameters are summarized in Table 4.1.

ensemble
a (fm)
L3
Lt
×
(MeV)
M val
π
(MeV)
M val
ηs

Ncfg

a09m310
0.0888(8)
323
96
×
313.1(13)
698.0(7)
1009

a12m310
0.1207(11)
243
64
×
309.0(11)
684.1(6)
1013

a15m310
0.1510(20)
163
48
×
319.1(31)
687.3(13)
900

Table 4.1: Lattice spacing a in fm, lattice size L3
×
and M val
in MeV, and number of configurations Ncfg for the lattice ensembles analyzed in
ηs
Chapters 3 and 4. Ncfg applies only to the nucleon correlators analyzed in Chapter 3 and
not to the gluon functions analyzed in Chapter 4. Data taken from Ref. [7].

Lt, valence pion and ηs masses M val
π

Before any calculations are performed, the lattice gauge links undergo 5 steps of hyper-

cubic smearing (5-HYP) [11] because this is the smearing that was applied to the hadronic

matrix element to reduce statistical uncertainties [12]. The gluon fields are then gauge-fixed

26

in the Landau gauge by imposing ∂µAµ(x + ˆµ/2) = 0 with
(cid:21)
(cid:18)

(cid:19)

(cid:20)Uµ(x)

Uµ(x)†

Aµ

x +

ˆµ

=

1
2

−
2ig0a

,

traceless

(4.1)

where the term in brackets has been made traceless (over color). The gauge-fixed gluon fields

are eventually used to calculate Aµ(p) with the Fourier transform
(cid:19)

(cid:18)

Aµ(p) = a4 (cid:88)

eip·(x+ 1

2 ˆµ)Aµ

x +

ˆµ

.

1
2

x∈Λ

(4.2)

4.1 Regularization-independent Momentum Subtraction (RI/MOM)

Scheme

In my particular work, I calculated the necessary renormalization factors for converting

the bare gluon momentum fraction of the nucleon, denoted as

bare
g

x
⟨

⟩

(where x is the momen-

tum fraction), calculated from LQCD to the the modified minimal substraction scheme MS

(MS-bar) using the regularization-independent momentum-subtraction (RI/MOM) scheme

[8]. The gluon momentum fraction is equal to the overlap with the ground state of the

temporal component of the energy-momentum tensor (EMT), Og,tt, after kinematic factors

that arise from the normalization of the fields and the choice of projection are divided out.

That is,

x
⟨

⟩g ∝ ⟨

Og,tt
0
|

, where
0
⟩
|

Og,µν = FµρFνρ

1
4

−

gµνFρσFρσ,

where Fµν is the gluon field strength tensor and gµν is the metric. Thus, renormalizing

is equivalent to renormalizing the operator Og.

(4.3)

x

⟨

⟩g

The renormalization procedure outlined here is derived in Refs. [8, 13] and also utilized in

our work [7]. The renormalized gluon momentum fraction in the MS scheme can be written

as

x
⟨

MS
g = Z MS
Og
⟩

(cid:0)µ2, µ2

R

(cid:1)

g = RMS (cid:0)µ2, µ2
bare
⟩

R

x
⟨

(cid:1) Z RI

Og

(cid:1)

(cid:0)µ2

R

x
⟨

bare
g
⟩

,

(4.4)

where the Z’s are renormalization constants, µ is the matching scale, µR is the renormaliza-

tion scale, and RMS(µ, µR) is the one-loop perturbative matching ratio given by
(cid:19)

(cid:19)

(cid:21)

g2Nf
16π2

(cid:20) 2
3

ln

(cid:18) µ2
µ2
R

−

+

10
9

g2Nc
16π2

−

(cid:18) 4

2ξ +

3 −

ξ2
4

RMS (cid:0)µ2, µ2

(cid:1) = 1

R

.

(4.5)

27

In this equation, derived in Ref. [14], Nf is the number of flavors, Nc is the number of colors,

g2 = 4πα(µ) is the coupling strength (with α(µ) as the coupling constant), and ξ is the

parameter from Riemann’s zeta function. In this work, Nf = 4, Nc = 3, ξ = 0 in the Landau

gauge, and µ = 2 GeV. The RI/MOM renormalization constant Z RI

g (µ2

R) can be obtained

from the normalization condition

Zg

(cid:0)p2(cid:1) Z RI

Og

(cid:0)p2(cid:1) Λbare

Og

(cid:0)p2(cid:1) (cid:16)

Λtree
Og

(cid:0)p2(cid:1)(cid:17)−1(cid:12)

(cid:12)
(cid:12)
(cid:12)p2=µ2

R

= 1,

(4.6)

where Zg (p2) is the gluon-field renormalization constant, and ΛOg (p2)’s are amputated

Green’s functions for the operator Og for gluons in the Landau gauge.

In Refs. [8, 13],

from this condition, it is shown that

(cid:16)

(cid:17)−1 (cid:0)µ2

R

(cid:1) =

Z RI
Og

p2

(Og,µµ
⟨
2(p2

µ −

Og,νν) Tr[Aτ (p)Aτ (
−
−
p2
p)]
Tr[Aτ (p)Aτ (
ν)
⟩
⟨

−

p)]
⟩

(cid:12)
(cid:12)
(cid:12)
(cid:12)p2=µ2

R,τ ̸=µ̸=ν,pτ =0

,

(4.7)

where there is no summation over repeated indices. The quantity Tr[Aµ(p)Aν(

−

gluon propagator or gluon 2pt function and is also written as

Dg,µν(p) = Tr[Aµ(p)Aν(

−

p)],

such that the quantity used in Eq. 4.7 is Dg,τ τ .

4.2 Cluster-decomposition Error Reduction (CDER)

p)] is the

(4.8)

Naïve attempts to calculate Z MS
Og

proved difficult as a signal-to-noise ratio under 100%

was obtained [8, 7]. In order to improve the signal, a technique was developed in Ref. [8]

called cluster-decomposition error reduction (CDER). The motivation for this method is that

the correlators decay exponentially in the distance between operator insertions (e.g. gluon

fields and the EMT), so there is no point in integrating beyond the correlation length as it

will only pick up noise. This is implemented by imposing cutoffs in the spatial integrals used

28

to calculate the correlators. In this specific work, they become [7]

Dg,τ τ =

Tr[Aτ (p)Aτ (
⟨

C CDER
3

(p) =

Og,µν Tr[Aτ (p)Aτ (
⟨

p)]

−

⟩ ≈

p)]

−

⟩ ≈

(cid:28)(cid:90)

|r′|<r2

(cid:28) (cid:90)

(cid:90)

d4r′

(cid:90)

d4r

d4x eip·r′

Tr[Aτ (x)Aτ (x + r′)]

(cid:29)

(cid:90)

d4r′

d4x eip·r′

|r|<r1

|r′|<r2

Og,µν(x + r) Tr[Aτ (x)Aτ (x + r′)]

×

(cid:29)

.

(4.9)

In Ref. [8], it was found that r1 ≈
approximately a factor of 2V / (π2r2

0.9 fm and r2 =
≈
2) where V is the lattice volume.

1r2

1.3 fm reduces the uncertainty by a

Although the CDER method outlined above seems robust and reliable for improving

the signal-to-noise ratio, we took a slightly different approach that was motivated by these

results. Rather than using a hyperspherical cutoff, we used a hypercubic cutoff by truncating

our lattices into 16 L4
c

sub-lattices such that Lc

2r1, 2r2. This method also increases

≈

the number of measurements that we are able to make from Ncfg to approximately 16Ncfg

where Ncfg is the original number of lattice gauge configurations.1 The truncations were

implemented by selecting a random lattice site to serve as the center of the sub-lattice.

Then, we create the sub-lattice by selecting the Lc points in every dimension that surround

this center. For example, if lattice site (n1, n2, n3, n4) is the randomly chosen center, then

the Cartesian product

4
(cid:89)

µ=1

(cid:26)

nµ

−

Lc
2

, nµ

Lc
2

−

+ 1,

· · ·

, nµ, nµ + 1,

,

Lc
2 −

· · ·

(cid:27)
1

is the set of points in the sub-lattice. For lattice with L divisible by 4, the center is not

actually the true center because the Lc sites in a given dimension must be chosen such that

Lc/2 are on one of the center and Lc/2

−

1 are on the other side (because the center is

one of the Lc sites. However, this does not really matter; the lattice site that is randomly

selected could also be chosen to be one of the corners instead. My only recommendation is

that the randomly selected point have a consistent definition on each sub-lattice. On each

1It’s approximate because some data from some of the truncated lattices could not be successfully

transferred from the remote cluster to my local machine.

29

sub-lattice, lattice quantities like Dg,µν(p) and FµνFρσ are calculated as they would be for

the entire lattice for all combinations of spacetime indices and several p such that p2

100

≤

in lattice units.

Our work in Ref. [7] shows that this method is successful in increasing the signal-noise

ratio, and different choices of Lc produce consistent results on a given ensemble. This

was demonstrated by calculating

R) using the full lattice and using different
truncation sizes and plotting the results to show that the truncated lattices produce much

(µ2, µ2

Z MS
Og

(cid:16)

(cid:17)−1

smaller error. For the final NPR calculations, we decided on Lc = L/2 where L is the spatial

extent of the original lattice, which corresponds to aLc

[1.2, 1.4] fm and r1 ≈

∈

0.7 fm,

consistent with the choices of cutoff in Ref. [8]. This information is summarized in Table 4.2.

ensemble
L
Lc
aLc (fm)
Ncfg
Nmeas

a09m310
32
16
1.44
347
5552

a12m310
24
12
1.44
409
6543

a15m310
16
8
1.2
900
14,388

Table 4.2: The original spatial extent L, truncated lattice size Lc = L/2, physical size aLc
of each lattice, number of gauge configurations Ncfg, and number of measurements
Nmeas ≈

16Ncfg.

4.3 Renormalized Gluon Momentum Fraction

(cid:16)

(cid:17)−1

Z RI
Og

R) is calculated according to Eq. 4.7. Different combinations of τ , µ, and ν
satisfy the condition given in the equation, but in this work, we only used combinations

(µ2

such that µ = (τ mod 4) + 1 and ν = (µ mod 4) + 1. Thus, each value of

(cid:17)−1

(cid:16)

Z RI
Og

(µ2

R)

used for future calculations is the average of four different values that correspond to the

four different possible choices of τ (one for each spacetime dimension). Additionally, when

calculating

(Og,µµ
⟨
Tr[Aτ (p)Aτ (

−

Og,νν) Tr[Aτ (p)Aτ (

−

−

, we subtracted the means of (Og,µµ
p)]
⟩

−

Og,νν) and

p)] from themselves before multiplying them together. This is done because

the means of each quantity are supposed to be 0 anyway, so it improves the quality of the

end result.

30

In order to calculate the renormalized gluon momentum fraction

, which is calculated

x
⟨

⟩g

for 0 momentum transfer in the 3pt correlator, the renormalization constant Z MS

RMS (µ2, p2) Z RI

Og (p2) must be extrapolated to p2 = 0. The calculated

divided by RMS (µ2, p2) calculated by Eq. 4.5 to obtain

relatively small dependence of

(cid:16)

Z MS

Og (µ2, p2)

(cid:17)−1

(cid:16)

(cid:17)−1

Z MS
Og
on p2, we fit the data to

(cid:16)

Og (µ2, p2) =
(cid:17)−1
(p2) are

Z RI
Og
(µ2, p2). Because of the

(cid:16)

Z MS
Og

(cid:17)−1 (cid:0)µ2, p2(cid:1) =

(cid:16)

Z MS
Og

(cid:17)−1

+ c1p2 + c2p4,

(4.10)

where

(cid:17)−1

(cid:16)

Z MS
Og

is the reciprocal of the renormalization constant at p2 = 0. This functional

form is chosen because its constituents contain p2, as seen in Eqs. 4.5–4.7. Thus, any Taylor

expansions will only contain even powers of

. The relatively weak dependence of the data
p
|
|
points on p2, as seen in Fig. 4.1, justifies fitting only up to quadratic order in p2. It should

be noted that previous works only fitted to a constant [15] or to linear order in p2 [13].2

The smallest p used for each fit are 2.4, 2.0, and 1.5 GeV for a15m310, a12m310, and

a09m310, respectively, just like PNDME in their study on the same mixed-action [15] as in

this study, in order to avoid non-perturbative effects. As for the upper bound of the fitting

range, it is shown in our previous work [7] that

(cid:17)−1

(cid:16)

Z MS
Og

is mostly consistent for different

choices of the upper bound provided that the upper p2 is sufficiently large. Because of this,

I did not impose an upper bound. The results are shown in Fig. 4.1 and Table 4.3.

ensemble
(cid:17)−1
(cid:16)

Z MS
Og

a09m310

a12m310

a15m310

1.562(69)

1.403(45)

0.998(20)

Table 4.3: MS renormalization constants
by fitting our calculated data points to Eq. 4.10. These values correspond to the
reconstructed fit band at p2 = 0 in Fig. 4.1.

calculated with the RI/MOM scheme

Z MS
Og

(cid:16)

(cid:17)−1

2Other works [8, 13] used the dimensionless (ap)2 as the expansion parameter, but this does not change
the results since the 0th-order term is unaffected and the a2n can be absorbed into the fitting parameters
cn.

31

(cid:16)

(cid:17)−1

Z MS
Og

Figure 4.1:
(µ2, p2) vs. p2 for the a09m310 (upper left), a12m310 (upper right),
and a15m310 (lower) ensembles. The minimum p2 for each fit are 2.42, 22, and 1.52 GeV2
for a09m310, a12m310, and a15m310, respectively, as in Ref. [15]. Thus, the reconstructed
fit bands do not pass through data points for the smallest p2 for a09m310 and a12m310.
The data shows fairly weak dependence on p2, which justifies a fit that is quadratic in p2.

4.4 NPR Calculation Workflow

I will now explain the sequence of programs and calculations required to turn lattice

data to calculated renormalization constants in the RI/MOM scheme. In particular, I will

explain how the relevant equations that are stated earlier in this chapter were practically

implemented to calculate necessary quantities. There are two major sections of the code.

The first is a Chroma3 code that must be run on a supercomputer cluster. The second is a

Mathematica4 notebook that can be run on any modern personal computer with a multi-core

CPU, although around 100 GB of available disk space is recommended to store the necessary

input files. The workflow for each of these sections is outlined in the flowcharts in Figs. 4.2

3https://jeffersonlab.github.io/chroma/
4https://www.wolfram.com/mathematica/

32

a09m310,Lc=16,MS010203040506070012345p2(GeV2)(ZOMS_)-1a12m310,Lc=12,MS0102030405060700.00.51.01.52.02.53.0p2(GeV2)(ZOMS_)-1a15m310,Lc=8,MS010203040500.00.51.01.52.0p2(GeV2)(ZOMS_)-1and 4.3.

Figure 4.2: Outline of the steps for running the first major section of the NPR code, i.e.
the Chroma section. We repeated the truncation steps (last two steps) 16 times for each
lattice gauge configuration.

4.4.1 Chroma

All of the following calculations were done with the Chroma C++ software library running

on the Fermilab National Accelerator Laboratory (FNAL) lattice QCD supercomputer clus-

ter, abbreviated here as FNAL. The process of calculating the NPR constant begins with

33

 Input lattice gauge configurationApply 5-HYP gauge link smearingFix to Landau gaugeFull lattice NPR calculation of Tr𝐴!𝑝𝐴!−𝑝and 𝐹"#𝐹!$for various momenta, for every combination of spacetime indicesPick random lattice site for center of truncation, create sub-lattice from 𝐿%contiguous lattice points in each dimension, including the centerCalculate Tr𝐴!𝑝𝐴!−𝑝and 𝐹"#𝐹!$for various momenta, for every combination of spacetime indices16× Figure 4.3: Outline of the steps for running the second major section of the NPR code, i.e.
the Mathematica section. I organized the notebook such that after the necessary
parameters are entered, the rest of the notebook can be run at once and treated as a black
box. These parts of the code are denoted by the black rectangle.

34

 Make sure Parameters section of Mathematica notebook is accurateUse 𝐹!"𝐹#$to calculate 𝑂%,!"Use 𝑂%,!", Tr𝐴#𝑝𝐴#−𝑝, and 𝑝to calculate 𝑍’!()*+𝑝,, average Tr𝐴#𝑝𝐴#−𝑝with identical 𝑝,Calculate the normalized 2pt correlator, calculate the 3pt correlator, then jackknife bothCalculate 𝑅-.𝜇,,𝑝,, average over identical 𝑝,,combine with 𝑍’!()*+𝑝,to calculate 𝑍’!-.*+𝜇,,𝑝,Fit 𝑍’!-.*+𝜇,,𝑝,vs. 𝑝,to calculate 𝑍’!-.*+at 𝑝,=0, renormalize 𝑥%using a lattice gauge configuration as input. Chroma also requires XML files as input that

specify what is to be calculated. One of the elements in the XML file is the path of the

lattice file. All of the XML files are generated by Perl scripts that require only a couple

of command line arguments to specify either the lattice location and/or the output location

(if the particular step outputs another lattice file). They print to console/stdout, so the

output must be redirected into an XML file. Each time the Chroma executable is run, it

sends output to stdout, but it is not all necessary for future calculations.

The FNAL cluster uses Slurm as its job manager, so our work was completed by submit-

ting sbatch scripts. There are a few variables in the script that need to be set first. They

are listed below:

• INDIR: The submission directory for the script containing the Perl scripts and direc-

tories that will contain the output XML files. Works best if set to SLURM_SUBMIT_DIR.

• OUTDIR: Set equal to INDIR.

• HYP1LOC: Directory containing the input lattice files. These files already have one step

of hypercubic smearing (1-HYP) applied to them, so we need to apply four more.

• HYP5LOC: Directory that will contain the output lattices that have 5-HYP smearing

applied.

• LATFILENAME: The portion of the lattice file names up to and not including the first

dot, e.g. l48144f211b672m0048m024m286a if the lattice file is

l48144f211b672m0048m024m286a.hyp.2100

• Ls: Integer spatial extent of the lattice.

• Lt: Integer temporal extent of the lattice.

• exe: Path to Chroma executable.

• exe1: Path to lattice truncation part 1 executable.

35

• exe2: Path to lattice truncation part 2 executable.

• trunc_out: Directory that will contain the truncated lattice files.

• OMP_NUM_THREADS: Number of threads used for thread parallelism. Should be set to 1.

• GEOM: String of four integers, each separated by a space. The four integers must multiply

to the number of cores/tasks requested. The first three should ideally be identical and

must divide the spatial extent of the lattice. The last integers must divide the temporal

extent. For example, if 512 tasks are requested for a 483

144 lattice, we could have

×

GEOM="4 4 4 8". These numbers must also evenly divide the L4
c

sub-lattices.

In addition to initializing environment variables, certain directories must be created be-

fore the script can be run. They are:

• HYP5LOC: The user must manually create the directory whose path is the value of this

variable in advance.

• trunc_out: The user must manually create the directory whose path is the value of

this variable in advance.

• INDIR/xmls

• INDIR/xmls.GTG

• For each choice of truncated size Lc, there must be a folder INDIR/xmls.LX.GT where

X is replaced by the truncation size with no leading zeroes. For example, Lc = 8 and

Lc = 12 require INDIR/xmls.L8.GT and INDIR/xmls.L12.GT, respectively.

The first step that we perform is applying the 5-HYP smearing to the gauge links. This is

done by running the HYP.pl script where the first argument is the path to the lattice gauge

configuration (i.e. the ones produced by the MILC collaboration) and the second argument

36

is the path where the output lattice file that has been smeared should be created. Run exe

with the XML created by HYP.pl as the input.

The second step is applying the gauge-fixing. This is done by running the GF.pl script

where the first argument is the path to the newly-created smeared lattice from the previous

step and the second argument is the path where the output lattice file that has been gauge-

fixed will be created. Run exe with the XML created by GF.pl.

Third, the full lattice NPR calculation is performed. This is done by running the

NPR.gluon.gf.pl script where the first argument is the path to the newly-created gauge-

fixed lattice file from the previous step and the second argument is the maximum p2 in

integer units (100 in our work). Run exe with the XML created by NPR.gluon.gf.pl. The

details of how this section of the code works and some of the difficulties with this step will

explained later. This step outputs both FµνFρσ and Dg,µν(p) for momenta up to the square

root of the second argument to the NPR.gluon.gf.pl. For FµνFρσ, only the entries with a

Wilson length z of 0 are needed, i.e. the local operator. z specifies the number of sites in

the t-direction between Fµν and Fρσ. The data is organized as such. For FµνFρσ, a line of

data is formatted as:

FF z µ ν

ρ σ Re [FµνFρσ]

Im [FµνFρσ] ,

and a line of Dg,µν data is formatted as:

G2ptp µ ν

px

py

pz

pt Re [Dg,µν]

Im [Dg,µν] .

The string literals FF and G2ptp appear on every line of data to indicate what kind of data

is on that line of output. For a nonzero Wilson length z, FµνFρσ is actually

FµνU4(z)Fρσ,zU4(z)† where the index 4 indicates that the Wilson line U4(z) is in the t-

direction and Fρσ,z indicates that Fρσ has been shifted z lattice units in the t-direction.

Next, the lattice data is restructured in preparation for truncation. This is done by run-

ning the trun_xmlpart1.pl script where the fist argument is the gauge-fixed lattice from be-

fore and the second argument is the path where the output lattice file will be generated. Run

37

exe1 with the XML created by trun_xmlpart1.pl. Then, the lattice is actually truncated.

The sbatch script loops over each desired Lc. Each truncation is performed 16 times at a

random part of the lattice each time. The scripts for this section are trun_xmlpart2_LX.pl

where X is the value of Lc. They require the following six command line arguments in this

order: the path to the restructured lattice file that was just created, the path to the re-

structured lattice file up to and excluding .vec (this is used to create the output lattices

that have been truncated), the x-coordinate of the truncation center, the y-coordinate of the

truncation center, the z-coordinate of the truncation center, the t-coordinate of the trunca-

tion center, each of which is an integer. Each of the 16 randomly-chosen center coordinates

produces a different XML file. Run exe2 once with each XML file.

Finally, the NPR calculation on the truncated lattices is performed. The sbatch script

searches for the truncated lattice files of a given Lc and loops over them. Then the script

NPR.gluon.LX.pl where X is the value of Lc is run for each lattice file where the first argu-

ment is the path to the truncated lattice and the second argument is the maximum p2 (100 in

our work). Then, run exe with these XML files that were created. As with the full lattice cal-

culation, the output sent to stdout has data that is needed for the next steps. Place the lines

beginning with FF with a Wilson length of 0 (z = 0) in files named FF_CFG_TRUNC_LC.txt

and the lines of data that begin with G2ptp in files named G2pt_CFG_TRUNC_LC.txt, where

CFG is the gauge configuration number, TRUNC is the truncation number (between 1 and 16),

and LC is Lc. The rest of the lines should not be added to these files. Also, although this

step outputs Dg,µν(p), recall that only the entries with µ = ν

τ are needed. This will be

≡

taken care of later. These files are required for the Mathematica code that will be run next.

4.4.1.1

Implementation of the NPR Calculation

The NPR calculation itself has a few different steps to obtain FµνFρσ and Dg,µν(p) from

the gauge-fixed or truncated input lattices. Because it is the most complicated of the above

steps, the source code deserves further explanation. The following omits any details of the

setup and initialization required for all Chroma codes and instead focuses on the calculations

38

performed.

The NPR code begins by calculating all of the plaquettes with the gauge links (variable

u) and stores them as the variable plane_plaq. This variable has 3 indices. The first index

determines the directions of the plaquette while the other two indices are ν and µ (Eq. 2.18).

The first index ranges from 0–3 and determines the directions of the plaquette by taking

µ and/or ν to be negative and potentially swapping them. Given the plaquette Uµν, the

indices and corresponding directions are as follows:

• 0: Uµν

• 1: Uν−µ

• 2: U−µ−ν

• 3: U−νµ

These maintain the same orientation (clockwise or counterclockwise) as Uµν. All four of

these are needed to calculate Fµν with the plaquette definition (Eq. 4.12). The plaquettes

are calculated with a function called plaquette. Its arguments, in order, are µ, ν, m, n, Uµ,

Uν, dir+1. dir is determined by the above list. This calculates the m

n plaquette

×

m−1
(cid:89)

Uµ(x + iˆµ)

n−1
(cid:89)

j=0

Uν(x + mˆµ + j ˆν)

i=0
m
(cid:89)

k=1

×

Uµ(x + (m

−

k)ˆµ + nˆν)†

n
(cid:89)

ℓ=1

Uν(x + (n

ℓ)ˆν)†,

−

(4.11)

subject to the specified dir. The individual products must be expanded from left to right

in order for the gauge links to be ordered correctly. For plane_plaq, m = n = 1, as in

Eq. 2.18.

Next, the plaquette function is used to calculate 1

1, 1

2, 1

3, 1

4, 2

2, and 3

×
plaquettes in the standard directions (dir=0), stored in variables named plane_plaq_mn,

×

×

×

×

×

2

where m and n are replaced with their numerical values. Then each plaquette is averaged

over color by taking the trace of the real part and dividing by Nc = 3 and averaged over

39

lattice sites, and the results are stored in tr_plane_plaq_mn where m and n are replaced

with their numerical values. These traced, averaged plaquettes are written to the output

XML file. However, they are not used in the NPR calculation.

Next, Fµν(x) is calculated according to the plaquette definition

Fµν(x) =

i
8a2g

(cid:0)U[µν](x) + U[ν−µ](x) + U[−µ−ν](x) + U[−νµ](x)(cid:1) ,

(4.12)

where

U[µν] = Uµν

Uνµ = Uµν

−

U †

µν.

−

(4.13)

Fµν is calculated except for the multiplication by i/(a2g) and stored in F. This is passed

to a function that calculates FµνU4(z)Fρσ,zU4(z)† up to a specified maximum Wilson length

and outputs the results to stdout as described earlier. This function also calculates several

operators that are not needed for the NPR calculation.

is removed by subtracting Tr(cid:2)Uµ(x)

Then Aµ(x) is calculated according to Eq. 4.1 and stored in the variable Ax. The trace
Uµ(x)†(cid:3)/Nc from each diagonal matrix element of
Uµ(x)†. Afterward, Aµ(x) is Fourier transformed to Aµ(p) by Eq. 4.2. To assist

Uµ(x)

−

−

with this, the SftMom class provided by Chroma stores the phases and available momentum

4-vectors that appear in the complex exponentials eip·x in the Fourier transform in the

variable phases. Another variable shift_phase stores the additional shift eip·ˆµ/2. The

Fourier transformed gluon field is stored in Ap. This is finally used to calculate the gluon 2pt

function Dg,µν(p) by calculating Tr[Aµ(p)Aν(

p)]. The results are then printed to stdout

−
as described earlier. Note that this section of the code calculates a few other quantities that

contain Dg,µν(p) and prints them, but they are not needed for any future calculations. It

appears that some calculations that are performed in the Mathematica notebook may have

been attempted in this code, but they were likely scrapped in favor of the Mathematica

notebook.

40

4.4.1.2 Outstanding Issues

Much like the Mathematica notebook that will be discussed next, the original author

of this NPR calculation code left room for improvements to the efficiency and readability.

This led to some issues when I tried to run the NPR calculation on the a06m310 ensemble,

which is a 483

×

144 lattice with lattice spacing a

≈

0.06 fm. This is noticeably larger (in

terms of total number of sites, and this is needed to avoid finite volume effects) than the

other ensembles used in our work, so this lattice has a larger memory requirement than our

previous runs. This caused us to run out of memory on the FNAL cluster even when running

the full lattice NPR calculation on 16 nodes. We are unable to request more nodes than

this, so I attempted to improve the efficiency of the code.

One change that I implemented was to ensure that vector-like objects that are passed as

parameters to functions are passed by reference, but this did not resolve the issue.5 Despite

this, I left this modification in the code because this is proper practice in C++ anyway. I

tried to profile the heap memory usage with the massif tool in valgrind, but this failed

to provide meaningful information. It is likely that the section of the code that causes the

program to run out of memory executes quickly enough that massif does not have time to

log the increased memory usage before the entire program crashes.

All of this occurred on the full lattice calculation, so I tried to run the calculation on the

truncated lattices because they would be smaller in memory. This also did not work because

part of the code that performs the truncation also had issues with memory usage. At the

time of writing, I am at an impasse with the a06m310 ensemble.

4.4.2 Mathematica

Once the Chroma step is completed, its output is used in the Mathematica notebook. The

most important requirement is to set everything correctly in the “Parameters" section at the

very beginning of the notebook. This includes lattice parameters, input directories that

contain the required FµνFρσ and Dg,µν(p), output directories, the matching scale µ (used to

5It’s possible that this is necessary at some point, but if it is, the program crashes before then.

41

calculate the previously-mentioned perturbative matching ratio), and a few other parameters

that are explained in the notebook. The user must also create a text file containing a list of

the lattice gauge configuration numbers that were used and specify which truncations in a

given configuration are not usable, if any (this feature was implemented after I found files

that were empty or incomplete). This is implemented via a map that maps configuration

numbers to a list of truncation numbers that are unusable. Other lines do not need to be

modified by the user, such as certain variables that check whether the files created by the

notebook already exist. The subsequent sections that perform calculations will not run if

their output already exists. This is useful if, for example, one only has already run the

notebook but wants to use the output to make plots at the end. As long as this section is

modified correctly, the rest of the notebook can be run at once and treated as a black box

until the plots at the end of the notebook are created. Also, this section contains a variable

with many values of µ in GeV and the corresponding value of the running coupling constant

α(µ) hard-coded, and an interpolation is performed so that this quantity can be used to

calculate g2(µ) = 4πα(µ). This is needed to calculate RMS (µ2, µ2

R) (Eq. 4.5) later on.

The next section of the notebook calculates the gluon EMT Og,µµ (only the diagonals).

The Og,µµ

−

Og,νν in Eq. 4.7 term means that the FρσFρσ in Eq. 4.3 term does not need to

be calculated because it cancels when diagonal terms of Og are subtracted as in Eq. 4.7.

Thus, the notebook parses through the files containing FµνFρσ to find entries where µ = ρ

and ν = σ, and then it calculates Og. Only the real part of FµνFρσ is needed. Also, the

calculated Og needs to be divided by 2 because of how FρσFρσ was calculated in Chroma.

The next section of the notebook calculates the gluon propagator or 2pt function

Dg,τ τ (p). This section requires the most runtime by far because the files containing the

Dg,µν(p) data have over 300,000 lines. Thus, the notebook is bottlenecked by how quickly

these files can be read rather than by the calculations. However, only a small fraction of

these lines are needed for calculations, and the line numbers of the needed lines are the

same for every file. To improve efficiency, the notebook first parses through the first file and

42

makes a list of the line numbers that have µ = ν

τ and pτ = 0, as required by Eq. 4.7, but

≡

also p[(µ+1) mod 4]+1 = 0. Since we used µ = (τ mod 4) + 1 and ν = (µ mod 4) + 1, where

the indices correspond to Eq. 4.7, this corresponds to pν = 0. This is consistent with the

derivation in Ref. [8]. For example, if τ = t (index 4), then py = 0 is required (in addition to

pt = 0). Then, each file containing Dg,µν(p) is loaded and the gluon 2pt function is calculated

using just the relevant line numbers in the list. In addition to the 2pt function, a quantity

called the normalized 2pt function is also calculated, equal to

2p2

µDg,τ τ (p)
p2

,

−

(4.14)

where the additional terms correspond to the momenta in Eq. 4.7. The negative sign is re-

quired to make the renormalization constant positive and probably comes from how Dg,µν(p)

was calculated in Chroma.

In this process, the momentum 4-vector is mapped separately to the raw and normalized

2pt functions. Then, another map maps p2 in integer units to the 2pt function by averaging

over all values with the same p2. The mapping is done in these two separate steps because

there are multiple definitions used for the momentum. In particular, Refs. [8, 13] reference

a definition motivated by the poles of the lattice propagator, given by

a˜pµ = 2 sin

(cid:17)

(cid:16)apµ
2

= 2 sin

(cid:18) πnµ
Lµ

(cid:19)

,

(4.15)

where the second equality uses pµ = 2πnµ/(aLµ) where nµ

0, 1,

, Lµ

1
}

−

· · ·

∈ {

. Both this

definition and the standard Euclidean definition of the square of the momentum magnitude

are mapped to the 2pt function although analysis primary focuses on the latter in our work.

For this former, the map maps (aˆp)2 in exact symbolic form to the 2pt function. All of

these maps are saved as well as a map of the integer p2 to the p2 in physical units of GeV2

(Euclidean definition).

After the gluon 2pt propagator has been calculated and averaged over identical momenta,

the next section of the notebook calculates the quantity (Og,µµ

Og,νν) Dg,τ τ (p). This quan-

−

tity is called the gluon 3pt function. First, the Og,µµ calculated from the previous section

43

are used to calculate the required differences Ogµµ

Ogνν. Since each τ only uses µ = (τ

−

mod 4) + 1 and ν = (µ mod 4) + 1, for each µ, only one value of ν is required to perform

calculations. Thus, the data is stored in a Table as

Og,yy
{

−

Og,zz, Og,zz

Og,tt, Og,tt

Og,xx, Og,x

Og,yy

.

}

−

−

−

This means that the index of this list is τ , e.g. the first entry corresponds to τ = x (index

1) such that it contains Og,yy

−

Og,zz. The operator is combined with the 2pt function as

(cid:68)(cid:104)

(Og,µµ

Og,νν)

(Og,µµ

Og,νν)

−

−

−

(cid:105) (cid:104)

Dg,τ τ (p)

−

(cid:105)(cid:69)

Dg,τ τ (p)

,

where the overline denotes the mean over all configurations. This calculation subtracts the

mean of each quantity from itself as described earlier. Note that the non-normalized 2pt

function is used in this calculation. Once the 3pt function has been calculated, this quantity

is divided by the normalized 2pt function. The results are stored in maps that map p2 to

the ratios, one map for each momentum definition.

The last section of the notebook begins by averaging the ratios from the previous section

At this point,

over τ as the previous calculation produces a different value for each spacetime coordinate.
(cid:17)−1 (cid:0)µ2 = 4 GeV2, p2(cid:1) is finally
calculated. First, RMS (cid:0)µ2 = 4 GeV2, p2(cid:1) is calculated using the interpolated α(µ) in the

(p2) has been calculated. Then,

Z MS
Og

Z RI
Og

(cid:17)−1

(cid:16)

(cid:16)

“Parameters" section at the beginning of the notebook. This is combined with

(cid:17)−1

(cid:16)

Z RI
Og

(p2)

to calculate

(cid:16)

Z MS
Og

(cid:17)−1 (cid:0)µ2 = 4 GeV2, p2(cid:1) as

(cid:16)

Z MS
Og

(cid:17)−1 (cid:0)µ2 = 4 GeV2, p2(cid:1) =

(cid:16)

Z RI
Og

(cid:17)−1 (cid:0)p2(cid:1) /RMS (cid:0)µ2 = 4 GeV2, p2(cid:1) .

(4.16)

Then, the data points are fitted to Eq. 4.10 using data points that satisfy p2 > p2fitmin

defined in the “Parameters" section. This section also makes a plot of the data points and

the reconstructed fit band on separate plots. To combine them, use the Show function. This

can be done by copying each image and pasting them as the arguments of Show or by saving

each plot as a variable and using those variables as the arguments. Separate plots are also

produced for the alternate definition of momentum in Eq. 4.15.

44

4.4.2.1 My Contribution to the Notebook

The plots in Fig. 4.1 were produced using Mathematica notebooks that were written as

described in this subsection. I did not originally write them, but I made significant modifica-

tions and improvmenets to them. My work involved improving the readability and efficiency

of the code, adding additional comments, double-checking calculations and formulas, and

reducing the required runtime. I improved efficiency and runtime in various ways. I utilized

more efficient data structures for the required information that we needed to store (namely

maps, implemented as the Association data type in Mathematica). This was sensible com-

pared to the list, which functions as a map but with keys restricted to sequential integers,

because the p2 in integer units are not necessarily sequential and not necessarily in ascending

order. Thus, the map also kept the data ordered in terms of ascending momentum. I also

implemented the methods for parsing text files by using the fact that the needed line num-

bers are the same for each type of object (FµνFρσ or Dg,µν(p)). This improved the runtime

of the notebooks from the order of days to the order of hours. In addition, the modularity

of the code allows for one to make additions or modifications, particularly in the last section

where the plots are produced.

45

CHAPTER 5

CONCLUSION

Lattice QCD provides a method for theoretical calculations of phenomenological parameters

and observables. I outlined the basics of LQCD to demonstrate that it reduces to continuum

QCD in the limit of zero spacing, and thus, it is a viable method for calculations of physical

and phenomenological quantities. Then I defined the 2pt and 3pt correlators and fitted them

in order to calculate the matrix element

Og,tt
0
|
⟨

, which is used to calculate
0
⟩
|

x
⟨

bare
g
⟩

. The

effective mass and correlator ratio plots combined with the reasonable χ2/dof indicate that

the fits are of good quality. Afterward, I calculated the necessary renormalization constants

using the RI/MOM method to convert the results to the MS scheme. Using truncated lattices

inspired by the CDER technique, I was able to obtain satisfactory data with a suitable signal-

noise ratio such that the calculated quantities have reasonable error. Lastly, I documented

the many steps involved in calculating the NPR constants and improved the efficiency of

the Mathematica notebook that was previously written for this purpose. This allows for our

results to be replicated and for our method to be applied to other future work, such as the

gluon momentum fraction of mesons.

46

REFERENCES

[1] Karen McNulty Walsh. Physicists zoom in on gluons’ contribution to proton spin, Feb

2016.

[2] Zhouyou Fan, Rui Zhang, and Huey-Wen Lin. Nucleon gluon distribution function from
2 + 1 + 1-flavor lattice qcd. International Journal of Modern Physics A, 36(13):2150080,
2021.

[3] Zhouyou Fan and Huey-Wen Lin. Gluon parton distribution of the pion from lattice

qcd. Physics Letters B, 823:136778, 2021.

[4] Tanjib Khan, Raza Sabbir Sufian, Joseph Karpie, Christopher J. Monahan, Colin
Egerer, Bálint Joó, Wayne Morris, Kostas Orginos, Anatoly Radyushkin, David G.
Richards, Eloy Romero, and Savvas Zafeiropoulos. Unpolarized gluon distribution in
the nucleon from lattice quantum chromodynamics. Phys. Rev. D, 104:094516, Nov
2021.

[5] Christof Gattringer and Christian B. Lang. Quantum Chromodynamics on the Lattice:

An Introductory Presentation. Lecture Notes in Physics. Springer, 2010.

[6] B.A. Berg. Markov Chain Monte Carlo Simulations and Their Statistical Analysis:

With Web-based Fortran Code. World Scientific, 2004.

[7] Zhouyou Fan, Huey-Wen Lin, and Matthew Zeilbeck. Nonperturbatively renormalized
nucleon gluon momentum fraction in the continuum limit of Nf = 2 + 1 + 1 lattice qcd.
Phys. Rev. D, 107:034505, Feb 2023.

[8] Yi-Bo Yang, Ming Gong, Jian Liang, Huey-Wen Lin, Keh-Fei Liu, Dimitra Pefkou,
and Phiala Shanahan. Nonperturbatively renormalized glue momentum fraction at the
physical pion mass from lattice qcd. Phys. Rev. D, 98:074506, Oct 2018.

[9] E. Follana, Q. Mason, C. Davies, K. Hornbostel, G. P. Lepage, J. Shigemitsu, H. Trot-
tier, and K. Wong. Highly improved staggered quarks on the lattice with applications
to charm physics. Phys. Rev. D, 75:054502, Mar 2007.

[10] A. Bazavov, C. Bernard, J. Komijani, C. DeTar, L. Levkova, W. Freeman, Steven
Gottlieb, Ran Zhou, U. M. Heller, J. E. Hetrick, J. Laiho, J. Osborn, R. L. Sugar,
D. Toussaint, and R. S. Van de Water. Lattice qcd ensembles with four flavors of highly
improved staggered quarks. Phys. Rev. D, 87:054505, Mar 2013.

[11] Anna Hasenfratz and Francesco Knechtli. Flavor symmetry and the static potential

with hypercubic blocking. Phys. Rev. D, 64:034504, Jul 2001.

[12] Zhou-You Fan, Yi-Bo Yang, Adam Anthony, Huey-Wen Lin, and Keh-Fei Liu. Gluon
quasi-parton-distribution functions from lattice qcd. Phys. Rev. Lett., 121:242001, Dec
2018.

[13] P. E. Shanahan and W. Detmold. Gluon gravitational form factors of the nucleon and

the pion from lattice qcd. Phys. Rev. D, 99:014511, Jan 2019.

47

[14] Yi-Bo Yang, Michael Glatzmaier, Keh-Fei Liu, and Yong Zhao. The 1-loop correction
of the qcd energy momentum tensor with the overlap fermion and hyp smeared iwasaki
gluon, 2016.

[15] Santanu Mondal, Rajan Gupta, Sungwoo Park, Boram Yoon, Tanmoy Bhattacharya,
and Huey-Wen Lin. Moments of nucleon isovector structure functions in 2 + 1 + 1-flavor
qcd. Phys. Rev. D, 102:054512, Sep 2020.

48