text:
# Next-to-leading order numerical calculations in Coulomb gauge

## I Introduction

QCD calculations at next-to-leading order can be done in a style in which all of the integrations over loop three-momenta are done numerically. In this method, in particular, the integrations for virtual loop diagrams are performed numerically instead of analytically. The method works well in Feynman gauge beowulfPRL ; beowulfPRD ; beowulfrho  when applied to three-jet-like observables in electron-positron annihilation (eg. the three jet cross section).

Feynman gauge is the simplest gauge from a calculational point of view. However, it is unphysical in that unphysical gluon polarizations propagate into the final state. The contributions from unphysical polarizations cancel when one sums over graphs. However, we have in mind applications in which one wants to go beyond a pure next-to-leading order calculation by incorporating, in an approximate way, some effects at all orders in $`\alpha _s`$. (We have in mind, for instance, renormalon and parton showering effects.) For such applications, one must approximate, and the presence of unphysical degrees of freedom propagating over long distances makes it difficult to see what approximations to apply. The remedy is simple: do the calculation in a physical gauge, such as Coulomb gauge.

In this paper, we develop the apparatus needed for applying the numerical integration method in Coulomb gauge.111 We choose Coulomb gauge over space-like axial gauge because we have chosen to treat cross sections in electron-positron annihilation, which has a natural symmetry under rotations in the electron-positron c.m. frame. The choice of Coulomb gauge maintains this symmetry. Time-like axial gauge might have been a good choice, but this choice complicates the structure of amplitudes as a function of the energy in a virtual loop. For the most part, this is straightforward: one should simply replace the Feynman gauge Feynman rules by the Coulomb gauge Feynman rules. However, two point functions and three point functions need a special treatment (in any gauge). Three point one loop virtual subgraphs need a special treatment because they are ultraviolet divergent. The modified minimal subtraction (   MS ) prescription to calculate in $`32ϵ`$ space dimensions and remove poles is not useful for numerical integrations. Thus one must convert the   MS  subtraction to an equivalent subtraction defined in exactly 3 dimensions. (Four point subgraphs with all gluon legs need renormalization too, but such subgraphs do not occur in the applications that we have in mind, so we omit consideration of them.) The two point one loop virtual subgraphs need a special treatment because they need renormalization. They also need a special treatment for another reason. Let $`\mathrm{\Sigma }(q)`$ be the one loop quark self-energy function. Then when the self-energy attaches to a quark line that enters the final state, we need $`\text{/}q\mathrm{\Sigma }(q)\text{/}q/q^2`$ evaluated at $`q^2=0`$. We need to express $`\mathrm{\Sigma }(q)`$ as a numerical integral, but we need to do it in such a way that the integrand for $`\text{/}q\mathrm{\Sigma }(q)\text{/}q/q^2`$ is finite at $`q^2=0`$. The integral will then have a logarithmic infrared divergence, but this infrared divergence will be cancelled by a corresponding divergence in the graph with a cut self-energy diagram. The same consideration applies to the gluon propagator.

This paper will also serve to document the methods needed to treat two-point subgraphs and three-point virtual subgraphs even in Feynman gauge. These methods were discussed briefly in beowulfPRL  but the details were left to unpublished notes beowulfnotes  that accompany the associated computer code beowulfcode .

We have implemented the methods described in this paper in a computer program beowulfcode . The code is based on the Feynman gauge code described in beowulfPRL ; beowulfPRD ; beowulfrho . In its default mode, the program acts as a next-to-leading order Monte Carlo event generator, generating three and four parton final states with corresponding weights. Suppose, for example, one wishes to calculate the expectation value of $`(1t)^n`$, where $`t`$ is the thrust of each event and $`n`$ is fixed. To do this, a separate subroutine calculates $`(1t)^n`$ for each event, multiplies by the corresponding weight, and averages over events. As in other programs of this type, the weights can be positive or negative. The user can specify either Feynman or Coulomb gauge for the calculation. Although the graph-by-graph contributions are very much gauge dependent, the net results are the same for the two gauges.

The plan of this paper is as follows. We begin in Sec. II with a brief review of the main structure of a next-to-leading order calculation by numerical integration. Then in Sec. III we explain the momentum space coordinates used in the analysis of two point functions. In Secs. IV, V and VI, we analyze the one loop gluon self-energy in Coulomb gauge. In Secs. VII, VIII, and IX we turn to the one loop quark self-energy. In Sec. X, we examine the renormalization of three point functions in Coulomb gauge. Finally, in Sec. XI, we present some results from the computer program that implements the formulas from this paper. Appendix A contains formulas for Feynman gauge that correspond to the Coulomb gauge formulas of the main body of the text.

## II Calculations by numerical integration

The idea of this paper is to make the numerical method for next-to-leading order QCD calculations work in Coulomb gauge. Before beginning this task, we need to outline the numerical method itself, which works in any gauge. We present a sketch only since the details can be found in Refs. beowulfPRL ; beowulfPRD ; beowulfrho . Besides the matter of the gauge choice, there is one difference between the computer algorithms presented here and those of Refs. beowulfPRL ; beowulfPRD ; beowulfrho : here we wish to calculate the complete observable at next-to-leading order, that is the sum of the contributions proportional to $`(\alpha _s/\pi )^1`$ and $`(\alpha _s/\pi )^2`$, rather than just the coefficient of $`(\alpha _s/\pi )^2`$. This is a rather trivial change, but it affects the way the problem is set up below.

We wish to calculate an observable $``$ with the following structure

$$=\frac{1}{\sigma _0\left(\sqrt{S}\right)}\underset{n}{}\frac{1}{n!}𝑑\stackrel{}{P}_1\mathrm{}𝑑\stackrel{}{P}_n\frac{d\sigma }{d\stackrel{}{P}_1\mathrm{}d\stackrel{}{P}_n}𝒮_n(\stackrel{}{P}_1,\mathrm{},\stackrel{}{P}_n).$$
(1)
We work in the $`e^+e^{}`$ c.m. frame and $`\sqrt{S}`$ is the c.m. energy.222We also average over the direction of the beam axis relative to the $`z`$ axis of our coordinate system. Thus we can calculate typical event shape variables like the cross section to make three jets, but not correlations between a jet direction and the beam axis. However, nothing in the general methods used here prevents one from removing this simplification. The observable $``$ has been normalized by dividing by the Born level cross section $`\sigma _0`$ for $`e^++e^{}\mathrm{ℎ𝑎𝑑𝑟𝑜𝑛𝑠}`$. The quantity $`d\sigma /[d\stackrel{}{P}_1\mathrm{}d\stackrel{}{P}_n]`$ is the cross section to make $`n`$ massless partons. In defining this cross section, we treat all $`n`$ partons as identical and thus divide by $`n!`$. The functions $`𝒮_n`$ are the measurement functions that define the observable KS . They are symmetric under interchange of any of their variables and have the property of infrared safety:

$$𝒮_{n+1}(\stackrel{}{P}_1,\mathrm{},\lambda \stackrel{}{P}_n,(1\lambda )\stackrel{}{P}_n)=𝒮_n(\stackrel{}{P}_1,\mathrm{},\stackrel{}{P}_n)$$
(2)
for $`0\lambda <1`$. Typically the functions $`𝒮_n`$ are dimensionless, but we do not assume that here. We are concerned with three-jet-like quantities, which means that $`𝒮_2=0`$, so that the smallest value of $`n`$ that contributes is $`n=3`$. We work in next-to-leading order perturbation theory, so that there are at most four partons in the final state. Thus the sum over $`n`$ runs over $`n=3`$ and $`n=4`$.

The parton level cross sections contain delta functions, which we make explicit by writing

$$\frac{d\sigma }{d\stackrel{}{P}_1\mathrm{}d\stackrel{}{P}_n}=\delta \left(\stackrel{}{P}_i\right)\delta \left(|\stackrel{}{P}_i|\sqrt{S}\right)F_n(\stackrel{}{P}_1,\mathrm{},\stackrel{}{P}_n,\sqrt{S},\alpha _s\left(C\sqrt{S}\right)).$$
(3)
The function $`F_n`$ depends on the momenta and on $`\alpha _s`$, which we evaluate at a scale $`\mu =C\sqrt{S}`$, where $`C`$ is a dimensionless parameter of order 1. The order $`\alpha _s^2`$ contributions contain logarithms of $`\mu `$ and thus of $`\sqrt{S}`$, so we have indicated a separate dependence on $`\sqrt{S}`$. The dependence of $`F`$ on the dimensionless parameter $`C`$ is left implicit. Thus we write $``$ as

$``$ $`=`$ $`{\displaystyle \frac{1}{\sigma _0\left(\sqrt{S}\right)}}{\displaystyle \underset{n}{}}{\displaystyle \frac{1}{n!}}{\displaystyle 𝑑\stackrel{}{P}_1\mathrm{}𝑑\stackrel{}{P}_n\delta \left(\stackrel{}{P}_i\right)\delta \left(|\stackrel{}{P}_i|\sqrt{S}\right)}`$ (4)
$`\times F_n(\stackrel{}{P}_1,\mathrm{},\stackrel{}{P}_n,\sqrt{S},\alpha _s\left(C\sqrt{S}\right))𝒮_n(\stackrel{}{P}_1,\mathrm{},\stackrel{}{P}_n).`$
The energy conserving delta function would create problems in a numerical integration, so we get rid of it by the following strategy. We introduce a factor 1 written as

$$1=\sqrt{S}_0^{\mathrm{}}𝑑th\left(t\sqrt{S}\right),$$
(5)
where $`t`$ has the dimensions of time and $`h`$ is any convenient smooth function whose integral is 1. We change integration variables in the momentum integrals to dimensionless variables

$$\stackrel{}{p}_i=t\stackrel{}{P}_i.$$
(6)
Dimensional analysis gives

$`F_n(\stackrel{}{p}_1/t,\mathrm{},\stackrel{}{p}_n/t,\sqrt{S},\alpha _s(C\sqrt{S}))`$ $`=`$ $`t^{3n2}F(\stackrel{}{p}_1,\mathrm{},\stackrel{}{p}_n,t\sqrt{S},\alpha _s\left(C\sqrt{S}\right)),`$
$`\sigma _0\left(\sqrt{S}\right)`$ $`=`$ $`t^2\sigma _0\left(t\sqrt{S}\right).`$ (7)
Thus

$``$ $`=`$ $`{\displaystyle \frac{1}{\sigma _0\left(t\sqrt{S}\right)}}{\displaystyle \underset{n}{}}{\displaystyle \frac{1}{n!}}{\displaystyle 𝑑\stackrel{}{p}_1\mathrm{}𝑑\stackrel{}{p}_n𝑑th\left(t\sqrt{S}\right)\delta \left(\stackrel{}{p}_i\right)\delta \left(|\stackrel{}{p}_i|/\sqrt{S}t\right)}`$ (8)
$`\times F_n(\stackrel{}{p}_1,\mathrm{},\stackrel{}{p}_n,t\sqrt{S},\alpha _s\left(C\sqrt{S}\right))𝒮_n(\stackrel{}{p}_1/t,\mathrm{},\stackrel{}{p}_n/t).`$
Now we can use the integration over $`t`$ to eliminate the energy-conserving delta function. Denoting

$$\sqrt{s}|\stackrel{}{p}_i|$$
(9)
we have

$``$ $`=`$ $`{\displaystyle \frac{1}{\sigma _0\left(\sqrt{s}\right)}}{\displaystyle \underset{n}{}}{\displaystyle \frac{1}{n!}}{\displaystyle 𝑑\stackrel{}{p}_1\mathrm{}𝑑\stackrel{}{p}_n\delta \left(\stackrel{}{p}_i\right)h\left(\sqrt{s}\right)}`$ (10)
$`\times F_n(\stackrel{}{p}_1,\mathrm{},\stackrel{}{p}_n,\sqrt{s},\alpha _s\left(C\sqrt{S}\right))𝒮_n(\stackrel{}{p}_1\sqrt{S/s},\mathrm{},\stackrel{}{p}_n\sqrt{S/s}).`$
Eq. (10) is implemented as an event generator. Events with $`n`$ partons with scaled momenta $`\stackrel{}{p}_i`$ are generated along with a weight equal to $`hF_n/\sigma _0`$ in Eq. (10) divided by the density of points in $`\stackrel{}{p}_i`$ space. A separate routine then multiplies the weights by the measurement function $`𝒮_n(\stackrel{}{p}_1[S/s]^{1/2},\mathrm{},\stackrel{}{p}_n[S/s]^{1/2})`$ and takes the average of these results over a large number of generated points.

Two features are especially worth noting in the main part of Eq. (10), that is in every part other than the measurement functions. First, the true momentum variables have been replaced by dimensionless variables $`\stackrel{}{p}_i`$. Second, the energy conserving delta function has been replaced by $`h(\sqrt{s})`$, so that the dimensionless energy $`\sqrt{s}`$ is not fixed. The true c.m. energy $`\sqrt{S}`$ appears in only one place, in the argument of $`\alpha _s`$. In this calculation, we are using the massless theory. However, if one wanted to add quark masses $`m_i`$, then the functions $`F_n`$ in Eq. (10) would depend on additional dimensionless parameters $`[s/S]^{1/2}m_i`$. Then the c.m. energy $`\sqrt{S}`$ would appear in the argument of $`\alpha _s`$ and in these dimensionless mass parameters.

The contribution $`^{(2)}`$ to $``$ proportional to $`\alpha _s^2`$ can be expressed in terms of cut Feynman diagrams, as in Fig. 1. (In this section, we consider diagrams that do not have self-energy subdiagrams, since self-energy diagrams require a special treatment.) The dots where the parton lines cross the cut represent the function $`𝒮_n`$. Each diagram is a three loop diagram, so we have integrations over loop momenta $`\stackrel{}{l}_1`$, $`\stackrel{}{l}_2`$ and $`\stackrel{}{l}_3`$. Eq. (10) lacks an energy conserving delta function, so we have integrations over four energies, which we might take to be loop energies $`l_1^0`$, $`l_2^0`$ and $`l_3^0`$ chosen in the same way as the loop momenta and the energy $`l_0^0`$ entering the graph on the vector boson line. We first perform the energy integrations. For the graphs in which four parton lines cross the cut, there are four mass-shell delta functions $`\delta (p_J^2)`$. These delta functions eliminate the three energy integrals over $`l_1^0`$, $`l_2^0`$, and $`l_3^0`$ as well as the integral over $`l_0^0`$. For the graphs in which three parton lines cross the cut, we can eliminate the integration over $`l_0^0`$ and two of the $`l_J^0`$ integrals. One integral over the energy $`E`$ in the virtual loop remains. To perform this integration, we close the integration contour in the lower half $`E`$ plane. This gives a sum of terms obtained from the original integrand by some simple algebraic substitutions in which $`E`$ is replaced by a location $`E_i`$ of one of the poles in the lower half $`E`$ plane. Then we do the same thing except that we close the integration contour in the upper half $`E`$ plane. Finally, we take the average of these results. (For well behaved integrands, these two contributions are the same, but in Coulomb gauge some of the integrands are not so well behaved, as we shall see.)

Having performed the energy integrations, we are left with an integral of the form

$$^{(2)}=\underset{G}{}𝑑\stackrel{}{l}_1𝑑\stackrel{}{l}_2𝑑\stackrel{}{l}_3\underset{C}{}g(G,C;\stackrel{}{l}_1,\stackrel{}{l}_2,\stackrel{}{l}_3).$$
(11)
Here there is a sum over graphs $`G`$ (of which one is shown in Fig. 1) and there is a sum over the possible cuts $`C`$ of a given graph. The problem of calculating $`^{(2)}`$ is now set up in a convenient form for calculation.

