# QCD CALCULATIONS BY NUMERICAL INTEGRATION

## 1 Introduction

This talk concerns a method for performing perturbative calculations in Quantum Chromodynamics and other quantum field theories.  The method may prove useful for cross sections in which one measures something about the hadronic final state, so that one cannot use the powerful but very special methods that apply to completely inclusive quantities like the total cross section for $`e^+e^{}\mathrm{ℎ𝑎𝑑𝑟𝑜𝑛𝑠}`$. Examples of such cross sections include jet cross sections in hadron-hadron and lepton-hadron scattering and in $`e^+e^{}\mathrm{ℎ𝑎𝑑𝑟𝑜𝑛𝑠}`$.

Many calculations of this kind have been carried out at next-to-leading order in perturbation theory, resulting in an impressive confrontation between theory and experiment. These calculations have been based on a method introduced by Ellis, Ross, and Terrano  in the context of $`e^+e^{}\mathrm{ℎ𝑎𝑑𝑟𝑜𝑛𝑠}`$. Stated in the simplest terms, one performs some integrations over momenta $`\stackrel{}{p}_i`$ analytically, others numerically. I shall argue that it is possible instead to do all of these integrations numerically. This has some advantages, principally in the flexibility that it allows.

In this talk, I address specifically the calculation of three-jet-like infrared safe observables in $`e^+e^{}\mathrm{ℎ𝑎𝑑𝑟𝑜𝑛𝑠}`$ at next-to-leading order, that is order $`\alpha _s^2`$. Examples of such observables include the thrust distribution, the fraction of events that have three jets, and the energy-energy correlation function. There are, of course, several existing computer programs that can calculate any such quantity, the earliest of these being that of Kunszt and Nason,  which was developed in 1989. What is new in this talk is not the result achieved, but the method of calculation.

## 2 The quantity to be calculated.

The order $`\alpha _s^2`$ contribution to the observable being calculated has the form

$`\sigma ^{[2]}`$ $`=`$ $`{\displaystyle \frac{1}{3!}}{\displaystyle 𝑑\stackrel{}{k}_1𝑑\stackrel{}{k}_2𝑑\stackrel{}{k}_3\frac{d\sigma _3^{[2]}}{d\stackrel{}{k}_1d\stackrel{}{k}_2d\stackrel{}{k}_3}𝒮_3(\stackrel{}{k}_1,\stackrel{}{k}_2,\stackrel{}{k}_3)}`$   
$`+{\displaystyle \frac{1}{4!}}{\displaystyle 𝑑\stackrel{}{k}_1𝑑\stackrel{}{k}_2𝑑\stackrel{}{k}_3𝑑\stackrel{}{k}_4\frac{d\sigma _4^{[2]}}{d\stackrel{}{k}_1d\stackrel{}{k}_2d\stackrel{}{k}_3d\stackrel{}{k}_4}𝒮_4(\stackrel{}{k}_1,\stackrel{}{k}_2,\stackrel{}{k}_3,\stackrel{}{k}_4)}.`$   
There are, of course, infrared divergences associated with Eq. (2); we understand that an infrared cutoff has been supplied. In Eq. (2), the $`d\sigma _n^{[2]}`$ are the order $`\alpha _s^2`$ contributions to the parton level cross section, calculated with zero quark masses. Each contains momentum and energy conserving delta functions. The $`d\sigma _n^{[2]}`$ include ultraviolet renormalization in the   MS  scheme. The functions $`𝒮`$ describe the measurable quantity to be calculated. Since we calculate a “three-jet” quantity, $`𝒮_2=0`$. The measurement, as specified  by the functions $`𝒮_n`$, is to be infrared safe: the $`𝒮_n`$ are smooth functions of the parton momenta and

$$𝒮_{n+1}(\stackrel{}{k}_1,\mathrm{},\lambda \stackrel{}{k}_n,(1\lambda )\stackrel{}{k}_n)=𝒮_n(\stackrel{}{k}_1,\mathrm{},\stackrel{}{k}_n)$$
(2)   
for $`0\lambda <1`$. That is, collinear splittings and soft particles do not affect the measurement.