If we were using the Ellis-Ross-Terrano method for doing next-to-leading order calculations ERT , we would put the sum over cuts outside of the integrals in Eq. (11). For those cuts $`C`$ that have three partons in the final state, there is a virtual loop. We can arrange that one of the loop momenta, say $`\stackrel{}{l}_1`$, goes around this virtual loop. The essence of the Ellis-Ross-Terrano method is to perform the integration over the virtual loop momentum analytically, while the remaining integrations are performed numerically. The integration over the virtual loop momentum is often ultraviolet divergent, but the ultraviolet divergence is easily removed by a renormalization subtraction. The integration is also typically infrared divergent. This divergence is regulated by working in $`32ϵ`$ space dimensions and then taking $`ϵ0`$ while dropping the $`1/ϵ^n`$ contributions (after proving that they cancel against other contributions). After the $`\stackrel{}{l}_1`$ integration has been performed analytically, the integrations over $`\stackrel{}{l}_2`$ and $`\stackrel{}{l}_3`$ can be performed numerically. For the cuts $`C`$ that have four partons in the final state, there are also $`1/ϵ^n`$ infrared divergences. One uses either a “phase space slicing” or a “subtraction” procedure to isolate these divergences, which cancel the $`1/ϵ^n`$ pieces from the virtual graphs. In the end, we are left with an integral $`𝑑\stackrel{}{l}_1𝑑\stackrel{}{l}_2𝑑\stackrel{}{l}_3`$ in exactly three space dimensions that can be performed numerically.

In the numerical method, we keep the sum over cuts $`C`$ inside the integrations. We take care of the ultraviolet divergences by simple renormalization subtractions on the integrand. We make certain deformations on the integration contours so as to keep away from poles of the form $`1/[E_FE_I\pm iϵ]`$, where $`E_F`$ is the energy of the final state and $`E_I`$ is the energy of an intermediate state. Then the integrals are all convergent and we calculate them by Monte Carlo numerical integration.

Let us now look at the contour deformation in a little more detail. We denote the momenta $`\{\stackrel{}{l}_1,\stackrel{}{l}_2,\stackrel{}{l}_3\}`$ collectively by $`l`$ whenever we do not need a more detailed description. Thus

$$^{(2)}=\underset{G}{}𝑑l\underset{C}{}g(G,C;l).$$
(12)
For cuts $`C`$ that leave a virtual loop integration, there are singularities in the integrand of the form $`E_FE_I+iϵ`$ (or $`E_FE_Iiϵ`$ if the loop is in the complex conjugate amplitude to the right of the cut). Here $`E_F`$ is the energy of the final state defined by the cut $`C`$ and $`E_I`$ is the energy of a possible intermediate state. For the purpose of this review, all we need to know is that $`E_FE_I=0`$ on a surface in the space of $`\stackrel{}{l}_1`$ for fixed $`\stackrel{}{l}_2`$ and $`\stackrel{}{l}_3`$ if we pick $`\stackrel{}{l}_1`$ to be the momentum that flows around the virtual loop. These singularities do not create divergences. The Feynman rules provide us with the $`iϵ`$ prescriptions that tell us what to do about the singularities: we should deform the integration contour into the complex $`\stackrel{}{l}_1`$ space so as to keep away from them. Thus we write our integral in the form

$$^{(2)}=\underset{G}{}𝑑l\underset{C}{}𝒥(G,C;l)g(G,C;l+i\kappa (G,C;l)).$$
(13)
Here $`i\kappa `$ is a purely imaginary nine-dimensional vector that we add to the real nine-dimensional vector $`l`$ to make a complex nine-dimensional vector. The imaginary part $`\kappa `$ depends on the real part $`l`$, so that when we integrate over $`l`$, the complex vector $`l+i\kappa `$ lies on a surface, the integration contour, that is moved away from the real subspace. When we thus deform the contour, we supply a jacobian $`𝒥=det((l+i\kappa )/l)`$. (See Ref. beowulfPRD  for details.)

The amount of deformation $`\kappa `$ depends on the graph $`G`$ and, more significantly, the cut $`C`$. For cuts $`C`$ that leave no virtual loop, each of the momenta $`\stackrel{}{l}_1`$, $`\stackrel{}{l}_2`$, and $`\stackrel{}{l}_3`$ flows through the final state. For practical reasons, we want the final state momenta to be real. Thus we set $`\kappa =0`$ for cuts $`C`$ that leave no virtual loop. On the other hand, when the cut $`C`$ does leave a virtual loop, we choose a non-zero $`\kappa `$. We must, however, be careful. When $`\kappa =0`$ there are singularities in $`g`$ on certain surfaces that correspond to collinear parton momenta. These singularities cancel between $`g`$ for one cut $`C`$ and $`g`$ for another. This cancellation would be destroyed if, for $`l`$ approaching the collinear singularity, $`\kappa =0`$ for one of these cuts but not for the other. For this reason, we insist that for all cuts $`C`$, $`\kappa 0`$ as $`l`$ approaches one of the collinear singularities. The details can be found in Ref. beowulfPRD . All that is important here is that $`\kappa 0`$ quadratically with the distance to a collinear singularity.

Much has been left out in this brief overview, but we should now have enough background to see what to do in Coulomb gauge. It might seem that all we have to do is use the Feynman rules in Coulomb gauge, but there are some questions associated with the two and three point subdiagrams that need a special analysis. We now turn to that analysis, beginning with a description of a momentum space coordinate system that is useful for the two point subdiagrams.

## III Elliptical coordinates

Most of our effort will be devoted to the gluon and quark one loop self-energy diagrams. In these diagrams, we begin by performing the integration over the loop energy by contour integration. This leaves an integration over the loop three-momentum. We will be greatly helped by expressing the components of the loop momentum in terms of three appropriately chosen variables $`\{\mathrm{\Delta },x,\varphi \}`$. These variables are defined by considering, instead of the virtual self-energy graph, the corresponding cut self-energy graph.

Consider an off-shell parton carrying momentum $`\overline{q}^\mu `$ that splits into two partons that carry momenta $`k_\pm ^\mu `$. Let these two partons be on-shell, $`k_\pm ^2=0`$. Thus $`k_\pm ^0=\omega _\pm `$ where

$$\omega _\pm |\stackrel{}{k}_\pm |.$$
(14)
We consider the space-part of $`\overline{q}^\mu `$ to be fixed, and we call it $`\stackrel{}{q}`$, while the energy is determined by energy conservation

$$\overline{q}=(\omega _++\omega _{},\stackrel{}{q}).$$
(15)
We define a loop momentum $`\stackrel{}{l}`$ by

$$\stackrel{}{k}_\pm ={\scriptscriptstyle \frac{1}{2}}\stackrel{}{q}\pm \stackrel{}{l}.$$
(16)
We will use elliptical coordinates $`\{\mathrm{\Delta },x,\varphi \}`$ defined as follows. First, define $`𝒬`$ by

$$𝒬=|\stackrel{}{q}|.$$
(17)
Now, define coordinates $`\{\mathrm{\Delta },x\}`$ by

$`\mathrm{\Delta }+1`$ $`=`$ $`{\displaystyle \frac{1}{𝒬}}\left(\omega _++\omega _{}\right),`$
$`2x1`$ $`=`$ $`{\displaystyle \frac{1}{𝒬}}\left(\omega _+\omega _{}\right).`$ (18)
Then

$$0<\mathrm{\Delta },0<x<1.$$
(19)
Finally, let $`\varphi `$ be the azimuthal angle of $`\stackrel{}{l}`$ in a coordinate system in which the $`z`$-axis lies in the direction of $`\stackrel{}{q}`$ and the direction of the $`x`$-axis is defined arbitrarily. Often, we will want to work in $`32ϵ`$ space dimensions. In this case, $`\varphi `$ stands for a point on the surface of a unit sphere in $`22ϵ`$ dimensions, with

$$𝑑\varphi S(22ϵ)=\frac{2\pi ^{1ϵ}}{\mathrm{\Gamma }(1ϵ)}.$$
(20)
The surfaces of constant $`\mathrm{\Delta }`$ are ellipsoids, while the surfaces of constant $`x`$ are paraboloids that are orthogonal to the constant $`\mathrm{\Delta }`$ surfaces. Both of these surfaces are orthogonal to the surfaces of constant $`\varphi `$.

We note immediately that

$$\overline{q}^2=(\omega _++\omega _{})^2𝒬^2=𝒬^2\left[(\mathrm{\Delta }+1)^21\right]=𝒬^2\mathrm{\Delta }(\mathrm{\Delta }+2).$$
(21)
The inverse relation is

$$\mathrm{\Delta }=\sqrt{1+\overline{q}^2/𝒬^2}1.$$
(22)
It will often be convenient to use $`\overline{q}^2`$ as an independent variable instead of $`\mathrm{\Delta }`$.

The part, $`\stackrel{}{l}_T`$, of $`\stackrel{}{l}`$ transverse to $`\stackrel{}{q}`$ is determined by a unit vector in $`22ϵ`$ dimensions specified by $`\varphi `$ and by the magnitude $`|\stackrel{}{l}_T|`$, which is

$$|\stackrel{}{l}_T|=\sqrt{x(1x)\overline{q}^2}.$$
(23)
The component of $`\stackrel{}{l}`$ along $`\stackrel{}{q}`$ is

$$\stackrel{}{l}\stackrel{}{q}/𝒬=(𝒬/2)(1+\mathrm{\Delta })(2x1).$$
(24)
A straightforward calculation shows that the jacobian of the transformation is given by

$$\frac{d^{32ϵ}\stackrel{}{l}}{2\omega _+2\omega _{}}=\frac{1}{8𝒬}\frac{1}{1+\mathrm{\Delta }}[x(1x)\overline{q}^2]^ϵd\overline{q}^2dxd^{12ϵ}\varphi .$$
(25)
With $`ϵ=0`$, the transformation from $`\{\mathrm{\Delta },x,\varphi \}`$ to $`\stackrel{}{l}`$ is

$$\stackrel{}{l}=|\stackrel{}{l}_T|\mathrm{cos}\varphi \stackrel{}{n}_x+|\stackrel{}{l}_T|\mathrm{sin}\varphi \stackrel{}{n}_y+(\stackrel{}{l}\stackrel{}{q}/𝒬)\stackrel{}{n}_z.$$
(26)
Here $`\stackrel{}{n}_z=\stackrel{}{q}/𝒬`$, while $`\stackrel{}{n}_x`$ and $`\stackrel{}{n}_y`$ are two unit vectors orthogonal to each other in the plane orthogonal to $`n_z`$.

## IV Structure of graphs with a cut gluon propagator

In this and the following sections, we analyze the gluon propagator. We denote by $`n^\mu `$ a unit vector in the time direction, $`n=(1,0,0,0)`$. We use Coulomb gauge. Information about the use of Coulomb gauge can be found in Refs. ChristLee  and leibrandt . We consider diagrams in which there is zero or one loop in the gluon propagator. Thus we deal with the gluon propagator at orders $`\alpha _s^0`$ and $`\alpha _s^1`$.

We will be interested in the factors in the cross section that arise from a cut gluon propagator when the cut propagator is a subgraph of a larger graph. To set the notation, we write the contribution from an order $`\alpha _s^0`$ cut propagator, illustrated in Fig. 2, as

$$[\mathrm{Born}]=𝑑\stackrel{}{q}\frac{D(q)^{\mu \nu }}{2𝒬}R_{\mu \nu }^{(0)},$$
(27)
where $`\stackrel{}{q}`$ is the three-momentum carried by the propagator, $`𝒬|\stackrel{}{q}|`$, and $`q=(𝒬,\stackrel{}{q})`$. Then $`d\stackrel{}{q}/[2𝒬]`$ is the standard Lorentz invariant integration over the gluon mass shell. The tensor $`R_{\mu \nu }^{(0)}`$ denotes the factors associated with the rest of the graph and with the final state measurement function $`𝒮`$; $`R_{\mu \nu }^{(0)}`$ depends on $`\stackrel{}{q}`$, but this dependence is suppressed in the notation. Finally, $`D^{\mu \nu }`$ is the numerator function for a bare gluon propagator in Coulomb gauge, which is (for both on-shell and off-shell gluons),

$$D(k)^{\mu \nu }=g^{\mu \nu }+\frac{1}{\omega ^2}\left[k^\mu \stackrel{~}{k}^\nu \stackrel{~}{k}^\mu k^\nu +k^\mu k^\nu \right],$$
(28)
where

$$\stackrel{~}{k}=(0,\stackrel{}{k})$$
(29)
consists of just the $`\mu =1,2`$ and 3 components of $`k^\mu `$ and where, as in the previous section,

$$\omega =|\stackrel{}{k}|=\sqrt{\stackrel{~}{k}^2}.$$
(30)
At order $`\alpha _s^1`$, we consider a gluon propagator with a one loop self-energy subgraph. We consider three cases. In the first case, the two bare propagators in the self-energy subgraph are cut and the neighboring bare propagators are virtual. This case is illustrated in Fig. 3. We write the contribution of a cut gluon self-energy graph in the form

$$[\mathrm{real}]=𝑑\stackrel{}{q}\frac{1}{4\pi 𝒬}\frac{d\overline{q}^2}{\overline{q}^2}𝑑x𝑑\varphi _g^{\mu \nu }(\overline{q}^2,x,\varphi )R_{\mu \nu }(\overline{q}^2,x,\varphi ),$$
(31)
ignoring the infrared divergence, which will cancel inside the integral after combining real and virtual contributions in Eq. (IV) below. The integration over the loop momentum $`\stackrel{}{l}`$ has been changed to an integration over $`\{\overline{q}^2,x,\varphi \}`$ and the factors coming from the cut self-energy graph and its adjacent virtual gluon propagators are included in $`/(4\pi 𝒬\overline{q}^2)`$. Then $`R_{\mu \nu }`$ represents the rest of the graph, including the measurement function $`𝒮`$. All of the factors in $`R_{\mu \nu }`$ depend on the virtuality $`\overline{q}^2`$. The measurement function depends also on $`x`$ and $`\varphi `$. The dependence of $`R_{\mu \nu }`$ on $`\stackrel{}{q}`$ is suppressed in the notation. The conventions chosen are such that in Eq. (27) for a cut bare propagator, $`R_{\mu \nu }^{(0)}`$ is $`R_{\mu \nu }(\overline{q}^2,x,\varphi )`$ evaluated with $`\overline{q}^2=0`$. Here we use the infrared safety property (2) of the measurement function.

In the second case to be considered below the self-energy loop is entirely virtual and the neighboring bare propagators are not cut, so that the incoming momentum $`q`$ is off-shell. This case is illustrated in Fig. 4. We consider $`q^2<0`$. In this case, we investigate the quantity

$$F_g^{\mu \nu }(q)=\frac{1}{q^2}D(q)_\alpha ^\mu \mathrm{\Pi }(q)^{\alpha \beta }D(q)_\beta ^\nu ,$$
(32)
where $`\mathrm{\Pi }^{\alpha \beta }`$ is the self-energy function, Eq. (44). There are two factors of $`1/q^2`$ in the propagator function; we include just one of them in the definition of $`F_g`$. We find that $`F_g^{\mu \nu }(q)`$ is given by an integral of the form

$$F_g^{\mu \nu }(q)=\frac{1}{2\pi }𝑑\overline{q}^2𝑑x𝑑\varphi \frac{1}{\overline{q}^2q^2}W_g^{\mu \nu }(\overline{q}^2,x,\varphi ).$$
(33)
The   MS  ultraviolet renormalization is expressed through certain subtraction terms included in $`W`$.

In the third case, the self-energy loop is entirely virtual and one of the neighboring bare propagators is cut, so that the incoming momentum satisfies $`q^2=0`$. This case is illustrated in Fig. 5. Then we need $`F_g^{\mu \nu }(q)`$ with $`q^2=0`$, which takes the form

$$F_g^{\mu \nu }(q)=\frac{1}{2\pi }\frac{d\overline{q}^2}{\overline{q}^2}𝑑x𝑑\varphi 𝒲_g^{\mu \nu }(\overline{q}^2,x,\varphi ),$$
(34)
where $`𝒲_g`$ is $`W_g`$ with $`q^2=0`$. Again, the integral is infrared divergent, but the divergence will cancel inside the integrand in Eq. (IV). The corresponding contribution to the cross section $``$ is

$`[\mathrm{virtual}]`$ $`=`$ $`{\displaystyle 𝑑\stackrel{}{q}\frac{1}{2𝒬}F_g^{\mu \nu }(q)R_{\mu \nu }^{(0)}}`$
$`=`$ $`{\displaystyle 𝑑\stackrel{}{q}\frac{1}{4\pi 𝒬}\frac{d\overline{q}^2}{\overline{q}^2}𝑑x𝑑\varphi 𝒲_g^{\mu \nu }(\overline{q}^2,x,\varphi )R_{\mu \nu }^{(0)}}.`$
We should note that $`[\mathrm{virtual}]`$ in Eq. (IV) is the total contribution from the two graphs in which one or the other of the neighboring bare propagators is cut. The contribution from either of these graphs is half of this.

When we add the contributions from graphs with a cut self-energy subdiagram and with a virtual self-energy subdiagram with a cut adjoining bare propagator, we get

$`[\mathrm{real}]+[\mathrm{virtual}]`$ $`=`$ $`{\displaystyle 𝑑\stackrel{}{q}\frac{1}{4\pi 𝒬}\frac{d\overline{q}^2}{\overline{q}^2}𝑑x𝑑\varphi }`$
$`\times \left\{_g^{\mu \nu }(\overline{q}^2,x,\varphi )R_{\mu \nu }(\overline{q}^2,x,\varphi )+𝒲_g^{\mu \nu }(\overline{q}^2,x,\varphi )R_{\mu \nu }^{(0)}\right\}.`$
As noted above, the integrals of the two terms separately would be infrared divergent and would not make sense by themselves. However, the $`\overline{q}^20`$ singularities in the integrand in Eq. (IV) cancel, so that the integral is finite and is suitable for calculation by numerical integration.

With this preparation, we are ready to turn to the calculations for the three cases.

## V Real gluon self-energy graph

It is straightforward to evaluate the function $`_g^{\mu \nu }`$ defined in Eq. (31). We find

$`_g^{\mu \nu }`$ $`=`$ $`{\displaystyle \frac{\alpha _s}{4\pi }}{\displaystyle \frac{1}{1+\mathrm{\Delta }}}`$ (37)
$`\times \{N_{TT}D(𝒬,\stackrel{}{q})^{\mu \nu }+N_{tt}[{\displaystyle \frac{l_T^\mu l_T^\nu }{\stackrel{}{l}_T^{\mathrm{\hspace{0.17em}2}}}}\frac{1}{2}D(𝒬,\stackrel{}{q})^{\mu \nu }]`$
$`+{\displaystyle \frac{\overline{q}^2}{𝒬^2}}N_{EE}n^\mu n^\nu +N_{Et}{\displaystyle \frac{\overline{q}n}{(1+\mathrm{\Delta })𝒬^2}}(l_T^\mu n^\nu +n^\mu l_T^\nu )\}.`$
The coefficients $`N_{TT},N_{tt},N_{EE},`$ and $`N_{Et}`$ are functions of $`\overline{q}^2`$ and $`x`$ and are given below. The tensor $`D(𝒬,\stackrel{}{q})^{\mu \nu }`$ is the numerator (28) for an on-shell gluon, with momentum $`q_{\mathrm{os}}=(𝒬,\stackrel{}{q})`$. In a reference frame in which $`\stackrel{}{q}`$ is aligned with the $`z`$-axis, the only non-zero components of $`D(𝒬,\stackrel{}{q})^{\mu \nu }`$ are $`D(𝒬,\stackrel{}{q})^{ij}=\delta ^{ij}`$ with $`i,j\{1,2\}`$. In the second term, using the notation of Sec. III, $`l_T^\mu `$ is the vector

$`l_T^0`$ $`=`$ $`0,`$
$`\stackrel{}{l}_T`$ $`=`$ $`\stackrel{}{l}{\displaystyle \frac{\stackrel{}{l}\stackrel{}{q}}{𝒬^2}}\stackrel{}{q},`$ (38)
or

$$l_T^\mu =\sqrt{x(1x)\overline{q}^2}\left(\mathrm{cos}\varphi n_x^\mu +\mathrm{sin}\varphi n_y^\mu \right).$$
(39)
Note that this term vanishes if we average over angles $`\varphi `$. The third term gives $`^{00}`$. It vanishes when $`\overline{q}^20`$. Finally, the fourth term gives $`^{0i}`$ and $`^{i0}`$. Note that it vanishes if we average over angles $`\varphi `$.

We have written Eq. (37) in a more elaborate form than might have seemed necessary: $`\overline{q}n=(1+\mathrm{\Delta })Q`$, so the coefficient in the fourth term could have been simplified. The reason for the more elaborate form is as follows. The tensor $`_g^{\mu \nu }`$ is a function of the momenta $`\overline{q}^\mu ,k_+^\mu ,k_{}^\mu `$ carried on the lines of the graph. In the derivation, we have understood that gluons with momenta $`k_+^\mu `$ and $`k_{}^\mu `$ enter the final state. However in a calculation, it might be convenient to use momenta with reversed signs: $`\overline{q}^\mu \overline{q}^\mu ,k_+^\mu k_+^\mu ,k_{}^\mu k_{}^\mu `$. Then the vector $`l_T^\mu `$ is also reversed, while $`𝒬`$, $`\mathrm{\Delta }`$, and $`x`$ are unchanged. We have written $`_g^{\mu \nu }`$ in Eq. (37) in a form that is unchanged under this reversal of momenta.

The coefficients are

$`N_{TT}`$ $`=`$ $`2C_A\{x(1x)+{\displaystyle \frac{8x(1x)[1x(1x)]}{\overline{q}^2/𝒬^2+4x(1x)}}`$
$`+{\displaystyle \frac{16x(1x)[14x(1x)+2x^2(1x)^2]}{[\overline{q}^2/𝒬^2+4x(1x)]^2}}\}`$
$`+N_F\left\{12x(1x)\right\},`$
$`N_{tt}`$ $`=`$ $`4C_A\left\{x(1x){\displaystyle \frac{8x^2(1x)^2}{\overline{q}^2/𝒬^2+4x(1x)}}+{\displaystyle \frac{32x^3(1x)^3}{[\overline{q}^2/𝒬^2+4x(1x)]^2}}\right\}`$
$`4N_Fx(1x),`$
$`N_{EE}`$ $`=`$ $`C_A\{[14x(1x)]{\displaystyle \frac{8x(1x)[14x(1x)]}{\overline{q}^2/𝒬^2+4x(1x)}}`$
$`+{\displaystyle \frac{32x^2(1x)^2[14x(1x)]}{[\overline{q}^2/𝒬^2+4x(1x)]^2}}\}`$
$`+4N_Fx(1x),`$
$`N_{Et}`$ $`=`$ $`2C_A(2x1)\left\{1{\displaystyle \frac{2[14x(1x)]}{\overline{q}^2/𝒬^2+4x(1x)}}+{\displaystyle \frac{16x(1x)[12x(1x)]}{[\overline{q}^2/𝒬^2+4x(1x)]^2}}\right\}`$ (40)
$`+2N_F(2x1).`$
(In this paper, we use the standard notation in which $`N_F`$ is the number of quark flavors, $`N_C`$ is the number of colors, $`C_A=N_C`$, and $`C_F=(N_C^21)/(2N_C)`$.)

The behavior of $`^{\mu \nu }`$ as $`\overline{q}^20`$ is important. Leaving out the $`N_{tt}`$ and $`N_{Et}`$ terms, which average to zero after integrating over angles, we obtain

$$^{\mu \nu }D(q)^{\mu \nu }\frac{\alpha _s}{2\pi }\left\{\frac{1}{2}\stackrel{~}{P}_{g/g}(x)+N_F\stackrel{~}{P}_{q/g}(x)\right\},$$
(41)
where

$`\stackrel{~}{P}_{q/g}(x)`$ $`=`$ $`\frac{1}{2}[12x(1x)],`$
$`\stackrel{~}{P}_{g/g}(x)`$ $`=`$ $`2C_A{\displaystyle \frac{[1x(1x)]^2}{x(1x)}}`$ (42)
are the one loop parton evolution kernels without their $`x1`$ regulation. Notice that $`\stackrel{~}{P}_{g/g}(x)`$ is singular at $`x0`$ and $`x1`$ but that $`^{\mu \nu }`$ is not singular at these points. The singularities emerge only when we take the $`\overline{q}^20`$ limit. Also notice that in front of $`\stackrel{~}{P}_{g/g}(x)`$ in Eq. (41) there is a symmetry factor $`1/2`$, which arises because the two final state gluons in $`ggg`$ are identical.

## VI Virtual gluon self-energy graph

In this section, we analyze the virtual gluon self-energy graph at spacelike momentum $`q^\mu `$. We consider the self-energy function $`\mathrm{\Pi }(q)^{\mu \nu }`$. Later, we also consider the quantity

$$F_g^{\mu \nu }(q)=\frac{1}{q^2}D(q)_\mu ^{}^\mu \mathrm{\Pi }(q)^{\mu ^{}\nu ^{}}D(q)_\nu ^{}^\nu ,$$
(43)
where $`D(q)_\mu ^{}^\mu `$ is the numerator (28) of the bare gluon propagator. In applications with a cut gluon propagator, as in Fig. 5, one needs $`F_g^{\mu \nu }(q)`$ evaluated with $`q^2=0`$.

We let the momenta of the two partons in the loop be $`k_\pm ^\mu `$, with $`k_+^\mu +k_{}^\mu =q^\mu `$. Then the Feynman rules for $`\mathrm{\Pi }(q)^{\mu \nu }`$ give

$`\mathrm{\Pi }(q)^{\mu \nu }`$ $`=`$ $`ig^2\stackrel{~}{\mu }^{2ϵ}{\displaystyle \frac{d^{42ϵ}k_+}{(2\pi )^{42ϵ}}\frac{1}{(k_+^2+iϵ)(k_{}^2+iϵ)}}`$ (44)
$`\times \{{\displaystyle \frac{C_A}{2}}V^{\mu \alpha \beta }(q,k_+,k_{})D_{\alpha \alpha ^{}}(k_+)D_{\beta \beta ^{}}(k_{})V^{\alpha ^{}\beta ^{}\nu }(k_+,k_{},q)`$
$`+{\displaystyle \frac{C_A}{2}}(\stackrel{~}{k}_+^\mu \stackrel{~}{k}_{}^\nu +\stackrel{~}{k}_{}^\mu \stackrel{~}{k}_+^\nu ){\displaystyle \frac{k_+^2k_{}^2}{\omega _+^2\omega _{}^2}}`$
$`+{\displaystyle \frac{1}{2}}N_F\mathrm{Tr}\left[\gamma ^\mu \text{/}k_{}\gamma ^\nu \text{/}k_+\right]\}`$
$`(\mathrm{pole}),`$
where

$$V^{\alpha \beta \gamma }(k_A,k_B,k_C)=g^{\alpha \beta }(k_A^\gamma k_B^\gamma )+g^{\beta \gamma }(k_B^\alpha k_C^\alpha )+g^{\gamma \alpha }(k_C^\beta k_A^\beta ).$$
(45)
Here the first line is for the gluon loop, the second line is for the ghost loop, symmetrized over the two possible momentum labellings, and the third line is for the quark loop. The gluon loop includes a symmetry factor $`1/2`$. The terms corresponding to gluon and ghost loops contain group factors $`C_A=3`$, while the term corresponding to a quark loop has a factor $`N_F`$ and a color factor $`T_F=1/2`$. We subtract the ultraviolet pole, as required by the   MS  prescription. (The parameter $`\stackrel{~}{\mu }`$ is related to the   MS  scale $`\mu `$ by $`\mu ^2=4\pi \stackrel{~}{\mu }^2e^\gamma `$.)

We choose a coordinate system in which the $`z`$-axis is aligned with $`\stackrel{}{q}`$. Then $`q=(q^0,0,0,𝒬)`$. We assume that $`q^2<0`$. Later we will take the limit $`q^20`$.

### VI.1 The energy integral

We begin by performing the $`k_+^0`$ integral. Consider a term in $`\mathrm{\Pi }^{\mu \nu }`$ of the form

$$ig^2\stackrel{~}{\mu }^{2ϵ}\frac{d^{42ϵ}k_+}{(2\pi )^{42ϵ}}\frac{Af(k_+,k_{})}{(k_+^2+iϵ)(k_{}^2+iϵ)},$$
(46)
where $`A`$ is a function of $`\stackrel{}{k}_+`$ and $`\stackrel{}{k}_{}`$ but is independent of $`k_+^0`$ and $`k_{}^0`$. We perform the $`k_+^0`$ integration, leaving an integral over $`\stackrel{}{k}_+`$. Then we change variables to the elliptical coordinates $`\{\overline{q}^2,x,\varphi \}`$ defined from $`\stackrel{}{k}_\pm `$ as in Sec. III, so that our term takes the form

$$\frac{\alpha _s}{8\pi ^2}𝑑\overline{q}^2𝑑xd^{12ϵ}\varphi \left[\frac{x(1x)\overline{q}^2}{4\pi ^2\stackrel{~}{\mu }^2}\right]^ϵ\frac{Ag(\overline{q}^2,x,\varphi )}{\overline{q}^2q^2}.$$
(47)
Thus we will need an integration table that translates $`f`$ into $`g`$.

There is, however, a problem. Some of the integrals over $`k_+^0`$ are divergent. Thus we need a definition. We elect to perform the integrals over loop energies inside the integrals over loop three-momenta. We calculate these integrals by closing the energy contours in the lower half plane and then calculate them again by closing the contours in the upper half plane. Finally, we average the two results. The radius of the large semicircles that close the contours are always to be big enough to enclose all poles. Thus our prescription is a simple algebraic prescription of adding, with the appropriate signs and a factor 1/2, the residues of all the poles in the complex energy plane. If we apply this prescription in Feynman gauge, then the energy integrations are all convergent and we get the usual answer. In Coulomb gauge, we make this prescription part of the definition of the gauge. The required integral table is given as Table 1.

Of these integrals, the most divergent is that for $`f=k_+^2k_{}^2`$. There are no poles in the complex loop-energy plane, so we have

$$𝑑E=0.$$
(48)
One might expect that if the self-energy diagram is embedded in a larger diagram as part of a gauge invariant calculation, then an $`𝑑E`$ contribution from some other virtual loop diagram would cancel this contribution. However, so far as we can see, this does not happen. Specifically, the gluon-quark-antiquark one loop three point function does not have an $`𝑑E`$ divergence. Instead, we look to the coefficient of $`𝑑E`$. This coefficient is independent of $`q^2`$ and is proportional to the integral