It is convenient to calculate a quantity that is dimensionless. Let us, therefore, define the functions $`𝒮_n`$ so that they are dimensionless and eliminate the remaining dimensionality in the problem by dividing by $`\sigma _0`$, the total $`e^+e^{}`$ cross section at the Born level. We also remove the factor of $`(\alpha _s/\pi )^2`$. Thus, we calculate

$$=\frac{\sigma ^{[2]}}{\sigma _0(\alpha _s/\pi )^2}.$$
(3)   
## 3 Setting up the calculation

We note that $``$ is a function of the c.m. energy $`\sqrt{s}`$ and the $`\overline{\mathrm{MS}}`$ renormalization scale $`\mu `$. We will choose $`\mu `$ to be proportional to $`\sqrt{s}`$: $`\mu =A_{UV}\sqrt{s}`$. Then $``$ depends on $`A`$. But, because it is dimensionless, it is independent of $`\sqrt{s}`$. This allows us to write

$$=_0^{\mathrm{}}𝑑\sqrt{s}h(\sqrt{s})(A_{UV},\sqrt{s}),$$
(4)   
where $`h`$ is any function with

$$_0^{\mathrm{}}𝑑\sqrt{s}h(\sqrt{s})=1.$$
(5)   
The integration over $`\sqrt{s}`$ eliminates the energy conserving delta function in $``$. This is important in allowing cancellations of infrared singularities between contributions with 3- and 4-particle intermediate states to take place. The physical meaning is that, by smearing in the energy $`\sqrt{s}`$, we force the time variables in the two current operators that create the hadronic state to be within $`1/\sqrt{s}`$ of each other. Thus we have a truly short distance problem.

## 4 The Ellis, Ross, and Terrano method

I now describe how one would calculate $``$ using the Ellis-Ross-Terrano method. Each partonic cross section in Eq. (2) can be expressed as an amplitude times a complex conjugate amplitude, as illustrated in Fig. 1. One must calculate the amplitudes in $`42ϵ`$ dimensions. (In the case of the process $`e^+e^{}\mathrm{ℎ𝑎𝑑𝑟𝑜𝑛𝑠}`$, this calculation was performed by Ellis, Ross, and Terrano. ) For tree diagrams, the calculation is straightforward. For loop diagrams, this involves an integration, which is performed analytically. The integrals are divergent in four dimensions, so one obtains divergent terms proportional to $`1/ϵ^2`$ and $`1/ϵ`$ in addition to terms that are finite as $`ϵ0`$. Having the amplitudes and complex conjugate amplitudes, one must now multiply by the functions $`𝒮_n`$ and integrate over the final state parton momenta. These integrations are too complicated to perform analytically, so one must use numerical methods. Unfortunately, the integrals are divergent at $`ϵ=0`$. Thus one must split the integrals into two parts. One part can be divergent at $`ϵ=0`$ but must be simple enough to calculate analytically. The other part can be complicated, but must be convergent at $`ϵ=0`$. One calculates the simple, divergent part and cancels the $`1/ϵ^2`$ and $`1/ϵ`$ pole terms against the pole terms coming from the virtual loop diagrams. This leaves the complicated, convergent integration to be performed numerically.

This method is a little bit cumbersome, but it works and has been enormously successful. However it has proven to be difficult to apply the method in the case of two virtual loops. Even with one virtual loop, the method is not very flexible. In modern implementions of the Ellis- Ross-Terrano method, one can easily choose the functions $`𝒮_n`$, but any other modification of the integrand requires one to recalculate the amplitudes, and the modification must be simple enough that one can calculate the amplitudes in closed form.

## 5 Calculation by numerical integration

Let us, therefore, inquire whether there is any other way that one might perform such a calculation. We note that the quantity $``$ can be expressed in terms of cut Feynman diagrams, as in Fig. 2. Each diagram is a three loop diagram, so we have integrations over loop momenta $`\mathrm{}_1^\mu `$, $`\mathrm{}_2^\mu `$ and $`\mathrm{}_3^\mu `$.