$$\stackrel{~}{\mu }^{2ϵ}\frac{d^{32ϵ}\stackrel{}{k}_+}{(2\pi )^{32ϵ}}\left\{\frac{1}{\stackrel{}{k}_+^2}+\frac{1}{\stackrel{}{k}_{}^2}\right\}.$$
(49)
This integral vanishes in dimensional regulation for any $`ϵ`$. Thus we have an ambiguous contribution of infinity times zero and we find it sensible that the contour integration prescription given above instructs us to discard this contribution.

The integral for $`f=(k_+^2)^2(k_{}^2)^2`$ is logarithmically divergent and can be obtained from the simple integral

$$𝑑E\frac{E}{E^2M^2+iϵ}=0.$$
(50)
Zero is a sensible result for this integral since the integrand is odd under $`EE`$. In our prescription, we get $`i\pi `$ when we close the contour in the lower half plane and $`+i\pi `$ when we close the contour in the upper half plane. Averaging these two results give zero.

Our prescription may be compared to that of Leibbrandt and Williams leibrandt . In that prescription, one imagines that the loop energy is integrated over the imaginary axis and that the integration is performed in $`12\beta `$ dimensions. This prescription gives the same result as the one adopted here: the integral in Eq. (48) becomes $`d^{12\beta }\stackrel{}{k}`$, which vanishes because it has no scale, while the integral in Eq. (50) becomes $`d^{12\beta }\stackrel{}{k}\stackrel{}{k}\stackrel{}{n}/(\stackrel{}{k}^2+M^2+iϵ)`$, which vanishes because it is odd under $`\stackrel{}{k}\stackrel{}{n}\stackrel{}{k}\stackrel{}{n}`$.

### VI.2 Components needed

We will calculate the individual components of $`\mathrm{\Pi }(q)^{\mu \nu }`$ that we need in the coordinate frame with $`q=(q^0,0,0,𝒬)`$. Since $`D(q)_3^\mu =D(q)_3^\nu =0`$, we need not consider the components $`\mathrm{\Pi }(q)^{3\nu }`$ or $`\mathrm{\Pi }(q)^{\mu 3}`$. In addition, $`\mathrm{\Pi }(q)^{0i}=\mathrm{\Pi }(q)^{i0}=0`$ for $`i\{1,2\}`$ because of rotational symmetry. Thus what we need are $`\mathrm{\Pi }(q)^{00}`$ and $`\mathrm{\Pi }(q)^{ij}`$ for $`i,j\{1,2\}`$. We define

$`\mathrm{\Pi }(q)^{ij}`$ $`=`$ $`q^2A_T(q^2)\delta ^{ij},i,j\{1,2\},`$
$`\mathrm{\Pi }(q)^{00}`$ $`=`$ $`𝒬^2A_E(q^2).`$ (51)
We turn first to $`A_T(q)`$.

### VI.3 Transverse components

We evaluate $`A_T(q^2)`$ using the integrals in Table 1 and find

$`A_T(q^2)`$ $`=`$ $`{\displaystyle \frac{\alpha _s}{8\pi ^2}}{\displaystyle 𝑑\overline{q}^2𝑑xd^{12ϵ}\varphi \left[\frac{x(1x)\overline{q}^2}{4\pi ^2\stackrel{~}{\mu }^2}\right]^ϵ\frac{1}{\overline{q}^2q^2}\underset{J=0}{\overset{2}{}}\frac{A_{T,J}^{\prime \prime }}{[\overline{q}^2/𝒬^2+4x(1x)]^J}}`$ (52)
$`(\mathrm{pole}),`$
where

$`A_{T,0}^{\prime \prime }`$ $`=`$ $`C_A{\displaystyle \frac{1}{q^2}}\{\overline{q}^2/23q^2/24𝒬^2+2\overline{q}^2x(1x)+24𝒬^2x(1x)`$
$`+{\displaystyle \frac{x(1x)}{(1ϵ)}}[\overline{q}^2+4𝒬^2+3q^220𝒬^2x(1x)]\}`$
$`N_F{\displaystyle \frac{1}{q^2}}\left\{\overline{q}^2+{\displaystyle \frac{x(1x)}{(1ϵ)}}2\overline{q}^2\right\},`$
$`A_{T,1}^{\prime \prime }`$ $`=`$ $`C_A{\displaystyle \frac{1}{q^2}}\{4q^2+2(q^2)^2/𝒬^2+48𝒬^2x(1x)224𝒬^2[x(1x)]^2`$
$`+{\displaystyle \frac{x(1x)}{(1ϵ)}}[4q^24(q^2)^2/𝒬^232𝒬^2x(1x)`$
$`28q^2x(1x)+144𝒬^2[x(1x)]^2]\},`$
$`A_{T,2}^{\prime \prime }`$ $`=`$ $`C_A{\displaystyle \frac{1}{q^2}}\{128𝒬^2[x(1x)]^2+512𝒬^2[x(1x)]^3`$ (53)
$`+{\displaystyle \frac{[x(1x)]^2}{(1ϵ)}}[16q^2+16(q^2)^2/𝒬^2+64𝒬^2x(1x)`$
$`+64q^2x(1x)256𝒬^2[x(1x)]^2]\}.`$
As indicated by the notation, we expect that $`q^2A_T(q^2)`$ vanishes at $`q^2=0`$, but it is not evident from the form above that it does so. Nevertheless, explicit analytical integration shows that $`q^2A_T(q^2)=0`$ at $`q^2=0`$ for any $`ϵ`$. In order to get a form in which the vanishing of $`q^2A_T(q^2)`$ at $`q^2=0`$ is manifest, we simply subtract $`q^2A_T(q^2)`$ at $`q^2=0`$, which is zero, from $`q^2A_T(q^2)`$. Then $`A_T(q^2)`$ still has the form (52) but with new coefficients $`A_{T,J}^{}`$:

$`A_{T,0}^{}`$ $`=`$ $`2C_A\left\{1+x(1x)+{\displaystyle \frac{2x(1x)}{(1ϵ)}}\right\}+N_F\left\{1{\displaystyle \frac{2x(1x)}{(1ϵ)}}\right\},`$
$`A_{T,1}^{}`$ $`=`$ $`2C_A\left\{q^2/𝒬^2+12x(1x){\displaystyle \frac{2x(1x)}{(1ϵ)}}\left[q^2/𝒬^2+12x(1x)\right]\right\},`$
$`A_{T,2}^{}`$ $`=`$ $`16C_Ax(1x)\left\{28x(1x)+{\displaystyle \frac{x(1x)}{(1ϵ)}}\left[q^2/𝒬^2+8x(1x)\right]\right\}.`$ (54)
We now examine the ultraviolet renormalization of $`A_T(q^2)`$. Define

$`\mathrm{\Delta }A_T(q^2)`$ $`=`$ $`{\displaystyle \frac{\alpha _s}{8\pi ^2}}{\displaystyle 𝑑\overline{q}^2𝑑xd^{12ϵ}\varphi \left[\frac{x(1x)\overline{q}^2}{4\pi ^2\stackrel{~}{\mu }^2}\right]^ϵ}`$ (55)
$`\times \left\{{\displaystyle \frac{2C_A+N_F}{\overline{q}^2+e^{\lambda _1}\mu ^2}}+x(1x){\displaystyle \frac{2C_A}{\overline{q}^2+e^{\lambda _2}\mu ^2}}+{\displaystyle \frac{x(1x)}{1ϵ}}{\displaystyle \frac{4C_A2N_F}{\overline{q}^2+e^{\lambda _3}\mu ^2}}\right\}`$
$`(\mathrm{pole}).`$
Here $`\mu ^2=4\pi \stackrel{~}{\mu }^2e^\gamma `$ is the   MS  scale, $`\gamma =0.577\mathrm{}`$ is the Euler constant, and $`\lambda _1`$, $`\lambda _2`$, and $`\lambda _3`$ are parameters that we will adjust. The integrand in $`\mathrm{\Delta }A_T`$ matches that of $`A_T`$ when $`\overline{q}^2\mathrm{}`$ but it has the opposite sign, so if we add $`\mathrm{\Delta }A_T`$ to $`A_T`$, we will obtain an ultraviolet convergent integral. Furthermore, since the integrands match for $`\overline{q}^2\mathrm{}`$, the ultraviolet pole terms that are included in the definitions are opposite and will cancel in the sum $`A_T+\mathrm{\Delta }A_T`$. We can easily perform the integral:

$`\mathrm{\Delta }A_T(q^2)`$ $`=`$ $`{\displaystyle \frac{\alpha _s}{4\pi }}\mathrm{\Gamma }(ϵ)e^{\gamma ϵ}\{(2C_AN_F){\displaystyle \frac{\mathrm{\Gamma }(1ϵ)^2}{\mathrm{\Gamma }(22ϵ)}}e^{\lambda _1ϵ}+2C_A{\displaystyle \frac{\mathrm{\Gamma }(2ϵ)^2}{\mathrm{\Gamma }(42ϵ)}}e^{\lambda _2ϵ}`$ (56)
$`+(4C_A2N_F){\displaystyle \frac{\mathrm{\Gamma }(2ϵ)^2}{(1ϵ)\mathrm{\Gamma }(42ϵ)}}e^{\lambda _3ϵ}\}(\mathrm{pole})`$
$`=`$ $`{\displaystyle \frac{\alpha _s}{4\pi }}\left\{(2C_AN_F)\left(2\lambda _1\right)+{\displaystyle \frac{2C_A}{6}}\left({\displaystyle \frac{5}{3}}\lambda _2\right)+{\displaystyle \frac{4C_A2N_F}{6}}\left({\displaystyle \frac{8}{3}}\lambda _3\right)\right\}`$
$`+𝒪(ϵ)\text{}.`$
We set

$$\lambda _1=2,\lambda _2=\frac{5}{3},\lambda _3=\frac{8}{3}.$$
(57)
Then

$$\mathrm{\Delta }A_T(q^2)=0+𝒪(ϵ).$$
(58)
Since $`\mathrm{\Delta }A_T`$ is zero, we can add it to $`A_T`$ to obtain, after setting $`ϵ`$ to zero,

$$A_T(q^2)=\frac{\alpha _s}{8\pi ^2}𝑑\overline{q}^2𝑑x𝑑\varphi \frac{1}{\overline{q}^2q^2}\underset{J=0}{\overset{2}{}}\frac{A_{T,J}}{[\overline{q}^2/𝒬^2+4x(1x)]^J},$$
(59)
where

$`A_{T,0}`$ $`=`$ $`(2C_AN_F){\displaystyle \frac{q^2+e^2\mu ^2}{\overline{q}^2+e^2\mu ^2}}+2C_Ax(1x){\displaystyle \frac{q^2+e^{5/3}\mu ^2}{\overline{q}^2+e^{5/3}\mu ^2}}`$
$`+(4C_A2N_F)x(1x){\displaystyle \frac{q^2+e^{8/3}\mu ^2}{\overline{q}^2+e^{8/3}\mu ^2}},`$
$`A_{T,1}`$ $`=`$ $`2C_A\left\{q^2/𝒬^2+12x(1x)2x(1x)\left[q^2/𝒬^2+12x(1x)\right]\right\},`$
$`A_{T,2}`$ $`=`$ $`16C_Ax(1x)\left\{28x(1x)+x(1x)\left[q^2/𝒬^2+8x(1x)\right]\right\}.`$ (60)
### VI.4 Timelike components

We evaluate $`A_E(q^2)`$ using the integrals in Table 1 and find

$`A_E(q^2)`$ $`=`$ $`{\displaystyle \frac{\alpha _s}{8\pi ^2}}{\displaystyle 𝑑\overline{q}^2𝑑xd^{12ϵ}\varphi \left[\frac{x(1x)\overline{q}^2}{4\pi ^2\stackrel{~}{\mu }^2}\right]^ϵ\frac{1}{\overline{q}^2q^2}\underset{J=0}{\overset{2}{}}\frac{A_{E,J}^{\prime \prime }}{[\overline{q}^2/𝒬^2+4x(1x)]^J}}`$ (61)
$`(\mathrm{pole}),`$
where

$`A_{E,0}^{\prime \prime }`$ $`=`$ $`C_A\left\{24x(1x)+(1ϵ)\left[(q^2\overline{q}^2)/𝒬^2+14x(1x)\right]\right\}+4N_Fx(1x),`$
$`A_{E,1}^{\prime \prime }`$ $`=`$ $`8C_Ax(1x)\left\{1+3q^2/𝒬^2+28x(1x)\right\},`$
$`A_{E,2}^{\prime \prime }`$ $`=`$ $`32C_A[x(1x)]^2\left\{13q^2/𝒬^216x(1x)\right\}.`$ (62)
We note that $`A_{E,0}^{\prime \prime }`$ contains a term $`C_A(1ϵ)(q^2\overline{q}^2)/𝒬^2`$. This term is undesirable for numerical integration because it leads to a quadratically ultraviolet divergent integral. Furthermore, this term gives a contribution to the integral that vanishes for all $`ϵ`$. For these reasons, we eliminate the term. Then $`A_E(q^2)`$ still has the form (61) but with new coefficients $`A_{E,J}^{}`$:

$`A_{E,0}^{}`$ $`=`$ $`C_A\left\{24x(1x)+(1ϵ)\left[14x(1x)\right]\right\}+4N_Fx(1x),`$
$`A_{E,1}^{}`$ $`=`$ $`8C_Ax(1x)\left\{1+3q^2/𝒬^2+28x(1x)\right\},`$
$`A_{E,2}^{}`$ $`=`$ $`32C_A[x(1x)]^2\left\{13q^2/𝒬^216x(1x)\right\}.`$ (63)
We now examine the ultraviolet renormalization of $`A_E(q^2)`$. Define

$`\mathrm{\Delta }A_E(q^2)`$ $`=`$ $`{\displaystyle \frac{\alpha _s}{8\pi ^2}}{\displaystyle 𝑑\overline{q}^2𝑑xd^{12ϵ}\varphi \left[\frac{x(1x)\overline{q}^2}{4\pi ^2\stackrel{~}{\mu }^2}\right]^ϵ}`$ (64)
$`\times \left\{(1ϵ){\displaystyle \frac{C_A}{\overline{q}^2+e^{\lambda _4}\mu ^2}}x(1x){\displaystyle \frac{24C_A4N_F}{\overline{q}^2+e^{\lambda _5}\mu ^2}}x(1x)(1ϵ){\displaystyle \frac{4C_A}{\overline{q}^2+e^{\lambda _6}\mu ^2}}\right\}`$
$`(\mathrm{pole}).`$
Here as before $`\mu ^2=4\pi \stackrel{~}{\mu }^2e^\gamma `$ is the   MS  scale and $`\lambda _4`$, $`\lambda _5`$, and $`\lambda _6`$ are parameters that we will adjust. The integrand in $`\mathrm{\Delta }A_E`$ matches that of $`A_E`$ when $`\overline{q}^2\mathrm{}`$ but it has the opposite sign, so if we add $`\mathrm{\Delta }A_E`$ to $`A_E`$, we will obtain an ultraviolet convergent integral. Furthermore, since the integrands match for $`\overline{q}^2\mathrm{}`$, the ultraviolet pole terms that are included in the definitions are opposite and will cancel in the sum $`A_E+\mathrm{\Delta }A_E`$. We can easily perform the integration:

$`\mathrm{\Delta }A_E(q^2)`$ $`=`$ $`{\displaystyle \frac{\alpha _s}{4\pi }}\mathrm{\Gamma }(ϵ)e^{\gamma ϵ}\{C_A{\displaystyle \frac{(1ϵ)\mathrm{\Gamma }(1ϵ)^2}{\mathrm{\Gamma }(22ϵ)}}e^{\lambda _4ϵ}`$ (65)
$`(24C_A4N_F){\displaystyle \frac{\mathrm{\Gamma }(2ϵ)^2}{\mathrm{\Gamma }(42ϵ)}}e^{\lambda _5ϵ}4C_A{\displaystyle \frac{(1ϵ)\mathrm{\Gamma }(2ϵ)^2}{\mathrm{\Gamma }(42ϵ)}}e^{\lambda _6ϵ}\}(\mathrm{pole})`$
$`=`$ $`{\displaystyle \frac{\alpha _s}{4\pi }}\left\{C_A(1\lambda _4)(24C_A4N_F){\displaystyle \frac{1}{6}}({\displaystyle \frac{5}{3}}\lambda _5)4C_A{\displaystyle \frac{1}{6}}({\displaystyle \frac{2}{3}}\lambda _6)\right\}`$
$`+𝒪(ϵ).\text{}`$
We set

$$\lambda _4=1,\lambda _5=\frac{5}{3},\lambda _6=\frac{2}{3}.$$
(66)
Then

$$\mathrm{\Delta }A_E(q^2)=0+𝒪(ϵ).$$
(67)
Since $`\mathrm{\Delta }A_E`$ is zero, we can add it to $`A_E`$ to obtain, after setting $`ϵ`$ to zero,

$$A_E(q^2)=\frac{\alpha _s}{8\pi ^2}𝑑\overline{q}^2𝑑x𝑑\varphi \frac{1}{\overline{q}^2q^2}\underset{J=0}{\overset{2}{}}\frac{A_{E,J}}{[\overline{q}^2/𝒬^2+4x(1x)]^J},$$
(68)
where

$`A_{E,0}`$ $`=`$ $`C_A{\displaystyle \frac{q^2+e^1\mu ^2}{\overline{q}^2+e^1\mu ^2}}(24C_A4N_F)x(1x){\displaystyle \frac{q^2+e^{5/3}\mu ^2}{\overline{q}^2+e^{5/3}\mu ^2}}4C_Ax(1x){\displaystyle \frac{q^2+e^{2/3}\mu ^2}{\overline{q}^2+e^{2/3}\mu ^2}},`$
$`A_{E,1}`$ $`=`$ $`8C_Ax(1x)\left\{1+3q^2/𝒬^2+28x(1x)\right\},`$
$`A_{E,2}`$ $`=`$ $`32C_A[x(1x)]^2\left\{13q^2/𝒬^216x(1x)\right\}.`$ (69)
### VI.5 The tensor $`F_g^{\mu \nu }`$

We can assemble this information to obtain $`F_g^{\mu \nu }(q)`$ defined in Eq. (43). With $`\stackrel{}{q}`$ directed along the $`z`$ axis, the only non-zero components of $`F_g^{\mu \nu }`$ are $`F_g^{00}`$ and $`F_g^{ij}`$ for $`i,j\{1,2\}`$:

$`F_g^{ij}(q)`$ $`=`$ $`A_T(q^2)\delta ^{ij},i,j\{1,2\},`$
$`F_g^{00}(q)`$ $`=`$ $`{\displaystyle \frac{q^2}{𝒬^2}}A_E(q^2).`$ (70)
We can write this in a form valid for any direction of $`\stackrel{}{q}`$ as

$$F_g^{\mu \nu }(q)=D(q)^{\mu \nu }A_T(q^2)+n^\mu n^\nu \frac{q^2}{𝒬^2}\left[A_E(q^2)A_T(q^2)\right].$$
(71)
Here $`D(q)^{\mu \nu }`$ is the numerator for the bare gluon propagator, Eq. (28), with $`D^{00}=q^2/𝒬^2`$. Note that this result is physically sensible. The complete one loop contribution to the gluon propagator is $`F^{\mu \nu }(q)/q^2`$. It has a simple pole times logarithms at $`q^2=0`$. The pole multiplies a projection onto transverse polarizations.

### VI.6 $`W_g^{\mu \nu }`$ and its supplementary terms

From Eq. (71) we obtain the tensor $`W_g^{\mu \nu }`$ defined in Eq. (33),

$`W_g^{\mu \nu }`$ $`=`$ $`{\displaystyle \frac{\alpha _s}{4\pi }}D(q)^{\mu \nu }{\displaystyle \underset{J}{}}{\displaystyle \frac{A_{T,J}}{[\overline{q}^2/𝒬^2+4x(1x)]^J}}`$ (72)
$`{\displaystyle \frac{\alpha _s}{4\pi }}n^\mu n^\nu {\displaystyle \frac{q^2}{𝒬^2}}{\displaystyle \underset{J}{}}{\displaystyle \frac{A_{E,J}A_{T,J}}{[\overline{q}^2/𝒬^2+4x(1x)]^J}}.`$
We can now set $`q^2`$ to 0 in $`W_g^{\mu \nu }`$ to obtain $`𝒲_g^{\mu \nu }`$. We have

$$𝒲_g^{\mu \nu }[\mathrm{simple}]=\frac{\alpha _s}{2\pi }D(q)^{\mu \nu }𝒫_g(\overline{q}^2,x),$$
(73)
where

$$𝒫_g(\overline{q}^2,x)=\frac{1}{2}\underset{J}{}\frac{A_{T,J}}{[\overline{q}^2/𝒬^2+4x(1x)]^J}.$$
(74)
Here the coefficients $`A_{T,J}`$ from Eq. (60) are evaluated with $`q^2=0`$.

We see immediately that there will be a numerical problem for the cancellation of $`𝒲_g^{\mu \nu }`$ with $`_g^{\mu \nu }`$ at the small $`\overline{q}^2`$ endpoint of the integration (IV). In Eq. (37) for $`_g^{\mu \nu }`$, there are terms with four tensor structures. There is a good cancellation for the $`D(q)^{\mu \nu }`$ structure (as we will see below) and the $`n^\mu n^\nu `$ structure is not important because it multiplies a factor $`\overline{q}^2`$ in Eq. (37) for $`_g^{\mu \nu }`$. For the other two terms, involving $`l_T^\mu `$, the singularity in $`_g^{\mu \nu }`$ will vanish if we integrate over the angle $`\varphi `$ of $`l_T^\mu `$ before letting $`\overline{q}^2`$ become small. However, this is not suitable for a numerical integration. Therefore, we modify $`𝒲_g^{\mu \nu }`$ to

$`𝒲_g^{\mu \nu }`$ $`=`$ $`{\displaystyle \frac{\alpha _s}{2\pi }}D(q)^{\mu \nu }𝒫_g(\overline{q}^2,x)`$ (75)
$`{\displaystyle \frac{\alpha _s}{4\pi }}\left[{\displaystyle \frac{l_T^\mu l_T^\nu }{\stackrel{}{l}_T^{\mathrm{\hspace{0.17em}2}}}}\frac{1}{2}D(q)^{\mu \nu }\right]{\displaystyle \frac{N_{tt}}{(1+\overline{q}^2/𝒬^2)}}`$
$`{\displaystyle \frac{\alpha _s}{4\pi }}{\displaystyle \frac{qn}{(1+\mathrm{\Delta })𝒬^2}}\left(l_T^\mu n^\nu +n^\mu l_T^\nu \right){\displaystyle \frac{N_{Et}}{(1+\overline{q}^2/𝒬^2)}},`$
where $`N_{tt}`$ and $`N_{Et}`$ are functions of $`\overline{q}^2`$ and $`x`$ and are given in Eq. (40). The integral of the extra terms vanishes (because the integration over $`\varphi `$ gives zero), so we are adding zero to $`F_g^{\mu \nu }`$. However the cancellation in Eq. (IV) now works point by point in $`\{\overline{q}^2,x,\varphi \}`$ space. We have inserted factors $`1/(1+\overline{q}^2/𝒬^2)`$ in the extra terms so as not to create problems at $`\overline{q}^2\mathrm{}`$ at the same time as we were alleviating problems at $`\overline{q}^20`$. As in Eq. (37), we have written $`𝒲_g^{\mu \nu }`$ in Eq. (75) in a form that is invariant under the replacements $`q^\mu q^\mu ,l^\mu l^\mu `$.

Let us examine the cancellation for $`\overline{q}^20`$ in Eq. (IV). The terms with tensor structures involving $`l_T^\mu `$ cancel the corresponding terms in $`_g^{\mu \nu }`$ by construction. For the $`D(q)^{\mu \nu }`$ term, the $`\overline{q}^20`$ limit is

$$𝒲_g^{\mu \nu }\frac{\alpha _s}{2\pi }D(q)^{\mu \nu }𝒫_g(0,x).$$
(76)
We find for $`𝒫_g(\overline{q}^2,x)`$ at $`\overline{q}^2=0`$,

$$𝒫_g(0,x)=\frac{1}{2}\stackrel{~}{P}_{g/g}(x)+N_F\stackrel{~}{P}_{q/g}(x),$$
(77)
where the parton evolution kernels $`\stackrel{~}{P}_{g/g}(x)`$ and $`\stackrel{~}{P}_{q/g}(x)`$ are given in Eq. (42). Using Eq. (41), we see that this is just the behavior we needed to make the cancellation work.

## VII Structure of graphs with a cut quark propagator

In this and the following sections, we analyze the quark propagator in Coulomb gauge. To set the notation, we write the contribution from an order $`\alpha _s^0`$ cut propagator as

$$[\mathrm{Born}]=𝑑\stackrel{}{q}\mathrm{Tr}\left\{\frac{\text{/}q}{2𝒬}R_0\right\},$$
(78)
where, as in Sec. IV, $`R_0`$ denotes the factors associated with the rest of the graph and with the final state measurement function $`𝒮`$. The function $`R_0`$ carries hidden Dirac indices and there is a trace over the Dirac indices of $`\text{/}qR_0`$.

Following the notation employed for the gluon propagator, we write the contribution of a cut quark self-energy graph as

$$[\mathrm{real}]=𝑑\stackrel{}{q}\mathrm{Tr}\left\{\frac{1}{4\pi 𝒬}\frac{d\overline{q}^2}{\overline{q}^2}𝑑x𝑑\varphi _q(\overline{q}^2,x,\varphi )R(\overline{q}^2,x,\varphi )\right\}.$$
(79)
With this notation, the contribution from the virtual self-energy graph with an adjoining cut bare propagator is

$$[\mathrm{virtual}]=𝑑\stackrel{}{q}\mathrm{Tr}\left\{\frac{1}{2𝒬}F_q(q)R_0\right\},$$
(80)
where

$$F_q(q)=\frac{\text{/}q\mathrm{\Sigma }(q)\text{/}q}{q^2}.$$
(81)
We should note that $`[\mathrm{virtual}]`$ in Eq. (80) is the total contribution from the two graphs in which one or the other of the neighboring bare propagators is cut. The contribution from either of these graphs is half of this.

If both of the adjoining bare propagators are uncut, the corresponding expression is

$$[\mathrm{all}\mathrm{uncut}]=d^4q\mathrm{Tr}\left\{\frac{1}{q^2}F_q(q)R(q^2)\right\},$$
(82)
where, here, $`R`$ is defined so that the corresponding Born contribution is $`R(q^2)\text{/}q/q^2`$.

We investigate $`F_q(q)`$ in the following section. First, we take $`q^2<0`$. We write $`F_q(q)`$ as an integral

$$F_q(q)=\frac{1}{2\pi }𝑑\overline{q}^2𝑑x𝑑\varphi \frac{1}{\overline{q}^2q^2}W_q(\overline{q}^2,x,\varphi ).$$
(83)
Then the contribution to the graph when the adjoining bare propagators are uncut is

$$[\mathrm{all}\mathrm{uncut}]=d^4q\mathrm{Tr}\left\{\frac{1}{2\pi q^2}𝑑\overline{q}^2𝑑x𝑑\varphi \frac{1}{\overline{q}^2q^2}W_q(\overline{q}^2,x,\varphi )R(q^2)\right\}.$$
(84)
Now, taking $`q^2`$ to zero, we write $`F_q(q)`$ as

$$F_q(q)=\frac{1}{2\pi }\frac{d\overline{q}^2}{\overline{q}^2}𝑑x𝑑\varphi 𝒲_q(\overline{q}^2,x,\varphi ),$$
(85)
where $`𝒲_q`$ is $`W_q`$ with $`q^2=0`$. Then

$`[\mathrm{real}]+[\mathrm{virtual}]`$ $`=`$ $`{\displaystyle }d\stackrel{}{q}\mathrm{Tr}\{{\displaystyle \frac{1}{4\pi 𝒬}}{\displaystyle }{\displaystyle \frac{d\overline{q}^2}{\overline{q}^2}}dxd\varphi `$
$`\times [_q(\overline{q}^2,x,\varphi )R(\overline{q}^2,x,\varphi )+𝒲_q(\overline{q}^2,x,\varphi )R_0]\}.`$
As we shall see, the $`\overline{q}^20`$ singularities in the integrand in Eq. (VII) cancel, so that the integral is finite and is suitable for calculation by numerical integration.

## VIII Real quark self-energy graph

It is straightforward to evaluate the function $`_q`$ defined in Eq. (79). We find

$$_q=\frac{\alpha _s}{4\pi }\frac{1}{1+\mathrm{\Delta }}\left\{N_L\text{/}q_{\mathrm{os}}+N_E\mathrm{\Delta }nq_{\mathrm{os}}\text{/}n+N_t\text{/}l_T\right\}.$$
(87)
The coefficients $`N_L,N_E`$, and $`N_t`$ are functions of $`\overline{q}^2`$ and $`x`$ and are given below. The momentum in the first two terms is $`q_{\mathrm{os}}=(\overline{q}^0/(1+\mathrm{\Delta }),\stackrel{}{q})=((k_+^0+k_{}^0)/(1+\mathrm{\Delta }),\stackrel{}{q})`$ so that $`q_{\mathrm{os}}^2=0`$. In the third term, using the notation of Sec. III, $`l_T^\mu `$ is the part of the loop momentum orthogonal to $`n^\mu `$ and $`q_{\mathrm{os}}^\mu `$. Note that this term vanishes if we average over angles $`\varphi `$.

The coefficients are

$`N_L`$ $`=`$ $`C_F\{12x(1x)+(2x1)(2x+\mathrm{\Delta })`$
$`{\displaystyle \frac{16x(1x)(2x1)}{2x+\mathrm{\Delta }}}+{\displaystyle \frac{16x(1x)[12x(1x)]}{(2x+\mathrm{\Delta })^2}}\},`$
$`N_E`$ $`=`$ $`2C_F{\displaystyle \frac{(1x)(4x^2+\mathrm{\Delta }^2)}{(2x+\mathrm{\Delta })^2}},`$
$`N_t`$ $`=`$ $`2C_F\left\{1{\displaystyle \frac{2(2x1)}{2x+\mathrm{\Delta }}}{\displaystyle \frac{8x(1x)}{(2x+\mathrm{\Delta })^2}}\right\}.`$ (88)
The behavior of $`𝒲_q^{\mu \nu }`$ as $`\overline{q}^20`$ is important. Leaving out the $`\text{/}l_T`$ term since $`l_T0`$ as $`\overline{q}^20`$, we obtain

$$_q^{\mu \nu }\text{/}q\frac{\alpha _s}{2\pi }\stackrel{~}{P}_{g/q}(x),$$
(89)
where