### 5.1 Performing the energy integrations

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 (k_J^2)`$. These delta functions eliminate the three energy integrals over $`\mathrm{}_1^0`$, $`\mathrm{}_2^0`$, and $`\mathrm{}_3^0`$ as well as the integral (3) over $`\sqrt{s}`$. For the graphs in which three parton lines cross the cut, we can eliminate the integration over $`\sqrt{s}`$ and two of the $`\mathrm{}_J^0`$ integrals. One integral over the energy $`E`$ in the virtual loop remains. The integrand contains a product of factors

$$\frac{i}{EQ_J^0\omega _J+iϵ}\frac{1}{EQ_J^0+\omega _Jiϵ},$$
(6)   
where $`EQ_J^0`$ is the energy carried by the $`J`$th propagator around the loop and $`\omega _J`$ is the absolute value of the momentum carried on that propagator. We perform the integration by closing the integration contour in the lower half $`E`$ plane. This leads to $`n`$ terms for a virtual $`n`$ point subgraph. In the $`J`$th term, the propagator energy $`EQ_J^0`$ is set equal to the corresponding $`\omega _J`$ and there is a factor $`1/(2\omega _J)`$. Note that the entire process of performing the energy integrals amounts to some simple algebraic substitutions.

### 5.2 Cancellations of pinch singularities

Let us denote the contribution to $``$ from the cut $`C`$ of graph $`G`$ by $`(G,C)`$. This contribution has the form

$$(G,C)=d^3\stackrel{}{\mathrm{}}_1d^3\stackrel{}{\mathrm{}}_2d^3\stackrel{}{\mathrm{}}_3g(G,C;\mathrm{}),$$
(7)   
where $`\mathrm{}`$ denotes the set of loop momenta $`\{\stackrel{}{\mathrm{}}_1,\stackrel{}{\mathrm{}}_2,\stackrel{}{\mathrm{}}_3\}`$. The functions $`g`$ have some singularities, called pinch singularities, that cannot be avoided by deforming the integration contour and some non-pinch singularities that can be avoided by a contour deformation.

Let us consider the pinch singularities, following the analysis of Sterman.  The pinch singularities occur when one parton branches into two partons with collinear momenta or when one parton momentum goes to zero. These singularities can lead to logarithmic divergences in the corresponding integral. (There is a complication associated with the $`q^\mu q^\nu `$ terms in gluon self-energy subgraphs but I ignore this complication in this talk.) If we were to calculate the total cross section by using measurement functions $`𝒮_n=1`$ in $`g`$, then the singularities would cancel between the functions $`g`$ associated with the various cuts $`C`$ of the same graph $`G`$. The underlying reason is unitarity. Now in our case of measured shape variables, the values of $`𝒮_n`$ corresponding to different cuts $`C`$ are different. This would ruin the cancellation, except that just at the collinear or soft points the functions $`𝒮_n`$ match. Thus the singularities present in the individual $`g(G,C;\mathrm{})`$ cancel in the sum, $`_Cg(G,C;\mathrm{})`$. To be precise, the cancellation reduces the strength of the singularity from a strength sufficient to give a logarithmically divergent integral to one that gives a convergent integral.

### 5.3 Avoiding the non-pinch singularities

The function $`g`$ also has singularities that can be avoided by deforming the integration contours, the non-pinch singularities. These singularities do not cancel, but the Feynman rules provide an $`iϵ`$ prescription that tells us that we should deform the $`\stackrel{}{\mathrm{}}`$ integration contour into the complex $`\stackrel{}{\mathrm{}}`$ plane so as to avoid the singularity. Here deforming the contour means replacing $`\stackrel{}{\mathrm{}}`$ by a complex vector $`\stackrel{}{\mathrm{}}+i\stackrel{}{\kappa }`$. Then one simply chooses the imaginary part, $`\stackrel{}{\kappa }`$, of the loop momentum as a function of the real part, $`\stackrel{}{\mathrm{}}`$, and supplies the appropriate jacobian $`𝒥`$. Then

$$(G,C)=d^3\stackrel{}{\mathrm{}}_1d^3\stackrel{}{\mathrm{}}_2d^3\stackrel{}{\mathrm{}}_3𝒥(G,C;\mathrm{})g(G,C;\mathrm{}+i\kappa (G,C;\mathrm{})).$$
(8)   
### 5.4 The integration

With the contours appropriately chosen, the integral

$$(G)=d^3\stackrel{}{\mathrm{}}_1d^3\stackrel{}{\mathrm{}}_2d^3\stackrel{}{\mathrm{}}_3\underset{C}{}𝒥(G,C;\mathrm{})g(G,C;\mathrm{}+i\kappa (G,C;\mathrm{}))$$
(9)   
is finite. One can simply compute it by Monte Carlo integration after removing from the integration tiny regions near the collinear and soft singular points, where roundoff errors spoil the cancellation of the individual contributions.

Note the significance of putting the summation over cuts inside the integral. When we sum over cuts for a given point in the space of loop momenta, the soft and collinear divergences cancel because the cancellation is built into the Feynman rules. If we were to put the sum over cuts outside the integration, as in the Ellis-Ross-Terrano method, then the individual integrals would be divergent. The calculation would thereby be rendered more difficult.

## 6 Additional points

Eq. (9) represents the main point of this talk. There are two additional points of principle that I discuss below: renormalization and the special treatment required for self-energy subdiagrams.

### 6.1 Renormalization

Some of the virtual loop subgraphs are ultraviolet divergent. Normally, one would renormalize these subgraphs by performing loop integrations in $`42ϵ`$ space-time dimensions and subtracting the resulting pole term, $`c/ϵ`$. Clearly that is not appropriate in a numerical integration. However, one can subtract instead an integral in 4 space-time dimensions such that, in the region of large loop momenta, the integrand of the subtraction term matches the integrand of the subdiagram in question. The integrand of the subtraction term can depend on a mass parameter $`\mu `$ in such a way that the subtraction term has no infrared singularities. (For example, in the simplest case one can have denominators of the form $`1/[\stackrel{}{\mathrm{}}^2+\mu ^2]`$.) Then, one can easily arrange the definition so that this subtraction has exactly the same effect as the conventional   MS  subtraction with scale parameter $`\mu `$.

### 6.2 Self-energy subgraphs

Virtual self-energy subgraphs, as in the right hand diagram in Fig. 3, require a special treatment. Consider a quark self-energy subdiagram $`i\mathrm{\Sigma }`$ with one adjoining virtual propagator and one adjoining cut propagator. This combination really represents a field strength renormalization for the quark field, and is interpreted as

$$\frac{1}{2}\left[\frac{\text{/}q}{q^2}\mathrm{\Sigma }(q)\text{/}q\right]_{q^2=0}\mathrm{\hspace{0.17em}2}\pi \delta (q^2).$$
(10)   
In order to take the $`q^20`$ limit here while at the same time maintaining the cancellation of collinear divergences, we write

$$\frac{\text{/}q\mathrm{\Sigma }(q)\text{/}q}{q^2}=\frac{g^2}{(2\pi )^3}C_Fd^3\stackrel{}{\mathrm{}}\frac{q^0|\stackrel{}{k}_+|\gamma ^0(|\stackrel{}{k}_+|+|\stackrel{}{k}_{}|)\stackrel{}{k}_+\stackrel{}{\gamma }}{|\stackrel{}{k}_+||\stackrel{}{k}_{}|\left[(|\stackrel{}{k}_+|+|\stackrel{}{k}_{}|)^2(q^0)^2\right]},$$
(11)   
where $`\stackrel{}{k}_\pm =\frac{1}{2}\stackrel{}{q}\pm \stackrel{}{\mathrm{}}`$. This expression is obtained by writing the left-hand side as a dispersive integral with the cut self-energy graph appearing as the discontinuity. When the virtual self-energy is written in this form, the $`q^20`$ limit is smooth and, in addition, the cancellation between the two cut diagrams in Fig. 3 in the collinear limit $`\stackrel{}{\mathrm{}}\stackrel{}{q}`$ is manifest. It should be noted that the integral in Eq. (11) is ultraviolet divergent and requires renormalization, which can be performed with an ad hoc subtraction as described above.

One may expect that the representation of virtual propagator corrections in terms of the cut propagator will prove to be convenient in future modifications of the method described here. One may want to make modifications to the gluon propagator, in particular, in order to implement a running strong coupling and to insert models for the long distance propagation of gluons. In addition, one may want to modify propagators so that partons can enter the final state with $`q^2>0`$ so as to mesh an order $`\alpha _s^2`$ perturbative calculation with a parton shower Monte Carlo program. In either case, the modified virtual gluon propagator should be related by a dispersive representation to the modified final state created by an outgoing gluon.

## 7 An example at order $`\alpha _s`$

I have argued that a completely numerical calculation of the $`e^+e^{}`$ observables considered here is possible in principle. The question now arises whether the principle can be put into operation in an actual computer program.

Consider first a simple example at one lower order in $`\alpha _s`$: the total cross section for $`e^+e^{}hadrons`$ at order $`\alpha _s`$. The cut diagrams are shown in Fig. 4.

Let us define $`A_1(\mu )`$ to be the sum of the top graphs in Fig. 4 and $`A_2(\mu )`$ to be the sum of the bottom graphs, so that

$$\sigma =\sigma _0\left\{1+\left[A_1(\mu )+A_2(\mu )\right]\frac{\alpha _s}{\pi }\right\}.$$
(12)   
Then a simple analytical calculation shows that

$$A_1(\mu )=\frac{5}{3}\frac{4}{3}\mathrm{ln}(\mu /\sqrt{s}),A_2(\mu )=\frac{8}{3}+\frac{4}{3}\mathrm{ln}(\mu /\sqrt{s}).$$
(13)   
I have constructed a demonstration computer program to calculate $`A_1(\mu )`$ and $`A_2(\mu )`$ numerically along the lines described above. This example has UV renormalization, propagators to be treated dispersively, singularities to be avoided by deforming the contour, soft singularities that must cancel, and collinear singularities that must cancel. Thus it is not trivial from a numerical point of view. The results are shown in Fig. 5. Evidently, the method works.

## 8 An example at order $`\alpha _s^2`$

In order to further test the methods described here, I have constructed a demonstration computer program  to calculate three-jet-like observables for $`e^+e^{}`$ annihilation at order $`\alpha _s^2`$. I have used the program to calculate $`d\sigma ^{[2]}/dT`$, the order $`\alpha _s^2`$ contribution to the thrust distribution. More precisely, I have calculated the ratio $`R(T)`$ of $`d\sigma ^{[2]}/dT`$ to $`[d\sigma ^{[2]}/dT]_{\mathrm{KN}}`$, where where $`[d\sigma ^{[2]}/dT]_{\mathrm{KN}}`$ is a fit to the tabulated results for $`d\sigma ^{[2]}/dT`$ as given by Kunszt and Nason.  In the range $`0.71<T<0.95`$, the function $`[d\sigma ^{[2]}/dT]_{\mathrm{KN}}`$ varies by about a factor of 30. The ratio $`R(T)`$ should be 1. The results are reported in Fig. 6. Again, it appears that the numerical method of calculation works in a practical sense.

## 9 Outlook

I am working to test the code of the demonstration program. Indeed, between the time of the RADCOR98 conference and the time of this writing, I have found and corrected a small error in the program logic which resulted in an error in the calculation of about 1%. With this bug corrected, the results produced by the current version of the program are better than those shown in Fig. 6. After testing the code, I will have to document the code and publish a detailed description of the algorithms used.

With this code in hand, or with improved code from other authors, I anticipate attacks on more difficult problems than the problem discussed here. It remains to be seen for what problems the completely numerical method will prove to be more powerful than the analytical/numerical method that has served us so well up to now.

## Acknowledgment

This work was supported by the U. S. Department of Energy.

## References