$$\stackrel{~}{P}_{g/q}(x)=C_F\frac{1+(1x)^2}{x}$$
(90)
is the one loop parton evolution kernel for $`qg`$. Notice that $`\stackrel{~}{P}_{q/g}(x)`$ is singular at $`x0`$ but that $`_q`$ is not singular at this point. The singularity emerges only when we take the $`\overline{q}^20`$ limit.

For the integrand of a cut antiquark self-energy graph, we have $`_{\overline{q}}=_q`$. We note that $`_q`$ in Eq. (87) is odd under the interchange $`\overline{q}^\mu \overline{q}^\mu `$, $`k_+^\mu k_+^\mu `$, $`k_{}^\mu k_{}^\mu `$ (so that also $`l^\mu l^\mu `$). Thus for an antiquark, we can simply use Eq. (87) and reverse the momenta, so that $`\overline{q}^\mu `$ flows in the direction of the fermion arrow in the graph. We use the same principle for the analysis that follows of the virtual quark self-energy graph.

## IX Virtual quark self-energy graph

In this subsection, we analyze the virtual self-energy graph at space-like momentum $`q^\mu `$. We consider the quantity

$$F_q(q)=\frac{\text{/}q\mathrm{\Sigma }(q)\text{/}q}{q^2}.$$
(91)
The Feynman rules give

$$F_q(q)=\frac{i}{q^2}g^2C_F\stackrel{~}{\mu }^{2ϵ}\frac{d^{42ϵ}k_+}{(2\pi )^{42ϵ}}\frac{\text{/}q\gamma _\mu \text{/}k_{}\gamma _\nu \text{/}q}{(k_+^2+iϵ)(k_{}^2+iϵ)}D(k_+)^{\mu \nu }(\mathrm{pole}).$$
(92)
The function $`F_q(q)`$ must have the form

$$F_q(q)=\text{/}qB_L(q^2)+qn\text{/}n(q^2/𝒬^2)B_E(q^2).$$
(93)
The functions $`B_L`$ and $`B_E`$ can be extracted using

$`B_L(q^2)`$ $`=`$ $`{\displaystyle \frac{1}{4𝒬^2}}\mathrm{Tr}\left\{F_q(q)\text{/}\stackrel{~}{q}\right\},`$
$`B_E(q^2)`$ $`=`$ $`{\displaystyle \frac{1}{4qnq^2}}\mathrm{Tr}\left\{F_q(q)[q^2\text{/}nqn\text{/}q]\right\}.`$ (94)
### IX.1 Space-like part

We evaluate $`B_L(q^2)`$ using the integrals in Table 1 and dropping terms that are odd under $`x(1x)`$, which will integrate to zero. We find

$`B_L(q^2)`$ $`=`$ $`{\displaystyle \frac{\alpha _s}{8\pi ^2}}{\displaystyle 𝑑\overline{q}^2𝑑xd^{12ϵ}\varphi \left[\frac{x(1x)\overline{q}^2}{4\pi ^2\stackrel{~}{\mu }^2}\right]^ϵ\frac{1}{\overline{q}^2q^2}\underset{J=0}{\overset{2}{}}\frac{B_{L,J}^{\prime \prime }}{[\overline{q}^2/𝒬^2+4x(1x)]^J}}`$ (95)
$`(\mathrm{pole}),`$
where

$`B_{L,0}^{\prime \prime }`$ $`=`$ $`C_F\left\{4{\displaystyle \frac{𝒬^2}{q^2}}+24{\displaystyle \frac{𝒬^2}{q^2}}x(1x)1ϵ+12x(1x)\right\},`$
$`B_{L,1}^{\prime \prime }`$ $`=`$ $`2C_F\{24{\displaystyle \frac{𝒬^2}{q^2}}x(1x)112{\displaystyle \frac{𝒬^2}{q^2}}[x(1x)]^2+2+8x(1x)56[x(1x)]^2`$
$`+{\displaystyle \frac{q^2}{𝒬^2}}2{\displaystyle \frac{q^2}{𝒬^2}}x(1x)\},`$
$`B_{L,2}^{\prime \prime }`$ $`=`$ $`64C_F[x(1x)]^2\left\{2{\displaystyle \frac{𝒬^2}{q^2}}+8{\displaystyle \frac{𝒬^2}{q^2}}x(1x)1+4x(1x)\right\}.`$ (96)
We expect that $`B_L(q^2)`$ is finite at $`q^2=0`$ when $`ϵ`$ is small and negative, but it is not evident from the form above that this is so. Nevertheless, explicit analytical integration shows that $`q^2B_L(q^2)=0`$ at $`q^2=0`$ for any $`ϵ`$. In order to get a form in which the vanishing of $`q^2B_L(q^2)`$ at $`q^2=0`$ is manifest, we simply subtract $`q^2B_L(q^2)`$ at $`q^2=0`$, which is zero, from $`q^2B_L(q^2)`$. Then $`B_L(q^2)`$ still has the form (95) but with new coefficients $`B_{L,J}^{}`$:

$`B_{L,0}^{}`$ $`=`$ $`C_F\{1ϵ+12x(1x)\},`$
$`B_{L,1}^{}`$ $`=`$ $`2C_F\left\{20x(1x)56[x(1x)]^2+{\displaystyle \frac{q^2}{𝒬^2}}2{\displaystyle \frac{q^2}{𝒬^2}}x(1x)\right\},`$
$`B_{L,2}^{}`$ $`=`$ $`32C_Fx(1x)\left\{16x(1x)+8[x(1x)]^2\right\}.`$ (97)
We now examine the ultraviolet renormalization of $`B_L(q^2)`$. We replace the subtraction of the ultraviolet pole in $`42ϵ`$ dimensions by a modification of the integrand in four dimensions, just as we did in the case of the gluon propagator. The result is given by Eq. (95) with $`ϵ=0`$ and no pole term to subtract,

$$B_L(q^2)=\frac{\alpha _s}{8\pi ^2}𝑑\overline{q}^2𝑑x𝑑\varphi \frac{1}{\overline{q}^2q^2}\underset{J=0}{\overset{2}{}}\frac{B_{L,J}}{[\overline{q}^2/𝒬^2+4x(1x)]^J},$$
(98)
and with new coefficients $`B_{L,J}`$:

$`B_{L,0}`$ $`=`$ $`C_F\left\{{\displaystyle \frac{q^2+e^3\mu ^2}{\overline{q}^2+e^3\mu ^2}}+12x(1x){\displaystyle \frac{q^2+e^{5/3}\mu ^2}{\overline{q}^2+e^{5/3}\mu ^2}}\right\},`$
$`B_{L,1}`$ $`=`$ $`2C_F\left\{20x(1x)56[x(1x)]^2+{\displaystyle \frac{q^2}{𝒬^2}}[12x(1x)]\right\},`$
$`B_{L,2}`$ $`=`$ $`32C_Fx(1x)\left\{16x(1x)+8[x(1x)]^2\right\}.`$ (99)
### IX.2 Timelike part

We evaluate $`B_E(q^2)`$ using the integrals in Table 1 and dropping terms that are odd under $`x(1x)`$, which will integrate to zero. We find

$`B_E(q^2)`$ $`=`$ $`{\displaystyle \frac{\alpha _s}{8\pi ^2}}{\displaystyle 𝑑\overline{q}^2𝑑xd^{12ϵ}\varphi \left[\frac{x(1x)\overline{q}^2}{4\pi ^2\stackrel{~}{\mu }^2}\right]^ϵ\frac{1}{\overline{q}^2q^2}\underset{J=0}{\overset{2}{}}\frac{B_{E,J}^{}}{[\overline{q}^2/𝒬^2+4x(1x)]^J}}`$ (100)
$`(\mathrm{pole}),`$
where

$`B_{E,0}^{}`$ $`=`$ $`2C_F{\displaystyle \frac{𝒬^2}{q^2}}\{16x(1x)\},`$
$`B_{E,1}^{}`$ $`=`$ $`2C_F\left\{12x(1x){\displaystyle \frac{𝒬^2}{q^2}}+56[x(1x)]^2{\displaystyle \frac{𝒬^2}{q^2}}1+2x(1x)\right\},`$
$`B_{E,2}^{}`$ $`=`$ $`64C_F[x(1x)]^2{\displaystyle \frac{𝒬^2}{q^2}}\{14x(1x)\}.`$ (101)
We expect that $`B_E(q^2)`$ is finite at $`q^2=0`$ when $`ϵ`$ is small and negative. It is, however, not evident from the form above that this is so. Nevertheless, explicit analytical integration shows that $`q^2B_E(q^2)`$ vanishes at $`q^2=0`$ for any $`ϵ`$. In order to get a form in which this vanishing is manifest, we simply subtract $`q^2B_E(q^2)`$ at $`q^2=0`$, which is zero, from $`q^2B_E(q^2)`$. Then $`B_E(q^2)`$ still has the form (95) but with new coefficients $`B_{E,J}`$:

$`B_{E,0}`$ $`=`$ $`0,`$
$`B_{E,1}`$ $`=`$ $`8C_Fx(1x),`$
$`B_{E,2}`$ $`=`$ $`16C_Fx(1x)\{1+4x(1x)\}.`$ (102)
The pole term vanishes, so we can simply set $`ϵ`$ to zero, obtaining

$$B_E(q^2)=\frac{\alpha _s}{8\pi ^2}𝑑\overline{q}^2𝑑x𝑑\varphi \frac{1}{\overline{q}^2q^2}\underset{J=0}{\overset{2}{}}\frac{B_{E,J}}{[\overline{q}^2/𝒬^2+4x(1x)]^J}$$
(103)
with the same coefficients (102).

### IX.3 $`W_q`$ and its supplementary terms

We can summarize our results so far by writing

$$F_q(q)=\frac{1}{2\pi }\frac{d\overline{q}^2}{\overline{q}^2q^2}𝑑x𝑑\varphi W_q(\overline{q}^2,x,\varphi ),$$
(104)
where $`W_q`$ has the form

$$W_q=\frac{\alpha _s}{4\pi }\left\{\text{/}qU_L(x,\overline{q}^2,q^2)+\frac{q^2}{𝒬^2}qn\text{/}nU_E(x,\overline{q}^2,q^2)\right\}.$$
(105)
with

$`U_L(x,\overline{q}^2,q^2)`$ $`=`$ $`{\displaystyle \underset{J=0}{\overset{2}{}}}{\displaystyle \frac{B_{L,J}}{[\overline{q}^2/𝒬^2+4x(1x)]^J}},`$
$`U_E(x,\overline{q}^2,q^2)`$ $`=`$ $`{\displaystyle \underset{J=0}{\overset{2}{}}}{\displaystyle \frac{B_{E,J}}{[\overline{q}^2/𝒬^2+4x(1x)]^J}}.`$ (106)
The coefficients $`B_{L,J}`$ and $`B_{E,J}`$ are given in Eqs. (99) and (102). Then the contribution to the full graph from the quark self-energy graph when the adjoining bare propagators are uncut is given in terms of $`W_q`$ by Eq. (84).

We can now set $`q^2`$ to 0 in $`W_q`$ to obtain $`𝒲_q`$. The second term in Eq. (105) vanishes and we have

$$𝒲_q[\mathrm{simple}]=\frac{\alpha _s}{4\pi }\text{/}qU_L(x,\overline{q}^2,0).$$
(107)
We see immediately that there will be a numerical problem for the cancellation of $`𝒲_q`$ with $`_q`$ at the small $`\overline{q}^2`$ endpoint of the integration (VII). In Eq. (87) for $`_q`$, there are terms with three Dirac matrix structures. The $`\text{/}n`$ term is not important because it multiplies a factor $`\mathrm{\Delta }`$ and $`\mathrm{\Delta }\overline{q}^2`$ for small $`\overline{q}^2`$. In the term proportional to $`\text{/}l_T`$, the singularity in $`_g^{\mu \nu }/\overline{q}^2`$ will vanish if we integrate over the angle $`\varphi `$ of $`l_T^\mu `$ before letting $`\overline{q}^2`$ become small. However, this is not suitable for a numerical integration. Finally, the coefficient of $`\text{/}q`$ in $`𝒲_q`$ contains only terms that are even under $`x(1x)`$, while the coefficient of $`\text{/}q`$ in $`_q`$ contains both even and odd terms. The small $`\overline{q}^2`$ singularity will be cancelled if we integrate over $`x`$ before letting $`\overline{q}^2`$ become small. But, again, this is not suitable for a numerical integration.

We can make the cancellation happen point by point in $`x`$ and $`\varphi `$ by adding two terms to $`𝒲_q`$, so that it becomes

$$𝒲_q=\frac{\alpha _s}{4\pi }\left\{\text{/}qU_L(x,\overline{q}^2,0)+\text{/}q\frac{V_L(x,\overline{q}^2)}{(1+\overline{q}^2/𝒬^2)}+\text{/}l_T\frac{N_t(x,\overline{q}^2)}{(1+\overline{q}^2/𝒬^2)}\right\}.$$
(108)
Here the function $`V_L(x,\overline{q}^2)`$ is related to the coefficient $`N_L(x,\overline{q}^2)`$ of $`\text{/}q`$ in the cut self-energy, Eq. (87), by

$$V_L(x,\overline{q}^2)=\frac{1}{2(1+\mathrm{\Delta })}\left[N_L(x,\overline{q}^2)N_L(1x,\overline{q}^2)\right].$$
(109)
(The factor $`1/(1+\mathrm{\Delta })`$ here makes the calculated expression for $`V_L`$ simpler.) The function $`N_t(x,\overline{q}^2)`$ is the coefficient of $`\text{/}l`$ in the cut self-energy, Eq. (87), and is given in Eq. (88).

The integrals of the extra terms vanish (because the integrations over $`x`$ and $`\varphi `$ respectively give zero), so we are adding zero to $`𝒲_q`$. However the cancellation in Eq. (VII) can now work point by point in $`\{\overline{q}^2,x,\varphi \}`$ space. (We check this for the $`\text{/}q`$ terms below.) We have inserted factors $`1/(1+\overline{q}^2/𝒬^2)`$ in the extra terms so as not to create problems at $`\overline{q}^2\mathrm{}`$ at the same time as we were alleviating problems at $`\overline{q}^20`$.

Let us summarize. The expression for $`𝒲_q`$ in Eq. (85) is now

$$𝒲_q=\frac{\alpha _s}{2\pi }\text{/}q𝒫_q(\overline{q}^2,x)\frac{\alpha _s}{4\pi }\text{/}l_T\frac{N_t(x,\overline{q}^2)}{(1+\overline{q}^2/𝒬^2)},$$
(110)
with

$`𝒫_q(\overline{q}^2,x)`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left\{U_L(x,\overline{q}^2,0)+{\displaystyle \frac{V_L(x,\overline{q}^2)}{(1+\overline{q}^2/𝒬^2)}}\right\}`$
$`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{J=0}{\overset{2}{}}}{\displaystyle \frac{B_{L,J}}{[\overline{q}^2/𝒬^2+4x(1x)]^J}}`$
$`+{\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{(1+\overline{q}^2/𝒬^2)}}{\displaystyle \underset{J=0}{\overset{2}{}}}{\displaystyle \frac{N_{L,J}}{[\overline{q}^2/𝒬^2+4x(1x)]^J}},`$
$`N_t(x,\overline{q}^2)`$ $`=`$ $`{\displaystyle \underset{J=0}{\overset{2}{}}}{\displaystyle \frac{N_{t,J}}{(\mathrm{\Delta }+2x)^J}}.`$ (111)
Notice that the denominators in the $`\text{/}l_T`$ term are different from those in the $`\text{/}q`$ terms. The coefficients $`B_{L,J}`$ from Eq. (99) are here evaluated at $`q^2=0`$:

$`B_{L,0}`$ $`=`$ $`C_F\left\{{\displaystyle \frac{e^3\mu ^2}{\overline{q}^2+e^3\mu ^2}}+12x(1x){\displaystyle \frac{e^{5/3}\mu ^2}{\overline{q}^2+e^{5/3}\mu ^2}}\right\},`$
$`B_{L,1}`$ $`=`$ $`8C_Fx(1x)\left\{514x(1x)\right\},`$
$`B_{L,2}`$ $`=`$ $`32C_Fx(1x)\left\{16x(1x)+8[x(1x)]^2\right\}.`$ (112)
The coefficients $`N_{L,J}`$ are computed from the coefficients for $`N_L(x,\overline{q}^2)`$ in Eq. (88) and are

$`N_{L,0}`$ $`=`$ $`C_F(2x1),`$
$`N_{L,1}`$ $`=`$ $`16C_F(2x1)x(1x),`$
$`N_{L,2}`$ $`=`$ $`32C_F(2x1)x(1x)[12x(1x)].`$ (113)
The coefficients $`N_{t,J}`$ for $`N_t(x,\overline{q}^2)`$ are given in Eq. (88):

$`N_{t,0}`$ $`=`$ $`2C_F,`$
$`N_{t,1}`$ $`=`$ $`4C_F(2x1),`$
$`N_{t,2}`$ $`=`$ $`16C_Fx(1x).`$ (114)
If we take the $`\overline{q}^20`$ limit of $`𝒲_q`$, the $`\text{/}l_T`$ term does not contribute since $`l_T0`$ as $`\overline{q}^20`$. We are left with

$$𝒲_q\frac{\alpha _s}{2\pi }\text{/}q𝒫_q(0,x).$$
(115)
Evaluating $`𝒫_q(\overline{q}^2,x)`$ at $`\overline{q}^2=0`$ we find

$$𝒫_q(0,x)=\stackrel{~}{P}_{g/q}(x),$$
(116)
where $`\stackrel{~}{P}_{g/q}(x)`$ is the parton evolution kernel given in Eq. (90). Thus $`𝒲_q`$ properly cancels $`_q`$ in Eq. (VII).

## X Renormalization of three point functions

In this section, we consider how to renormalize the divergent one loop virtual three point functions in Coulomb gauge using numerical integration.

### X.1 Quark-antiquark-boson vertices

In this subsection, we construct the renormalization counter term as an integral over the four dimensional space of loop momenta. We begin with the corresponding integrals over a $`42ϵ`$ dimensional space, since we want to match the renormalization to standard   MS  renormalization. We first study the quark-antiquark-gluon vertex. Then we extend the result to the quark-antiquark-photon vertex, which has a somewhat simpler structure.

There are two contributions to the quark-antiquark-gluon vertex $`\mathrm{\Gamma }_a^\mu (k_1,k_2)`$ at one loop. Each of them has the form

$$\mathrm{\Gamma }_a^\mu (k_1,k_2)=ig^2Ct_a\stackrel{~}{\mu }^{2ϵ}\frac{d^{42ϵ}l}{(2\pi )^{42ϵ}}\frac{N(l,k_1,k_2)_\nu ^\mu \gamma ^\nu }{[(k_2l)^2+iϵ][(k_1l)^2+iϵ][l^2+iϵ]},$$
(117)
where $`C`$ is the color factor for that graph and the scale factor $`\stackrel{~}{\mu }`$ is related to the   MS  scale factor $`\mu `$ by $`\mu ^2=4\pi \stackrel{~}{\mu }^2e^\gamma `$. When the loop momentum is large, the function $`N`$ has the form

$`N_\nu ^\mu `$ $`=`$ $`A_1(ϵ)l^2g_\nu ^\mu +A_2(ϵ){\displaystyle \frac{l^2}{\stackrel{~}{l}^2}}l^2g_\nu ^\mu `$ (118)
$`+A_3(ϵ)l^\mu l_\nu +A_4(ϵ){\displaystyle \frac{l^2}{\stackrel{~}{l}^2}}(l^\mu \stackrel{~}{l}_\nu +\stackrel{~}{l}^\mu l_\nu l^\mu l_\nu )`$
$`+A_5(ϵ){\displaystyle \frac{l^2}{\stackrel{~}{l}^2}}(l^\mu \stackrel{~}{l}_\nu \stackrel{~}{l}^\mu l_\nu )+A_6(ϵ)\left({\displaystyle \frac{l^2}{\stackrel{~}{l}^2}}\right)^2\stackrel{~}{l}^\mu (l_\nu \stackrel{~}{l}_\nu )`$
$`+𝒪(l).`$
Here the omitted terms are suppressed by one or more powers of $`k_1/l`$ or $`k_2/l`$.

We subtract a suitably chosen quantity $`\stackrel{~}{\mathrm{\Gamma }}_a^\mu `$ from $`\mathrm{\Gamma }_a^\mu (k_1,k_2)`$:

$$\stackrel{~}{\mathrm{\Gamma }}_a^\mu =ig^2Ct_a\stackrel{~}{\mu }^{2ϵ}\frac{d^{42ϵ}l}{(2\pi )^{42ϵ}}C(l)_\nu ^\mu \gamma ^\nu .$$
(119)
Using $`D(l^2)=l^2M^2+iϵ`$ and $`n=(1,0,0,0)`$ we define

$`C_\nu ^\mu `$ $`=`$ $`{\displaystyle \frac{A_1(ϵ)}{D(l^2)^2}}g_\nu ^\mu +{\displaystyle \frac{A_2(ϵ)}{D(l^2)D(\stackrel{~}{l}^2)}}g_\nu ^\mu `$ (120)
$`+{\displaystyle \frac{A_3(ϵ)}{D(l^2)^3}}l^\mu l_\nu +{\displaystyle \frac{A_4(ϵ)}{D(l^2)^2D(\stackrel{~}{l}^2)}}(l^\mu \stackrel{~}{l}_\nu +\stackrel{~}{l}^\mu l_\nu l^\mu l_\nu )`$
$`+{\displaystyle \frac{A_5(ϵ)}{D(l^2)^2D(\stackrel{~}{l}^2)}}(l^\mu \stackrel{~}{l}_\nu \stackrel{~}{l}^\mu l_\nu )+{\displaystyle \frac{A_6(ϵ)}{D(l^2)D(\stackrel{~}{l}^2)^2}}\stackrel{~}{l}^\mu (l_\nu \stackrel{~}{l}_\nu )`$
$`+{\displaystyle \frac{B_1(ϵ)M^2}{D(l^2)^3}}g_\nu ^\mu +{\displaystyle \frac{B_2(ϵ)M^2}{D(l^2)^3}}n^\mu n_\nu .`$
The functions $`B_1(ϵ)`$ and $`B_2(ϵ)`$ are defined in terms of the functions $`A_i(ϵ)`$ in such a way that $`\stackrel{~}{\mathrm{\Gamma }}`$ has a very simple dependence on $`ϵ`$. We will give the definition below.

Two properties of $`\stackrel{~}{\mathrm{\Gamma }}_a^\mu `$ are important. First, the leading $`(l)^4`$ behavior of the integrand of $`\stackrel{~}{\mathrm{\Gamma }}_a^\mu `$ matches that for $`\mathrm{\Gamma }_a^\mu `$ for large $`l^\mu `$. Second, it is easy to compute $`\stackrel{~}{\mathrm{\Gamma }}_a^\mu `$.

To begin the computation, we recognize that the integral of various terms in $`C_\nu ^\mu `$ will be proportional to combinations of $`g_\nu ^\mu `$ and $`n^\mu n_\nu `$. This allows us to replace the integrand by

$`C_\nu ^\mu `$ $``$ $`{\displaystyle \frac{A_1(ϵ)}{D(l^2)^2}}g_\nu ^\mu +{\displaystyle \frac{A_2(ϵ)}{D(l^2)D(\stackrel{~}{l}^2)}}g_\nu ^\mu +{\displaystyle \frac{A_3(ϵ)}{(42ϵ)}}g_\nu ^\mu \left({\displaystyle \frac{1}{D(l^2)^2}}+{\displaystyle \frac{M^2}{D(l^2)^3}}\right)`$ (121)
$`+A_4(ϵ){\displaystyle \frac{g_\nu ^\mu +(22ϵ)n^\mu n_\nu }{(32ϵ)}}\left({\displaystyle \frac{1}{D(l^2)^2}}+{\displaystyle \frac{M^2}{D(l^2)^2D(\stackrel{~}{l}^2)}}\right)`$
$`A_4(ϵ)n^\mu n_\nu \left({\displaystyle \frac{1}{D(l^2)D(\stackrel{~}{l}^2)}}+{\displaystyle \frac{M^2}{D(l^2)^2D(\stackrel{~}{l}^2)}}\right)`$
$`+{\displaystyle \frac{B_1(ϵ)M^2}{D(l^2)^3}}g_\nu ^\mu +{\displaystyle \frac{B_2(ϵ)M^2}{D(l^2)^3}}n^\mu n_\nu .`$
Note that the contributions from $`A_5`$ and $`A_6`$ vanish after integration because the tensors that multiply them in Eq. (120) vanish when contracted with $`g_\mu ^\nu `$ or $`n_\mu n^\nu `$.

Now we can perform the integration using

$$I_{J,K}i\stackrel{~}{\mu }^{2ϵ}\frac{d^{42ϵ}l}{(2\pi )^{42ϵ}}\frac{(M^2)^{J+K2}}{[l^2M^2+iϵ]^J[\stackrel{~}{l}^2M^2]^K}=\lambda (J,K,ϵ)\frac{\mathrm{\Gamma }(ϵ)}{16\pi ^2}\left(\frac{M^2}{4\pi \stackrel{~}{\mu }^2}\right)^ϵ,$$
(122)
where

$$\lambda (J,K,ϵ)=(1)^{J+K+1}\frac{\mathrm{\Gamma }(J1/2)\mathrm{\Gamma }(J+K2+ϵ)}{\mathrm{\Gamma }(J)\mathrm{\Gamma }(J+K1/2)\mathrm{\Gamma }(ϵ)}.$$
(123)
Specifically

$$\begin{array}{cc}\lambda (2,0,ϵ)=1,& \lambda (1,1,ϵ)=2,\\ \\ \lambda (3,0,ϵ)=ϵ/2,& \lambda (2,1,ϵ)=2ϵ/3.\end{array}$$
(124)
We find

$$\stackrel{~}{\mathrm{\Gamma }}_a^\mu =\frac{\alpha _s}{4\pi }Ct_aG_\nu ^\mu \gamma ^\nu \mathrm{\Gamma }(ϵ)\left(\frac{M^2}{4\pi \stackrel{~}{\mu }^2}\right)^ϵ$$
(125)
with (after taking cancellations among the $`ϵ`$ dependent factors into account)

$`G_\nu ^\mu `$ $`=`$ $`g_\nu ^\mu \left\{A_1(ϵ)+2A_2(ϵ)+{\displaystyle \frac{1}{4}}A_3(ϵ)+{\displaystyle \frac{1}{3}}A_4(ϵ){\displaystyle \frac{ϵ}{2}}B_1(ϵ)\right\}`$ (126)
$`+n^\mu n_\nu \left\{{\displaystyle \frac{4}{3}}A_4(ϵ)+{\displaystyle \frac{ϵ}{2}}B_2(ϵ)\right\}.`$
We set

$`B_1(ϵ)`$ $`=`$ $`{\displaystyle \frac{2}{ϵ}}\left\{[A_1(ϵ)A_1(0)]+2[A_2(ϵ)A_2(0)]+{\displaystyle \frac{1}{4}}[A_3(ϵ)A_3(0)]+{\displaystyle \frac{1}{3}}[A_4(ϵ)A_4(0)]\right\},`$
$`B_2(ϵ)`$ $`=`$ $`{\displaystyle \frac{8}{3ϵ}}[A_4(ϵ)A_4(0)].`$ (127)
Then $`G_\nu ^\mu `$ takes the $`ϵ`$-independent form

$$G_\nu ^\mu =g_\nu ^\mu \left\{A_1(0)+2A_2(0)+\frac{1}{4}A_3(0)+\frac{1}{3}A_4(0)\right\}+n^\mu n_\nu \frac{4}{3}A_4(0).$$
(128)
Expanding Eq. (125) about $`ϵ=0`$, we have

$$\stackrel{~}{\mathrm{\Gamma }}_a^\mu =\left[\stackrel{~}{\mathrm{\Gamma }}_a^\mu \right]_{\mathrm{pole}}+\left[\stackrel{~}{\mathrm{\Gamma }}_a^\mu \right]_\mathrm{R}+𝒪(ϵ),$$
(129)
with

$`\left[\stackrel{~}{\mathrm{\Gamma }}_a^\mu \right]_{\mathrm{pole}}`$ $`=`$ $`{\displaystyle \frac{\alpha _s}{4\pi }}Ct_aG_\nu ^\mu \gamma ^\nu {\displaystyle \frac{1}{ϵ}},`$
$`\left[\stackrel{~}{\mathrm{\Gamma }}_a^\mu \right]_\mathrm{R}`$ $`=`$ $`{\displaystyle \frac{\alpha _s}{4\pi }}Ct_aG_\nu ^\mu \gamma ^\nu \mathrm{ln}\left({\displaystyle \frac{\mu ^2}{M^2}}\right).`$ (130)
Define $`[\mathrm{\Gamma }_a^\mu ]_R`$ by

$$\mathrm{\Gamma }_a^\mu =\left[\mathrm{\Gamma }_a^\mu \right]_{\mathrm{pole}}+\left[\mathrm{\Gamma }_a^\mu \right]_\mathrm{R}+𝒪(ϵ).$$
(131)
We recognize that

$$\left[\mathrm{\Gamma }_a^\mu \right]_{\mathrm{pole}}=\left[\stackrel{~}{\mathrm{\Gamma }}_a^\mu \right]_{\mathrm{pole}}$$
(132)
because the $`l^4`$ behaviors of the two integrands match. Thus we can use our results for $`\stackrel{~}{\mathrm{\Gamma }}_a^\mu `$ to write

$$\left[\mathrm{\Gamma }_a^\mu \right]_\mathrm{R}=\left\{\mathrm{\Gamma }_a^\mu \stackrel{~}{\mathrm{\Gamma }}_a^\mu \right\}_{ϵ=0}+\left[\stackrel{~}{\mathrm{\Gamma }}_a^\mu \right]_\mathrm{R}.$$
(133)
The first term can be evaluated by numerical integration in 4 dimensions. The second term vanishes if we set

$$M^2=\mu ^2.$$
(134)
Now we need the coefficients $`A_J(ϵ)`$ in Eqs. (118) and (120). For the graph in which the gluon connects to the gluon line, one finds

$`C=C_A/2,`$
$`A_1(ϵ)=2,A_2(ϵ)=2,A_3(ϵ)=4(1ϵ),`$
$`A_4(ϵ)=2,A_5(ϵ)=0,A_6(ϵ)=2,`$
$`B_1(ϵ)=2,B_2(ϵ)=0.`$ (135)
Here the values for the $`B_i(ϵ)`$ have been calculated from the values for the $`A_i(ϵ)`$. For the graph in which the gluon connects to the quark line, one finds

$`C=C_FC_A/2=1/(2N_C),`$
$`A_1(ϵ)=2ϵ,A_2(ϵ)=1,A_3(ϵ)=4(1ϵ),`$
$`A_4(ϵ)=0,A_5(ϵ)=2,A_6(ϵ)=0,`$
$`B_1(ϵ)=2,B_2(ϵ)=0.`$ (136)
If we change $`C`$ from $`C_FC_A/2`$ to $`C_F`$, this same result holds for the quark-antiquark-photon vertex, with the appropriate change in the color structure from $`Ct_a`$ to $`C`$.

### X.2 Performing the energy integrations

We have seen how to renormalize the virtual three point functions in a fashion that works in four dimensions but is equivalent to   MS  renormalization. Here, we deal with the implementation of the loop integrals as numerical integrals. We need to write the counter terms for the three point function as integrals over the space components of the loop momentum:

$$\stackrel{~}{\mathrm{\Gamma }}_a^\mu =g^2t_a\frac{d^3\stackrel{}{l}}{(2\pi )^3}E(\stackrel{}{l})_\nu ^\mu \gamma ^\nu .$$
(137)
To do this, we perform the integrations over the loop energy analytically. Let us define the integrals

$$I_{N,J}=i\omega ^{2NJ1}\frac{dl^0}{2\pi }\frac{(l^0)^J}{D(l^2)^N},$$
(138)
where

$$\omega =\sqrt{\stackrel{}{l}^2+M^2}.$$
(139)
Then

$`E_\nu ^\mu `$ $`=`$ $`{\displaystyle \frac{CA_1(0)}{\omega ^3}}g_\nu ^\mu I_{2,0}{\displaystyle \frac{CA_2(0)}{\omega ^3}}g_\nu ^\mu I_{1,0}`$ (140)
$`+{\displaystyle \frac{CA_3(0)}{\omega ^5}}\left\{\stackrel{~}{l}^\mu \stackrel{~}{l}_\nu I_{3,0}+\omega ^2n^\mu n_\nu I_{3,2}\right\}{\displaystyle \frac{CA_4(0)}{\omega ^5}}\left\{\stackrel{~}{l}^\mu \stackrel{~}{l}_\nu I_{2,0}\omega ^2n^\mu n_\nu I_{2,2}\right\}`$
$`+{\displaystyle \frac{CB_1(0)M^2}{\omega ^5}}g_\nu ^\mu I_{3,0}+{\displaystyle \frac{CB_2(0)M^2}{\omega ^5}}n^\mu n_\nu I_{3,0}.`$
Here we have used the fact that the $`I_{N,J}`$ for odd $`J`$ vanish.

The integrands of $`I_{N,J}`$ have singular factors of the form

$$\frac{1}{D(l^2)^N}\frac{1}{(l^2M^2+iϵ)^N}=\frac{1}{(l^0\omega +iϵ)^N}\frac{1}{(l^0+\omega iϵ)^N}.$$
(141)
To perform the integrals, we use the prescription in Sec. VI.1. We close the integration contour in the lower half plane and then in the upper half plane and take the average of the results. We get $`I_{N,1}=0`$ and

$`I_{1,0}=1/2,`$ $`I_{1,2}=1/2,`$ (142)
$`I_{2,0}=1/4,`$ $`I_{2,2}=1/4,`$
$`I_{3,0}=3/16,`$ $`I_{3,2}=1/16.`$
Thus

$`E_\nu ^\mu `$ $`=`$ $`{\displaystyle \frac{CA_1(0)}{4\omega ^3}}g_\nu ^\mu {\displaystyle \frac{CA_2(0)}{2\omega ^3}}g_\nu ^\mu `$ (143)
$`+{\displaystyle \frac{CA_3(0)}{16\omega ^5}}\left\{3\stackrel{~}{l}^\mu \stackrel{~}{l}_\nu \omega ^2n^\mu n_\nu \right\}+{\displaystyle \frac{CA_4(0)}{4\omega ^5}}\left\{\stackrel{~}{l}^\mu \stackrel{~}{l}_\nu +\omega ^2n^\mu n_\nu \right\}`$
$`+{\displaystyle \frac{3CB_1(0)M^2}{16\omega ^5}}g_\nu ^\mu +{\displaystyle \frac{3CB_2(0)M^2}{16\omega ^5}}n^\mu n_\nu .`$
We can use our explicit results for the coefficients $`A_J`$ and $`B_J`$ to obtain the net result for the counter term for the quark-antiquark-gluon one loop graphs (summed over the two graphs), expressed as an integral over $`\stackrel{}{l}`$. The counter term is given by Eq. (137) with

$$E_\nu ^\mu =\frac{C_F}{8\omega ^5}(4\omega ^2+3M^2)g_\nu ^\mu +\frac{3C_F4C_A}{4\omega ^5}\stackrel{~}{l}^\mu \stackrel{~}{l}_\nu \frac{C_F}{4\omega ^3}n^\mu n_\nu .$$
(144)
According to Eq. (133), we are to subtract $`\stackrel{~}{\mathrm{\Gamma }}_a^\mu `$ from the integral for $`\mathrm{\Gamma }_a^\mu `$ and set $`M`$ = $`\mu `$.

For the quark-antiquark-photon graph the counter term is given by Eq. (137) without the factor of $`t_a`$ and with

$$E_\nu ^\mu =\frac{C_F}{8\omega ^5}(4\omega ^2+3M^2)g_\nu ^\mu +\frac{3C_F}{4\omega ^5}\stackrel{~}{l}^\mu \stackrel{~}{l}_\nu \frac{C_F}{4\omega ^3}n^\mu n_\nu .$$
(145)
## XI Results

In this section, we look at some numerical results. As our example, we will consider one of the standard event shape variables, the thrust $`t`$. We examine the thrust distribution normalized to the total cross section $`\sigma (1+\alpha _s/\pi )\times \sigma _0`$ for $`e^++e^{}\mathrm{ℎ𝑎𝑑𝑟𝑜𝑛𝑠}`$,

$$(t)=\frac{1}{\sigma }\frac{d\sigma }{dt}.$$
(146)
For $`t<1`$, $`d\sigma /dt`$ has a contribution, $`d\sigma ^{(1)}/dt`$, of order $`\alpha _s^1`$ and a contribution, $`d\sigma ^{(2)}/dt`$, of order $`\alpha _s^2`$. Both contributions are included in the next-to-leading order results for $`(t)`$. We also isolate the second order term, $`d\sigma ^{(2)}/dt`$ and study the second order contributions to the moments of the thrust distribution,

$$I_n=\frac{1}{\sigma _0(\alpha _s/\pi )^2}_0^1𝑑t(1t)^n\frac{d\sigma ^{(2)}}{dt}.$$
(147)
The first question is whether the gauge choice makes any difference. There are ten topologies of Feynman diagrams that contribute order $`\alpha _s^2`$ terms to $`e^+e^{}\mathrm{ℎ𝑎𝑑𝑟𝑜𝑛𝑠}`$. These are shown in Fig. 6. For each topology, we calculate the corresponding contribution to the second moment of the thrust distribution, $`I_2`$. The results are shown in Table 2. We see that graph by graph, the results are completely different in Feynman and Coulomb gauges. However, the total $`I_2`$ summed over graphs is independent of the gauge.

Next, we test whether the Coulomb gauge calculation is working properly by checking whether $`I_n`$ calculated in Coulomb gauge matches $`I_n`$ calculated in Feynman gauge for several choices of $`n`$. (The results for Feynman gauge were checked against the program of Kunszt and Nason KN  in Ref. beowulfPRD .) The results are presented in Table 3. We see that the results are properly gauge invariant within the errors of the program.

Having seen that the program appears to be working properly, we exhibit a graph of the next-to-leading order thrust distribution $`(t)`$ versus t in Fig. 7. We also show the same distribution calculated at leading order and data from the Opal Collaboration thedata . The theoretical results are rather sensitive to the choice of the $`\overline{\mathrm{MS}}`$ renormalization scale $`\mu `$. We have chosen $`\mu `$ to be half of a typical jet energy in a three jet event, $`\sqrt{s}/3`$. That is, $`\mu =\sqrt{s}/6`$. The agreement between theory and data is not perfect, but this is to be expected in a strictly perturbative expansion that includes only the first two terms and no correction for effects beyond perturbation theory such as hadronization effects.

As mentioned in the introduction, one may wish to go beyond pure next-to-leading order calculations by incorporating, in an approximate way, some effects at all orders in $`\alpha _s`$. For instance, one may want to include renormalon effects by letting $`\alpha _s`$ run as a function of loop momenta inside graphs. Alternatively, one may also want to simulate realistic final states by adding parton showers to the next-to-leading order calculation. For such applications, one must approximate, and the presence of unphysical degrees of freedom propagating over long distances makes approximation difficult. A straightforward remedy is to do the calculation in a physical gauge, such as Coulomb gauge. In this paper we have seen how this goal can be accomplished.

## Appendix A Results for Feynman gauge

It is useful to have at hand the formulas for Feynman gauge that are analogous to what we have found in Coulomb gauge. We present the needed formulas in this appendix.

For the gluon propagator with a self-energy insertion, there is a problem that we treat as described in beowulfPRL ; beowulfPRD . The problem is most easily seen with the virtual self-energy graph. We know that $`\mathrm{\Pi }^{\alpha \beta }(q^2g^{\alpha \beta }q^\alpha q^\beta )`$. The first term is fine since it contains a factor $`q^2`$ that cancels the $`1/q^2`$ in the adjoining propagator. The next term does not have this good property, but its contribution will cancel in a sum over graphs because it is proportional to $`q^\alpha `$ and $`q^\beta `$. In order to make this cancellation happen in a single graph, we replace

$$(g_\alpha ^\mu )(g_\beta ^\nu )\mathrm{\Pi }^{\alpha \beta }$$
(148)
by

$$(g_\alpha ^\mu +q^\mu \stackrel{~}{q}_\alpha /\stackrel{~}{q}^2)(g_\beta ^\nu +q^\nu \stackrel{~}{q}_\beta /\stackrel{~}{q}^2)\mathrm{\Pi }^{\alpha \beta }.$$
(149)
This does not change the answer after summing over graphs because the added terms are proportional to $`q^\mu `$ or $`q^\nu `$. But now the $`q^\alpha q^\beta `$ term in $`\mathrm{\Pi }^{\alpha \beta }`$ gives zero in each graph because $`(g_\alpha ^\mu +q^\mu \stackrel{~}{q}_\alpha /\stackrel{~}{q}^2)q^\alpha =0`$. We make this replacement for both the virtual and real versions of $`\mathrm{\Pi }^{\alpha \beta }`$.

For the real gluon self-energy graph, Eq. (40) becomes

$`N_{TT}`$ $`=`$ $`2C_A\left\{1+x(1x)\right\}+N_F\left\{12x(1x)\right\},`$
$`N_{tt}`$ $`=`$ $`4C_Ax(1x)4N_Fx(1x),`$
$`N_{EE}`$ $`=`$ $`C_A\left\{1+4x(1x)\right\}+4N_Fx(1x),`$
$`N_{Et}`$ $`=`$ $`2C_A(2x1)+2N_F(2x1).`$ (150)
For the virtual gluon self-energy graph, Eq. (53) becomes

$`A_{T,0}^{\prime \prime }`$ $`=`$ $`C_A{\displaystyle \frac{1}{q^2}}\left\{{\displaystyle \frac{1}{2}}\overline{q}^2{\displaystyle \frac{5}{2}}q^2+2\overline{q}^2x(1x)\right\}+N_F{\displaystyle \frac{1}{q^2}}\left\{\overline{q}^2{\displaystyle \frac{x(1x)}{(1ϵ)}}2\overline{q}^2\right\},`$
$`A_{T,1}^{\prime \prime }`$ $`=`$ $`0,`$
$`A_{T,2}^{\prime \prime }`$ $`=`$ $`0.`$ (151)
After subtracting $`q^2A_T(q^2)`$ at $`q^2=0`$ from $`q^2A_T(q^2)`$, we are left with the revised version of Eq. (54),

$`A_{T,0}^{}`$ $`=`$ $`2C_A\left\{1+x(1x)\right\}+N_F\left\{1{\displaystyle \frac{2x(1x)}{(1ϵ)}}\right\},`$
$`A_{T,1}^{}`$ $`=`$ $`0,`$
$`A_{T,2}^{}`$ $`=`$ $`0.`$ (152)
After renormalization, we have the revised version of Eq. (60),

$`A_{T,0}`$ $`=`$ $`(2C_AN_F){\displaystyle \frac{q^2+e^2\mu ^2}{\overline{q}^2+e^2\mu ^2}}+2C_Ax(1x){\displaystyle \frac{q^2+e^{5/3}\mu ^2}{\overline{q}^2+e^{5/3}\mu ^2}}`$
$`2N_Fx(1x){\displaystyle \frac{q^2+e^{8/3}\mu ^2}{\overline{q}^2+e^{8/3}\mu ^2}},`$
$`A_{T,1}`$ $`=`$ $`0,`$
$`A_{T,2}`$ $`=`$ $`0.`$ (153)
In Feynman gauge, the second term in the decomposition Eq. (71) is not there, so that

$$F_g^{\mu \nu }(q)=D(q)^{\mu \nu }A_T(q^2).$$
(154)
That is, $`A_E(q^2)=A_T(q^2)`$ in Eq. (71).

For the quark propagator with a self-energy insertion, we simply change the gauge in the previous Coulomb gauge calculation. Then the coefficients for a quark propagator with a cut self-energy diagram are given by a revised form of Eq. (88),

$`N_L`$ $`=`$ $`C_F\left\{4x(1x)+(2x1)(2x+\mathrm{\Delta })\right\},`$
$`N_E`$ $`=`$ $`2C_F(1x),`$
$`N_t`$ $`=`$ $`2C_F.`$ (155)
For the virtual quark self-energy, we begin with the revised form of Eq. (96),

$`B_{L,0}^{\prime \prime }`$ $`=`$ $`C_F(1ϵ),`$
$`B_{L,1}^{\prime \prime }`$ $`=`$ $`0,`$
$`B_{L,2}^{\prime \prime }`$ $`=`$ $`0.`$ (156)
Then $`B_{L,J}^{}=B_{L,J}^{\prime \prime }`$. After renormalization, we have the revised form of Eq. (99),

$`B_{L,0}`$ $`=`$ $`C_F{\displaystyle \frac{q^2+e^1\mu ^2}{\overline{q}^2+e^1\mu ^2}},`$
$`B_{L,1}`$ $`=`$ $`0,`$
$`B_{L,2}`$ $`=`$ $`0.`$ (157)
In Feynman gauge, the second term in the decomposition (93) is not there, so that

$$F_q(q)=\text{/}qB_L(q^2).$$
(158)
That is, $`B_E(q^2)=0`$ in Eq. (93).

The coefficients needed in the renormalization of the quark-antiquark-vector boson three point functions change when we go from Coulomb gauge to Feynman gauge. Specifically, in place of Eq. (135) we have

$`C=C_A/2,`$
$`A_1(ϵ)=2,A_2(ϵ)=0,A_3(ϵ)=4(1ϵ),`$
$`A_4(ϵ)=0,A_5(ϵ)=0,A_6(ϵ)=0,`$
$`B_1(ϵ)=2,B_2(ϵ)=0.`$ (159)
In place of Eq. (136) we have

$`C=C_FC_A/2=1/(2N_C),`$
$`A_1(ϵ)=2(1ϵ),A_2(ϵ)=0,A_3(ϵ)=4(1ϵ),`$
$`A_4(ϵ)=0,A_5(ϵ)=0,A_6(ϵ)=0,`$
$`B_1(ϵ)=2,B_2(ϵ)=0.`$ (160)
###### Acknowledgements.

This work was supported in part by the U.S. Department of Energy, by the European Union under contract HPRN-CT-2000-00149, and by the British Particle Physics and Astronomy Research Council.