# Techniques for QCD calculations by numerical integration

## I Introduction

This paper concerns a method, which was introduced in , for performing perturbative calculations in quantum chromodynamics (QCD) and other quantum field theories. The method is intended for calculations of quantities in which one measures something about the hadronic final state produced in a collision and in which the observable is infrared safe – that is, insensitive to long-distance effects. Examples include jet cross sections in hadron-hadron and lepton-hadron scattering and in $`e^+e^{}\mathrm{ℎ𝑎𝑑𝑟𝑜𝑛𝑠}`$. There have been many calculations of this kind carried out at next-to-leading order in perturbation theory. These calculations are based on a method introduced by Ellis, Ross, and Terrano  in the context of $`e^+e^{}\mathrm{ℎ𝑎𝑑𝑟𝑜𝑛𝑠}`$. Stated in the simplest terms, the Ellis-Ross-Terrano method is to do some integrations over momenta $`\stackrel{}{\mathrm{}}_i`$ analytically, others numerically. In the method discussed here, one does all of these integrations numerically. Evidently, if one performs all of the integrations numerically, one gains flexibility to quite easily modify the integrand. There may be other advantages, as well as some disadvantages, to the numerical integration method compared to the numerical/analytical method.

In this paper, I address only the process $`e^+e^{}\mathrm{ℎ𝑎𝑑𝑟𝑜𝑛𝑠}`$. I discuss three-jet-like infrared safe observables at next-to-leading order, that is order $`\alpha _s^2`$. Examples of such observables include the thrust distribution and the fraction of events that have three jets.

The main techniques of the numerical integration method for $`e^+e^{}\mathrm{ℎ𝑎𝑑𝑟𝑜𝑛𝑠}`$ were presented briefly in . The principle purpose of this paper is to explain in detail some of the most important of these techniques. In the numerical/analytical method, one has to work hard to implement the cancellation of “collinear” and “soft” divergences that occur in the integrations. In the numerical method, as we will see, this cancellation happens automatically. On the other hand, in the completely numerical method one has the complication of having to deform some of the integration contours into the complex plane. We will see how to do this deformation. In both the numerical/analytical method and the completely numerical method, one must arrange that the density of integration points is singular near a soft gluon singularity of the integrand (even after cancellations). However, the precise behavior of the densities needed in the two cases is different. We will see what is needed in the numerical method.

These techniques are presented in Secs. II-VI. They are illustrated in Sec. VII with a numerical example. Although a full understanding of the example requires the preceding sections, the reader may want to look briefly at Sec. VII before starting on Secs. II-VI. A brief summary of techniques not presented in detail in this paper is given in Sec. VIII.

In , I presented results from a concrete implementation of the numerical method in computer code. Since then, one logical error in the code has been discovered and fixed and the performance of the program has been improved. Results from the improved code  are presented in Sec. IX.

Let us begin with a precise statement of the problem. We consider an observable such as a particular moment of the thrust distribution. The observable can be expanded in powers of $`\alpha _s/\pi `$,

$$\sigma =\underset{n}{}\sigma ^{[n]},\sigma ^{[n]}\left(\alpha _s/\pi \right)^n.$$
(1)   
The order $`\alpha _s^2`$ contribution has the form

$`\sigma ^{[2]}`$ $`=`$ $`{\displaystyle \frac{1}{2!}}{\displaystyle 𝑑\stackrel{}{k}_1𝑑\stackrel{}{k}_2\frac{d\sigma _2^{[2]}}{d\stackrel{}{k}_1d\stackrel{}{k}_2}𝒮_2(\stackrel{}{k}_1,\stackrel{}{k}_2)}`$ (4)   
$`+{\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)}.`$   
Here 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 $`\overline{\mathrm{MS}}`$ scheme. The functions $`𝒮`$ describe the measurable quantity to be calculated. We wish to calculate a “three-jet-like” quantity. That is, $`𝒮_2=0`$. The normalization is such that $`𝒮_n=1`$ for $`n=2,3,4`$ would give the order $`\alpha _s^2`$ perturbative contribution the the total cross section. There are, of course, infrared divergences associated with Eq. (4). For now, we may simply suppose that an infrared cutoff has been supplied.

The measurement, as specified by the functions $`𝒮_n`$, is to be infrared safe, as described in Ref. : 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)$$
(5)   
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 the functions $`𝒮_n`$ be dimensionless and eliminate the remaining dimensionality in the problem by dividing by $`\sigma _0`$, the total $`e^+e^{}`$ cross section at the Born level. Let us also remove the factor of $`(\alpha _s/\pi )^2`$. Thus, we calculate

$$=\frac{\sigma ^{[2]}}{\sigma _0(\alpha _s/\pi )^2}.$$
(6)   
Our problem is thus to calculate $``$. Let us now see how to set up this problem in a convenient form. 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}),$$
(7)   
where $`h`$ is any function with

$$_0^{\mathrm{}}𝑑\sqrt{s}h(\sqrt{s})=1.$$
(8)   
The quantity $``$ can be expressed in terms of cut Feynman diagrams, as in Fig. 1. The dots where the parton lines cross the cut represent the function $`𝒮_n(\stackrel{}{k}_1,\mathrm{},\stackrel{}{k}_n)`$. Each diagram is a three loop diagram, so we have integrations over loop momenta $`\mathrm{}_1^\mu `$, $`\mathrm{}_2^\mu `$ and $`\mathrm{}_3^\mu `$. 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 (8) 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. We perform this integration by closing the integration contour in the lower half $`E`$ plane. This gives a sum of terms obtained from the original integrand by some algebraic substitutions, as we will see in the following sections. Having performed the energy integrations, we are left with an integral of the form

$$=𝑑\stackrel{}{\mathrm{}}_1𝑑\stackrel{}{\mathrm{}}_2𝑑\stackrel{}{\mathrm{}}_3\underset{G}{}\underset{C}{}g(G,C;\stackrel{}{\mathrm{}}_1,\stackrel{}{\mathrm{}}_2,\stackrel{}{\mathrm{}}_3).$$
(9)   
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 of a given graph.

The problem of calculating $``$ is now set up in a convenient form for calculation. If we were using the Ellis-Ross-Terrano method, we would calculate some of the integrals in Eq. (9) numerically and some analytically. In the method described here, we first perform certain contour deformations, then calculate all of the integrals by Monte Carlo numerical integration. In the following sections, we will learn the main techniques for performing the integrations in Eq. (9). We will do this by studying a simple model problem that will enable us to see the essential features of the numerical method with as few extraneous difficulties as possible.

## II A simplified model

In the following sections, we consider a simplified model in which all complications that are not needed for a first understanding of the numerical method are stripped away. The model is represented by the graph shown in Fig. 2. There are contributions from all of the two and three parton cuts of this diagram, as shown in Fig. 3. Since QCD numerator functions do not play a major role, we consider this graph in $`\varphi ^3`$ theory. Thus, also, we can avoid the complications of ultraviolet renormalization. We consider the incoming momentum $`\stackrel{}{q}`$ to be fixed and nonzero. We calculate the integral of the graph over the incoming energy $`q^0`$. This is analogous to the technical trick of integrating over $`\sqrt{s}`$ in the full three loop QCD calculation (see Sec. I) and serves to provide three energy integrations to perform against three mass-shell delta functions for the three-parton cuts.

We need a nontrivial measurement function $`𝒮`$. As an example, we choose to measure the transverse energy in the final state normalized to the total energy:

$`𝒮_2(\stackrel{}{k}_1,\stackrel{}{k}_2)`$ $`=`$ $`(|\stackrel{}{k}_{T,1}|+|\stackrel{}{k}_{T,2}|)/(|\stackrel{}{k}_1|+|\stackrel{}{k}_2|)`$ (10)   
$`𝒮_3(\stackrel{}{k}_1,\stackrel{}{k}_2,\stackrel{}{k}_3)`$ $`=`$ $`(|\stackrel{}{k}_{T,1}|+|\stackrel{}{k}_{T,2}|+|\stackrel{}{k}_{T,3}|)/(|\stackrel{}{k}_1|+|\stackrel{}{k}_2|+|\stackrel{}{k}_3|),`$ (11)   
where $`\stackrel{}{k}_{T,j}`$ is the part of the momentum $`\stackrel{}{k}_j`$ of the $`j`$th final state particle that is orthogonal to $`\stackrel{}{q}`$.

There are two loops in our diagram. We choose the independent loop momenta to be $`\mathrm{}_2^\mu `$ and $`\mathrm{}_4^\mu `$. The other momenta are understood to be expressed in terms of $`\mathrm{}_2^\mu `$, $`\mathrm{}_4^\mu `$, and $`q^\mu `$.

Thus the example integral that we seek to calculate is

$$=\frac{g^4}{2}\frac{dq^0}{2\pi }\frac{d^4\mathrm{}_2}{(2\pi )^4}\frac{d^4\mathrm{}_4}{(2\pi )^4}𝒲.$$
(12)   
Here $`g`$ is the coupling, $`1/2`$ is the statistical factor for this graph, and the integrand $`𝒲`$ consists of four parts, one for each of the cuts in Fig. 3:

$$𝒲=𝒲_a+𝒲_b+𝒲_c+𝒲_d,$$
(13)   
where

$`𝒲_a`$ $`=`$ $`i𝒮_2(\stackrel{}{\mathrm{}}_4,\stackrel{}{\mathrm{}}_5){\displaystyle \frac{1}{\mathrm{}_1^2+iϵ}}{\displaystyle \frac{1}{\mathrm{}_2^2+iϵ}}{\displaystyle \frac{1}{\mathrm{}_3^2+iϵ}}(2\pi )\mathrm{\Delta }(\mathrm{}_4)(2\pi )\mathrm{\Delta }(\mathrm{}_5),`$ (14)   
$`𝒲_b`$ $`=`$ $`i𝒮_2(\stackrel{}{\mathrm{}}_1,\stackrel{}{\mathrm{}}_3)(2\pi )\mathrm{\Delta }(\mathrm{}_1){\displaystyle \frac{1}{\mathrm{}_2^2iϵ}}(2\pi )\mathrm{\Delta }(\mathrm{}_3){\displaystyle \frac{1}{\mathrm{}_4^2iϵ}}{\displaystyle \frac{1}{\mathrm{}_5^2iϵ}},`$ (15)   
$`𝒲_c`$ $`=`$ $`𝒮_3(\stackrel{}{\mathrm{}}_1,\stackrel{}{\mathrm{}}_2,\stackrel{}{\mathrm{}}_5)(2\pi )\mathrm{\Delta }(\mathrm{}_1)(2\pi )\mathrm{\Delta }(\mathrm{}_2){\displaystyle \frac{1}{\mathrm{}_3^2}}{\displaystyle \frac{1}{\mathrm{}_4^2}}(2\pi )\mathrm{\Delta }(\mathrm{}_5),`$ (16)   
$`𝒲_d`$ $`=`$ $`𝒮_3(\stackrel{}{\mathrm{}}_4,\stackrel{}{\mathrm{}}_2,\stackrel{}{\mathrm{}}_3){\displaystyle \frac{1}{\mathrm{}_1^2}}(2\pi )\mathrm{\Delta }(\mathrm{}_2)(2\pi )\mathrm{\Delta }(\mathrm{}_3)(2\pi )\mathrm{\Delta }(\mathrm{}_4){\displaystyle \frac{1}{\mathrm{}_5^2}}.`$ (17)   
Here we have used the notation

$$\mathrm{\Delta }(k)=\delta (k^2)\theta (k^0).$$
(18)   
## III The integration over energies

We begin by performing the integrals over the energies in Eq. (12). In the case of three partons in the final state, the three delta functions eliminate the three integrations. In the case of two partons in the final state, the two delta functions eliminate two of the energy integrations. This leaves one integration over the energy that circulates around the virtual loop. There are three poles in the upper half plane and three in the lower half plane. Closing the contour in one half plane or the other gives three contributions. Each of these contributions corresponds to putting one of the particles in the loop on shell. Thus altogether there are eight contributions to $``$, as indicated in Fig. 4. We write $``$ as

$$=\frac{g^4}{2(2\pi )^6}𝑑\stackrel{}{\mathrm{}}_4𝑑\stackrel{}{\mathrm{}}_2𝒢,$$
(19)   
where the integrand $`𝒢`$ has eight parts:

$$𝒢=𝒢_{a1}+𝒢_{a2}+𝒢_{a3}+𝒢_{b5}+𝒢_{b2}+𝒢_{b4}+𝒢_c+𝒢_d.$$
(20)   
The contributions to $`𝒢`$ are

$`𝒢_{a1}`$ $`=`$ $`𝒮_2(\stackrel{}{\mathrm{}}_4,\stackrel{}{\mathrm{}}_5){\displaystyle \frac{1}{2|\stackrel{}{\mathrm{}}_1|}}{\displaystyle \frac{1}{(|\stackrel{}{\mathrm{}}_1||\stackrel{}{\mathrm{}}_4|)^2\stackrel{}{\mathrm{}}_2^{\mathrm{\hspace{0.17em}2}}+iϵ}}{\displaystyle \frac{1}{(|\stackrel{}{\mathrm{}}_1||\stackrel{}{\mathrm{}}_4||\stackrel{}{\mathrm{}}_5|)^2\stackrel{}{\mathrm{}}_3^{\mathrm{\hspace{0.17em}2}}+iϵ}}{\displaystyle \frac{1}{2|\stackrel{}{\mathrm{}}_4|}}{\displaystyle \frac{1}{2|\stackrel{}{\mathrm{}}_5|}},`$ (21)   
$`𝒢_{a2}`$ $`=`$ $`𝒮_2(\stackrel{}{\mathrm{}}_4,\stackrel{}{\mathrm{}}_5){\displaystyle \frac{1}{(|\stackrel{}{\mathrm{}}_2|+|\stackrel{}{\mathrm{}}_4|)^2\stackrel{}{\mathrm{}}_1^{\mathrm{\hspace{0.17em}2}}+iϵ}}{\displaystyle \frac{1}{2|\stackrel{}{\mathrm{}}_2|}}{\displaystyle \frac{1}{(|\stackrel{}{\mathrm{}}_2||\stackrel{}{\mathrm{}}_5|)^2\stackrel{}{\mathrm{}}_3^{\mathrm{\hspace{0.17em}2}}+iϵ}}{\displaystyle \frac{1}{2|\stackrel{}{\mathrm{}}_4|}}{\displaystyle \frac{1}{2|\stackrel{}{\mathrm{}}_5|}},`$ (22)   
$`𝒢_{a3}`$ $`=`$ $`𝒮_2(\stackrel{}{\mathrm{}}_4,\stackrel{}{\mathrm{}}_5){\displaystyle \frac{1}{(|\stackrel{}{\mathrm{}}_3|+|\stackrel{}{\mathrm{}}_4|+|\stackrel{}{\mathrm{}}_5|)^2\stackrel{}{\mathrm{}}_1^{\mathrm{\hspace{0.17em}2}}+iϵ}}{\displaystyle \frac{1}{(|\stackrel{}{\mathrm{}}_3|+|\stackrel{}{\mathrm{}}_5|)^2\stackrel{}{\mathrm{}}_2^{\mathrm{\hspace{0.17em}2}}+iϵ}}{\displaystyle \frac{1}{2|\stackrel{}{\mathrm{}}_3|}}{\displaystyle \frac{1}{2|\stackrel{}{\mathrm{}}_4|}}{\displaystyle \frac{1}{2|\stackrel{}{\mathrm{}}_5|}},`$ (23)   
$`𝒢_{b5}`$ $`=`$ $`𝒮_2(\stackrel{}{\mathrm{}}_1,\stackrel{}{\mathrm{}}_3){\displaystyle \frac{1}{2|\stackrel{}{\mathrm{}}_1|}}{\displaystyle \frac{1}{(|\stackrel{}{\mathrm{}}_3|+|\stackrel{}{\mathrm{}}_5|)^2\stackrel{}{\mathrm{}}_2^{\mathrm{\hspace{0.17em}2}}iϵ}}{\displaystyle \frac{1}{2|\stackrel{}{\mathrm{}}_3|}}{\displaystyle \frac{1}{(|\stackrel{}{\mathrm{}}_1|+|\stackrel{}{\mathrm{}}_3|+|\stackrel{}{\mathrm{}}_5|)^2\stackrel{}{\mathrm{}}_4^{\mathrm{\hspace{0.17em}2}}iϵ}}{\displaystyle \frac{1}{2|\stackrel{}{\mathrm{}}_5|}},`$ (24)   
$`𝒢_{b2}`$ $`=`$ $`𝒮_2(\stackrel{}{\mathrm{}}_1,\stackrel{}{\mathrm{}}_3){\displaystyle \frac{1}{2|\stackrel{}{\mathrm{}}_1|}}{\displaystyle \frac{1}{2|\stackrel{}{\mathrm{}}_2|}}{\displaystyle \frac{1}{2|\stackrel{}{\mathrm{}}_3|}}{\displaystyle \frac{1}{(|\stackrel{}{\mathrm{}}_1|+|\stackrel{}{\mathrm{}}_2|)^2\stackrel{}{\mathrm{}}_4^{\mathrm{\hspace{0.17em}2}}iϵ}}{\displaystyle \frac{1}{(|\stackrel{}{\mathrm{}}_3||\stackrel{}{\mathrm{}}_2|)^2\stackrel{}{\mathrm{}}_5^{\mathrm{\hspace{0.17em}2}}iϵ}},`$ (25)   
$`𝒢_{b4}`$ $`=`$ $`𝒮_2(\stackrel{}{\mathrm{}}_1,\stackrel{}{\mathrm{}}_3){\displaystyle \frac{1}{2|\stackrel{}{\mathrm{}}_1|}}{\displaystyle \frac{1}{(|\stackrel{}{\mathrm{}}_1||\stackrel{}{\mathrm{}}_4|)^2\stackrel{}{\mathrm{}}_2^{\mathrm{\hspace{0.17em}2}}iϵ}}{\displaystyle \frac{1}{2|\stackrel{}{\mathrm{}}_3|}}{\displaystyle \frac{1}{2|\stackrel{}{\mathrm{}}_4|}}{\displaystyle \frac{1}{(|\stackrel{}{\mathrm{}}_1|+|\stackrel{}{\mathrm{}}_3||\stackrel{}{\mathrm{}}_4|)^2\stackrel{}{\mathrm{}}_5^{\mathrm{\hspace{0.17em}2}}iϵ}},`$ (26)   
$`𝒢_c`$ $`=`$ $`𝒮_3(\stackrel{}{\mathrm{}}_1,\stackrel{}{\mathrm{}}_2,\stackrel{}{\mathrm{}}_5){\displaystyle \frac{1}{2|\stackrel{}{\mathrm{}}_1|}}{\displaystyle \frac{1}{2|\stackrel{}{\mathrm{}}_2|}}{\displaystyle \frac{1}{(|\stackrel{}{\mathrm{}}_2|+|\stackrel{}{\mathrm{}}_5|)^2\stackrel{}{\mathrm{}}_3^{\mathrm{\hspace{0.17em}2}}}}{\displaystyle \frac{1}{(|\stackrel{}{\mathrm{}}_1|+|\stackrel{}{\mathrm{}}_2|)^2\stackrel{}{\mathrm{}}_4^{\mathrm{\hspace{0.17em}2}}}}{\displaystyle \frac{1}{2|\stackrel{}{\mathrm{}}_5|}},`$ (27)   
$`𝒢_d`$ $`=`$ $`𝒮_3(\stackrel{}{\mathrm{}}_4,\stackrel{}{\mathrm{}}_2,\stackrel{}{\mathrm{}}_3){\displaystyle \frac{1}{(|\stackrel{}{\mathrm{}}_2|+|\stackrel{}{\mathrm{}}_4|)^2\stackrel{}{\mathrm{}}_1^{\mathrm{\hspace{0.17em}2}}}}{\displaystyle \frac{1}{2|\stackrel{}{\mathrm{}}_2|}}{\displaystyle \frac{1}{2|\stackrel{}{\mathrm{}}_3|}}{\displaystyle \frac{1}{2|\stackrel{}{\mathrm{}}_4|}}{\displaystyle \frac{1}{(|\stackrel{}{\mathrm{}}_2|+|\stackrel{}{\mathrm{}}_3|)^2\stackrel{}{\mathrm{}}_5^{\mathrm{\hspace{0.17em}2}}}}.`$ (28)   
So far, the operations that we have performed have been purely algebraic. They are evidently of a sort that can be easily implemented in a computer program in an automatic fashion. We are left with an integral over the loop momenta $`\stackrel{}{\mathrm{}}_2`$ and $`\stackrel{}{\mathrm{}}_4`$. We seek to perform this integration numerically. However, the integrand $`𝒢`$ has singularities, so it is not completely self-evident how to proceed. It is to this question that we now turn.

## IV Cancellation of singularities

In this section, we discuss the cancellation of singularities in a numerical calculation of the integral in Eq. (19).

Let us concentrate to begin with on the cut shown in Fig. 3(a). Then there is a virtual loop consisting of the propagators with momentum labels $`\mathrm{}_1`$, $`\mathrm{}_2`$ and $`\mathrm{}_3`$. Recall that we are taking $`\stackrel{}{\mathrm{}}_2`$, and $`\stackrel{}{\mathrm{}}_4`$ as the independent loop momenta. Put the integration over $`\stackrel{}{\mathrm{}}_2`$ inside the integration over $`\stackrel{}{\mathrm{}}_4`$. Then we can consider $`\stackrel{}{\mathrm{}}_4`$ as fixed while $`\stackrel{}{\mathrm{}}_2`$ varies. Fig. 5 illustrates the space of the loop momentum $`\stackrel{}{\mathrm{}}_2`$ for a particular choice of $`\stackrel{}{q}`$ and at a particular point in the integration over $`\stackrel{}{\mathrm{}}_4`$. The origin of coordinates is at the point labeled $`\stackrel{}{\mathrm{}}_2=0`$. The vector $`\stackrel{}{\mathrm{}}_4`$ is indicated as an arrow with its head at $`\stackrel{}{\mathrm{}}_2=0`$. Then the point $`\stackrel{}{\mathrm{}}_1=0`$ is at the tail of this vector, as indicated. The vector $`\stackrel{}{\mathrm{}}_5=\stackrel{}{q}\stackrel{}{\mathrm{}}_4`$ is indicated as an arrow with its tail at $`\stackrel{}{\mathrm{}}_2=0`$. Then the point $`\stackrel{}{\mathrm{}}_3=0`$ is at the head of this vector, as indicated. Finally, the vector $`\stackrel{}{q}`$ is indicated as an arrow with its tail at $`\stackrel{}{\mathrm{}}_1=0`$.

Where are the singularities of the integrand for our graph?

There is, first of all, a singularity when the momentum of any propagator vanishes since there is always a contribution in which that propagator is put on-shell, with a singularity $`1/(2|\stackrel{}{\mathrm{}}|)`$. Since an integration $`𝑑\stackrel{}{\mathrm{}}/(2|\stackrel{}{\mathrm{}}|)`$ is convergent in the infrared by two powers, these singularities do not cause much difficulty. We simply have to choose a density of points with a matching $`1/|\stackrel{}{\mathrm{}}|`$ singularity, as described later in Sec. VI. We do not discuss these singularities further in this section.

The singularities of concern to us here are

1) A collinear singularity at $`\stackrel{}{\mathrm{}}_2=x\stackrel{}{\mathrm{}}_4`$ with $`0<x<1`$.

2) A collinear singularity at $`\stackrel{}{\mathrm{}}_2=x\stackrel{}{\mathrm{}}_5`$ with $`0<x<1`$.

3) A soft singularity at $`\stackrel{}{\mathrm{}}_2=0`$.

4) A scattering singularity at $`|\stackrel{}{\mathrm{}}_1|+|\stackrel{}{\mathrm{}}_3|=|\stackrel{}{\mathrm{}}_4|+|\stackrel{}{\mathrm{}}_5|`$.

The locations of these singularities are indicated in Fig. 6.

### A The collinear singularities

In this subsection, we examine the collinear singularity at $`\stackrel{}{\mathrm{}}_2=x\stackrel{}{\mathrm{}}_4`$ with $`0<x<1`$. The principles that we discover for this case will hold for the other collinear singularities as well.

The terms $`𝒢_{a1}`$ and $`𝒢_c`$ in the integrand $`𝒢`$, Eq. (20), are singular along the line $`\stackrel{}{\mathrm{}}_2=x\stackrel{}{\mathrm{}}_4`$, $`0<x<1`$. In order to examine this singularity, let us write $`𝒢_{a1}`$ and $`𝒢_c`$ as given in Eq. (28) in the form

$$𝒢_{a1}=\frac{1}{2|\stackrel{}{\mathrm{}}_2+\stackrel{}{\mathrm{}}_4|}\frac{1}{\left(E_2^{(a1)}\right)^2\stackrel{}{\mathrm{}}_2^{\mathrm{\hspace{0.17em}2}}}\frac{1}{2|\stackrel{}{\mathrm{}}_4|}(E_1,E_2^{(a1)},E_5,\stackrel{}{\mathrm{}}_2,\stackrel{}{\mathrm{}}_4)𝒮_2(\stackrel{}{\mathrm{}}_4,\stackrel{}{q}\stackrel{}{\mathrm{}}_4),$$
(29)   
$$𝒢_c=\frac{1}{2|\stackrel{}{\mathrm{}}_2+\stackrel{}{\mathrm{}}_4|}\frac{1}{2|\stackrel{}{\mathrm{}}_2|}\frac{1}{\left(E_1E_2^{(c)}\right)^2\stackrel{}{\mathrm{}}_4^{\mathrm{\hspace{0.17em}2}}}(E_1,E_2^{(c)},E_5,\stackrel{}{\mathrm{}}_2,\stackrel{}{\mathrm{}}_4)𝒮_3(\stackrel{}{\mathrm{}}_1,\stackrel{}{\mathrm{}}_2,\stackrel{}{q}\stackrel{}{\mathrm{}}_4).$$
(30)   
Here the first factors exhibit the denominators for the three propagators that carry collinear momenta at the singularity, $``$ denotes the rest of the Feynman graph, and the $`𝒮`$ functions are the measurement functions for the final state particles. The functions $``$ depend on the loop momenta $`\stackrel{}{\mathrm{}}_2`$ and $`\stackrel{}{\mathrm{}}_4`$ and on three loop energies, which we take to be $`E_1=\mathrm{}_1^0`$, $`E_2=\mathrm{}_2^0`$ and $`E_5=\mathrm{}_5^0`$. The energies are determined by the on-shell delta functions for the two contributions. For $`E_1`$ and $`E_5`$, the values are the same for the two contributions:

$`E_1`$ $`=`$ $`|\stackrel{}{\mathrm{}}_2+\stackrel{}{\mathrm{}}_4|,`$ (31)   
$`E_5`$ $`=`$ $`|\stackrel{}{q}\stackrel{}{\mathrm{}}_4|.`$ (32)   
For $`E_2`$, the values are different:

$`E_2^{(a1)}`$ $`=`$ $`|\stackrel{}{\mathrm{}}_2+\stackrel{}{\mathrm{}}_4||\stackrel{}{\mathrm{}}_4|,`$ (33)   
$`E_2^{(c)}`$ $`=`$ $`|\stackrel{}{\mathrm{}}_2|.`$ (34)   
In order to examine the behavior of $`𝒢_{a1}`$ and $`𝒢_c`$ near the singularity, let

$$\stackrel{}{\mathrm{}}_2=x\stackrel{}{\mathrm{}}_4+\stackrel{}{\mathrm{}}_T,$$
(35)   
where $`\stackrel{}{\mathrm{}}_T\stackrel{}{\mathrm{}}_4=0`$. The singularity is at $`\stackrel{}{\mathrm{}}_T0`$.

In $`𝒢_{a1}`$ the denominator $`\left(E_2^{(a1)}\right)^2\stackrel{}{\mathrm{}}_2^{\mathrm{\hspace{0.17em}2}}`$ vanishes at $`\stackrel{}{\mathrm{}}_T0`$:

$$\left(E_2^{(a1)}\right)^2\stackrel{}{\mathrm{}}_2^{\mathrm{\hspace{0.17em}2}}=\frac{\stackrel{}{\mathrm{}}_T^{\mathrm{\hspace{0.17em}2}}}{1x}\left(1+𝒪(\stackrel{}{\mathrm{}}_T^{\mathrm{\hspace{0.17em}2}})\right).$$
(36)   
Thus there is a $`1/\stackrel{}{\mathrm{}}_T^{\mathrm{\hspace{0.17em}2}}`$ singularity which would give a logarithmically divergent result for the integral of $`𝒢_{a1}`$ alone. Altogether, the denominator factors for $`𝒢_{a1}`$ are

$$\frac{1}{2|\stackrel{}{\mathrm{}}_2+\stackrel{}{\mathrm{}}_4|}\frac{1}{\left(E_2^{(a1)}\right)^2\stackrel{}{\mathrm{}}_2^{\mathrm{\hspace{0.17em}2}}}\frac{1}{2|\stackrel{}{\mathrm{}}_4|}=\frac{1}{4\stackrel{}{\mathrm{}}_4^{\mathrm{\hspace{0.17em}2}}}\frac{1}{\stackrel{}{\mathrm{}}_T^{\mathrm{\hspace{0.17em}2}}}\left(1+𝒪(\stackrel{}{\mathrm{}}_T^{\mathrm{\hspace{0.17em}2}})\right).$$
(37)   
Let us now look at the denominator factors for $`𝒢_c`$. The denominator $`\left(E_1E_2^{(c)}\right)^2\stackrel{}{\mathrm{}}_4^{\mathrm{\hspace{0.17em}2}}`$ takes the form

$$\left(\left(E_1E_2^{(c)}\right)^2\stackrel{}{\mathrm{}}_4^{\mathrm{\hspace{0.17em}2}}\right)^2=\frac{\stackrel{}{\mathrm{}}_T^{\mathrm{\hspace{0.17em}2}}}{x(1x)}\left(1+𝒪(\stackrel{}{\mathrm{}}_T^{\mathrm{\hspace{0.17em}2}})\right),$$
(38)   
so that the denominator factors together take the form

$$\frac{1}{2|\stackrel{}{\mathrm{}}_2+\stackrel{}{\mathrm{}}_4|}\frac{1}{2|\stackrel{}{\mathrm{}}_2|}\frac{1}{\left(E_1E_2^{(c)}\right)^2\stackrel{}{\mathrm{}}_4^{\mathrm{\hspace{0.17em}2}}}=\frac{1}{4\stackrel{}{\mathrm{}}_4^{\mathrm{\hspace{0.17em}2}}}\frac{1}{\stackrel{}{\mathrm{}}_T^{\mathrm{\hspace{0.17em}2}}}\left(1+𝒪(\stackrel{}{\mathrm{}}_T^{\mathrm{\hspace{0.17em}2}})\right).$$
(39)   
Again, we have a $`1/\stackrel{}{\mathrm{}}_T^{\mathrm{\hspace{0.17em}2}}`$ singularity.

Note, however, that the denominator factors in Eqs. (37) and (39) are equal except for their sign, up to corrections that are not singular as $`\stackrel{}{\mathrm{}}_T^{\mathrm{\hspace{0.17em}2}}0`$. Thus if the remaining factors $``$ and $`𝒮`$ were exactly the same for $`𝒢_{a1}`$ and $`𝒢_c`$ there would be no singularity in their sum.

We thus need to explore the matching of $``$ and $`𝒮`$. The two versions of $``$ are the same functions with the same arguments except for the fact that $`E_2^{(c)}E_2^{(a1)}`$. However,

$$E_2^{(c)}=E_2^{(a1)}+𝒪(\stackrel{}{\mathrm{}}_T^{\mathrm{\hspace{0.17em}2}}).$$
(40)   
Thus

$$(E_1,E_2^{(c)},E_5,\stackrel{}{\mathrm{}}_2,\stackrel{}{\mathrm{}}_4)=(E_1,E_2^{(a1)},E_5,\stackrel{}{\mathrm{}}_2,\stackrel{}{\mathrm{}}_4)+𝒪(\stackrel{}{\mathrm{}}_T^{\mathrm{\hspace{0.17em}2}}).$$
(41)   
For the functions $`𝒮`$ used in our example, we have

$$𝒮_3((1x)\stackrel{}{\mathrm{}}_4+\stackrel{}{\mathrm{}}_T,x\stackrel{}{\mathrm{}}_4\stackrel{}{\mathrm{}}_T,\stackrel{}{q}\stackrel{}{\mathrm{}}_4)=𝒮_2(\stackrel{}{\mathrm{}}_4,\stackrel{}{q}\stackrel{}{\mathrm{}}_4)+𝒪(\stackrel{}{\mathrm{}}_T^{\mathrm{\hspace{0.17em}2}}).$$
(42)   
Using these matching equations we find that

$$𝒢_{a1}+𝒢_c=𝒪(1)$$
(43)   
as $`\stackrel{}{\mathrm{}}_T0`$. There is no collinear singularity in $`𝒢`$.

How general is this result? First of all note that, in the part of the argument not involving the measurement functions $`𝒮`$, we used only the explicit structure of the denominators for three propagators that meet at a vertex. In the limit in which the momenta carried on these propagators become collinear, there is a cancellation of the collinear singularity arising from these denominators. The three propagators can be part of a much larger graph, and there can be non-trivial numerator factors, as in QCD. All of the other factors can be lumped into a function $``$ and treated as above. Thus this cancellation works in QCD as well as $`\varphi ^3`$ theory and it works for cut graphs with at most one virtual loop at any order of perturbation theory.

As for the measurement functions, in general we need to consider the difference between the measurement functions with $`n`$ and $`n+1`$ particles in the final state,

$$F(\stackrel{}{\mathrm{}}_T)=𝒮_{n+1}(\stackrel{}{k}_1,\mathrm{},\stackrel{}{k}_{n1},x\stackrel{}{k}_n\stackrel{}{\mathrm{}}_T,(1x)\stackrel{}{k}_n+\stackrel{}{\mathrm{}}_T)𝒮_n(\stackrel{}{k}_1,\mathrm{},\stackrel{}{k}_{n1},\stackrel{}{k}_n).$$
(44)   
Assuming that $`F`$ is an analytic function of $`\stackrel{}{\mathrm{}}_T`$, it will have an expansion around $`\stackrel{}{\mathrm{}}_T=0`$ of the form

$$F(\stackrel{}{\mathrm{}}_T)=a+b_i\mathrm{}_T^i+c_{ij}\mathrm{}_T^i\mathrm{}_T^j+\mathrm{}.$$
(45)   
Infrared safety requires that $`a=0`$. If $`b_i0`$ then $`F`$ vanishes on a surface that intersects the point $`\stackrel{}{\mathrm{}}_T=0`$. Measurement functions $`𝒮`$ with this property would define an infrared safe measurement, but I do not know of any example in common use. More typically, $`F`$ is non-zero in a neighborhood of $`\stackrel{}{\mathrm{}}_T=0`$ while vanishing at $`\stackrel{}{\mathrm{}}_T=0`$. Then both $`a`$ and the $`b_i`$ must vanish and the $`c_{ij}`$ should be a positive definite (or negative definite) matrix. Thus, for typical measurement functions,

$$F(\stackrel{}{\mathrm{}}_T)=𝒪(\stackrel{}{\mathrm{}}_T^{\mathrm{\hspace{0.17em}2}})$$
(46)   
as $`\stackrel{}{\mathrm{}}_T0`$. Then the integrand does not have collinear singularities.

For an atypical measurement function with $`b_i0`$, one would be left with an integrable singularity of the form $`\stackrel{}{b}\stackrel{}{\mathrm{}}_T/\stackrel{}{\mathrm{}}_T^{\mathrm{\hspace{0.17em}2}}`$. The current version of the computer code  has a mechanism to deal with this contingency, but I do not discuss it here since I know of no case in which it is needed.

### B The soft singularities

In this subsection, we examine the soft singularity at $`\stackrel{}{\mathrm{}}_2=0`$.

Let us concentrate to begin with on the cut graph shown in Fig. 3(a). When we perform the integration over the energy circulating in the virtual loop, there is a contribution from the term in which the propagator carrying momentum $`\mathrm{}_1^\mu `$ is put on shell, as in Fig. 4(a1). This contribution is $`𝒢_{a1}`$ in Eq. (28). Let us examine this contribution in the limit $`\stackrel{}{\mathrm{}}_20`$. Expanding in powers of $`\stackrel{}{\mathrm{}}_2`$, we have

$$\mathrm{}_1^0=|\stackrel{}{\mathrm{}}_1|=|\stackrel{}{\mathrm{}}_4+\stackrel{}{\mathrm{}}_2|=|\stackrel{}{\mathrm{}}_4|+|\stackrel{}{\mathrm{}}_2|\stackrel{}{u}_2\stackrel{}{u}_4+\mathrm{},$$
(47)   
where we adopt the notation

$$\stackrel{}{u}_J=\stackrel{}{\mathrm{}}_J/|\stackrel{}{\mathrm{}}_J|.$$
(48)   
Then

$$\mathrm{}_2^2=|\stackrel{}{\mathrm{}}_2|^2[1(\stackrel{}{u}_2\stackrel{}{u}_4)^2]+\mathrm{}$$
(49)   
and

$$\mathrm{}_3^2=2|\stackrel{}{\mathrm{}}_5||\stackrel{}{\mathrm{}}_2|\stackrel{}{u}_2(\stackrel{}{u}_5\stackrel{}{u}_4)+\mathrm{}.$$
(50)   
Thus

$`𝒢_{a1}`$ $``$ $`𝒮_2{\displaystyle \frac{1}{2|\stackrel{}{\mathrm{}}_4|}}{\displaystyle \frac{1}{|\stackrel{}{\mathrm{}}_2|^2[1(\stackrel{}{u}_2\stackrel{}{u}_4)^2]}}{\displaystyle \frac{1}{2|\stackrel{}{\mathrm{}}_5||\stackrel{}{\mathrm{}}_2|\stackrel{}{u}_2(\stackrel{}{u}_5\stackrel{}{u}_4)+iϵ}}{\displaystyle \frac{1}{2|\stackrel{}{\mathrm{}}_4|}}{\displaystyle \frac{1}{2|\stackrel{}{\mathrm{}}_5|}}`$ (51)   
$`=`$ $`{\displaystyle \frac{𝒮_2}{16|\stackrel{}{\mathrm{}}_4|^2|\stackrel{}{\mathrm{}}_5|^2}}{\displaystyle \frac{1}{|\stackrel{}{\mathrm{}}_2|^3}}{\displaystyle \frac{1}{1(\stackrel{}{u}_2\stackrel{}{u}_4)^2}}{\displaystyle \frac{1}{\stackrel{}{u}_2(\stackrel{}{u}_5\stackrel{}{u}_4)+iϵ}}.`$ (52)   
We proceed in this fashion to evaluate the contribution corresponding to Fig. 4(a2). Then we evaluate the contribution of Fig. 4(a3), but we find that this contribution is not singular as $`\stackrel{}{\mathrm{}}_20`$. Adding the three contributions, we obtain the net integrand for the cut graph of Fig. 3(a) in the soft limit $`\stackrel{}{\mathrm{}}_20`$:

$$𝒢_a\frac{𝒮_2}{32|\stackrel{}{\mathrm{}}_4|^2|\stackrel{}{\mathrm{}}_5|^2}\frac{1}{|\stackrel{}{\mathrm{}}_2|^3}\frac{1}{1+\stackrel{}{u}_2\stackrel{}{u}_4}\frac{1}{1\stackrel{}{u}_2\stackrel{}{u}_5}\frac{2\stackrel{}{u}_2(\stackrel{}{u}_5\stackrel{}{u}_4)}{\stackrel{}{u}_2(\stackrel{}{u}_5\stackrel{}{u}_4)+iϵ}.$$
(53)   
Some comments are in order here. First, we have included the leading term, with a $`1/|\stackrel{}{\mathrm{}}_2|^3`$ singularity, and dropped less singular terms. If we decompose the integration over $`\stackrel{}{\mathrm{}}_2`$ into $`𝑑\mathrm{\Omega }_2|\stackrel{}{\mathrm{}}_2|^2d|\stackrel{}{\mathrm{}}_2|`$, then a $`1/|\stackrel{}{\mathrm{}}_2|^3`$ singularity produces a logarithmic divergence in the integration over $`|\stackrel{}{\mathrm{}}_2|`$. The less singular terms will lead to a finite integration over $`|\stackrel{}{\mathrm{}}_2|`$, although the integration $`𝑑\mathrm{\Omega }`$ over the angles $`\stackrel{}{u}_2`$ can still be divergent. There are, in fact, singularities in the angular integration. The factor $`1/[1\stackrel{}{u}_2\stackrel{}{u}_5]`$ is singular when $`\stackrel{}{\mathrm{}}_2`$ is collinear with $`\stackrel{}{\mathrm{}}_5`$, while the factor $`1/[1+\stackrel{}{u}_2\stackrel{}{u}_4]`$ is singular when $`\stackrel{}{\mathrm{}}_2`$ is collinear with $`\stackrel{}{\mathrm{}}_4`$. These singularities produce logarithmically divergent integrations over $`\stackrel{}{u}_2`$. However, the analysis of the previous subsection shows that the collinear singularities cancel among the cuts of our graph. There is also a singularity on the plane $`\stackrel{}{u}_2(\stackrel{}{u}_5\stackrel{}{u}_4)=0`$. This is the scattering singularity on the ellipse $`|\stackrel{}{\mathrm{}}_1|+|\stackrel{}{\mathrm{}}_3|=|\stackrel{}{\mathrm{}}_4|+|\stackrel{}{\mathrm{}}_5|`$. This ellipse passes through the point $`\stackrel{}{\mathrm{}}_2=0`$ and the plane tangent to the ellipse at this point is the plane $`\stackrel{}{u}_2(\stackrel{}{u}_5\stackrel{}{u}_4)=0`$. I will have more to say about this singularity later. Here we note simply that it comes with an $`iϵ`$ prescription, which has been preserved in Eq. (53).

We now consider the cut graph shown in Fig. 3(b). Again, there are three contributions to consider, corresponding to the diagrams $`(b5)`$, $`(b2)`$ and $`(b4)`$ in Fig. 4. Adding the three contributions, we obtain the net integrand for the cut graph of Fig. 3(b) in the soft limit $`\stackrel{}{\mathrm{}}_20`$:

$$𝒢_b\frac{𝒮_2}{32|\stackrel{}{\mathrm{}}_4|^2|\stackrel{}{\mathrm{}}_5|^2}\frac{1}{|\stackrel{}{\mathrm{}}_2|^3}\frac{1}{1\stackrel{}{u}_2\stackrel{}{u}_4}\frac{1}{1+\stackrel{}{u}_2\stackrel{}{u}_5}\frac{2+\stackrel{}{u}_2(\stackrel{}{u}_5\stackrel{}{u}_4)}{\stackrel{}{u}_2(\stackrel{}{u}_5\stackrel{}{u}_4)+iϵ}.$$
(54)   
As in Eq. (53), there are a scattering singularity and two collinear singularities. However, the signs that indicate the location of the collinear singularities are reversed compared to Eq. (54). If we add $`𝒢_a`$ and $`𝒢_b`$ we obtain

$$𝒢_a+𝒢_b\frac{𝒮_2}{16|\stackrel{}{\mathrm{}}_4|^2|\stackrel{}{\mathrm{}}_5|^2}\frac{1}{|\stackrel{}{\mathrm{}}_2|^3}\frac{1+(\stackrel{}{u}_2\stackrel{}{u}_4)(\stackrel{}{u}_2\stackrel{}{u}_5)}{[1(\stackrel{}{u}_2\stackrel{}{u}_4)^2][1(\stackrel{}{u}_2\stackrel{}{u}_5)^2]}.$$
(55)   
Thus, the overall $`1/|\stackrel{}{\mathrm{}}_2|^3`$ singularity remains and the collinear singularities remain, but the scattering singularities cancel in the soft limit, $`\stackrel{}{\mathrm{}}_20`$, between the two cuts that leave virtual subgraphs.

There are two more cut graphs to consider. The graph shown in Fig. 3(c) gives

$$𝒢_c\frac{𝒮_3}{32|\stackrel{}{\mathrm{}}_4|^2|\stackrel{}{\mathrm{}}_5|^2}\frac{1}{|\stackrel{}{\mathrm{}}_2|^3}\frac{1}{[1+\stackrel{}{u}_2\stackrel{}{u}_4]}\frac{1}{[1+\stackrel{}{u}_2\stackrel{}{u}_5]}.$$
(56)   
The graph shown in Fig. 3(d) gives

$$𝒢_d\frac{𝒮_3}{32|\stackrel{}{\mathrm{}}_4|^2|\stackrel{}{\mathrm{}}_5|^2}\frac{1}{|\stackrel{}{\mathrm{}}_2|^3}\frac{1}{[1\stackrel{}{u}_2\stackrel{}{u}_4]}\frac{1}{[1\stackrel{}{u}_2\stackrel{}{u}_5]}.$$
(57)   
Adding these together, we find

$$𝒢_c+𝒢_d\frac{𝒮_3}{16|\stackrel{}{\mathrm{}}_4|^2|\stackrel{}{\mathrm{}}_5|^2}\frac{1}{|\stackrel{}{\mathrm{}}_2|^3}\frac{1+(\stackrel{}{u}_2\stackrel{}{u}_4)(\stackrel{}{u}_2\stackrel{}{u}_5)}{[1(\stackrel{}{u}_2\stackrel{}{u}_4)^2][1(\stackrel{}{u}_2\stackrel{}{u}_5)^2]}.$$
(58)   
We note that when we add the contributions of the cuts which leave virtual subgraphs to the contributions of the cuts which have no virtual subgraphs, the leading soft singularity cancels:

$$𝒢_a+𝒢_b+𝒢_c+𝒢_d0.$$
(59)   
That is, after cancellation, the overall singularity is at worst proportional to $`1/|\stackrel{}{\mathrm{}}_2|^2`$. It is thus an integrable singularity provided that all of the singularities of the angular integration over $`\stackrel{}{u}_2`$ cause no problems.

The cancellation of the leading soft singularity is built into the structure of Feynman diagrams, so that we do not have to do anything special to make it happen. However, there is a certain subtlety in arranging for the singularities in the angular integrations to be convergent in a Monte Carlo integration. Thus, we will return to the cancellation of the soft singularity after we have discussed contour deformations in the following section.

## V The scattering singularity and contour deformation

Consider the contribution from Fig. 4(a1), as given in Eq. (28). There is a factor

$$\frac{1}{(|\stackrel{}{\mathrm{}}_1||\stackrel{}{\mathrm{}}_4||\stackrel{}{\mathrm{}}_5|)^2\stackrel{}{\mathrm{}}_3^{\mathrm{\hspace{0.17em}2}}+iϵ}=\frac{1}{(|\stackrel{}{\mathrm{}}_3|+|\stackrel{}{\mathrm{}}_4|+|\stackrel{}{\mathrm{}}_5||\stackrel{}{\mathrm{}}_1|)(|\stackrel{}{\mathrm{}}_4|+|\stackrel{}{\mathrm{}}_5||\stackrel{}{\mathrm{}}_1||\stackrel{}{\mathrm{}}_3|+iϵ)},$$
(60)   
which has a singularity when $`|\stackrel{}{\mathrm{}}_1|+|\stackrel{}{\mathrm{}}_3|=|\stackrel{}{\mathrm{}}_4|+|\stackrel{}{\mathrm{}}_5|`$. In an analysis using time-ordered perturbation theory, the singular factor emerges from the energy denominator associated with the intermediate state consisting of partons 1 and 3,

$$E_FE(\stackrel{}{\mathrm{}}_2)+iϵ,$$
(61)   
where $`E_F=|\stackrel{}{\mathrm{}}_4|+|\stackrel{}{\mathrm{}}_5|`$ and

$$E(\stackrel{}{\mathrm{}}_2)=|\stackrel{}{\mathrm{}}_1|+|\stackrel{}{\mathrm{}}_3|=|\stackrel{}{\mathrm{}}_4+\stackrel{}{\mathrm{}}_2|+|\stackrel{}{\mathrm{}}_5+\stackrel{}{\mathrm{}}_2|.$$
(62)   
Thus the singularity appears when the momenta are right for particles 1 and 3 to be on-shell and scatter to produce the final state particles 4 and 5.

The contribution from Fig. 4(b4) has a scattering singularity at the same place as that from the cut diagram (a1). However, these singularities do not cancel in general because the functions $`𝒮_2(\stackrel{}{\mathrm{}}_4,\stackrel{}{\mathrm{}}_5)`$ and $`𝒮_2(\stackrel{}{\mathrm{}}_1,\stackrel{}{\mathrm{}}_3)`$ do not match. We thus have a problem if we would like to perform the integration numerically.

We notice, however, that the singularity is protected by an $`iϵ`$ prescription. The $`iϵ`$ in the denominator tells us what to do in an analytic calculation and it also tells us what to do in a numerical calculation: we need to deform the integration contour.

We are integrating over a loop momentum $`\stackrel{}{\mathrm{}}_2`$. Let us replace $`\stackrel{}{\mathrm{}}_2`$ by a complex momentum $`\stackrel{}{\mathrm{}}_{2,c}=\stackrel{}{\mathrm{}}_2+i\stackrel{}{\kappa }`$, where $`\stackrel{}{\kappa }`$ is a function, which remains to be determined, of $`\stackrel{}{\mathrm{}}_2`$. Then as we integrate over the real vector $`\stackrel{}{\mathrm{}}_2`$, we are integrating over a contour in the space of the complex vector $`\stackrel{}{\mathrm{}}_{2,c}`$. When we deform the original contour $`\stackrel{}{\mathrm{}}_{2,c}=\stackrel{}{\mathrm{}}_2`$ to the new contour $`\stackrel{}{\mathrm{}}_{2,c}=\stackrel{}{\mathrm{}}_2+i\stackrel{}{\kappa }`$, the integral does not change provided that we do not cross any points where the integrand is singular and provided that we include a jacobian

$$𝒥(\stackrel{}{\mathrm{}}_2)=det\left(\frac{\mathrm{}_{2,c}^i}{\mathrm{}_2^j}\right).$$
(63)   
There are some subtleties associated with this; the relevant theorem is proved in the Appendix.

We need to choose $`\stackrel{}{\kappa }`$ as a function of $`\stackrel{}{\mathrm{}}_2`$. Consider first the direction of $`\stackrel{}{\kappa }`$. On the deformed contour, the energy denominator (61) has the form

$$E_FE(\stackrel{}{\mathrm{}}_2+i\stackrel{}{\kappa })+iϵ.$$
(64)   
In order to fix the direction of deformation, it is useful to consider what happens when we deform the contour just a little way from the real $`\stackrel{}{\mathrm{}}_2`$ space. For small $`\kappa `$, we have

$$E(\stackrel{}{\mathrm{}}_2+i\stackrel{}{\kappa })|\stackrel{}{\mathrm{}}_1|+|\stackrel{}{\mathrm{}}_3|+i\stackrel{}{\kappa }\stackrel{}{w},$$
(65)   
where

$$\stackrel{}{w}=\frac{\stackrel{}{\mathrm{}}_1}{|\stackrel{}{\mathrm{}}_1|}+\frac{\stackrel{}{\mathrm{}}_3}{|\stackrel{}{\mathrm{}}_3|}.$$
(66)   
Thus the energy denominator is $`E_FE(\stackrel{}{\mathrm{}}_2)i\stackrel{}{\kappa }\stackrel{}{w}+iϵ`$ for small $`\stackrel{}{\kappa }`$. In order to keep on the proper side of the singularity, we want $`\stackrel{}{\kappa }\stackrel{}{w}`$ to be negative. The simplest way to insure this is to choose $`\stackrel{}{\kappa }`$ in the direction of $`\stackrel{}{w}`$. Thus we choose

$$\stackrel{}{\kappa }=D(\stackrel{}{\mathrm{}}_2)\stackrel{}{w},D(\stackrel{}{\mathrm{}}_2)0.$$
(67)   
Then the singular factor is approximately

$$\frac{1}{E_FE(\stackrel{}{\mathrm{}}_2)+iD(\stackrel{}{\mathrm{}}_2)\stackrel{}{w}^{\mathrm{\hspace{0.17em}2}}}$$
(68)   
for a small deformation. For a larger deformation, it not so simple to see that we stay on the correct side of the singularity, but it is easy to check numerically.

The next question is how should we choose $`D(\stackrel{}{\mathrm{}}_2)`$? We want $`D`$ not to be small when $`\stackrel{}{\mathrm{}}_2`$ is near the surface $`E(\stackrel{}{\mathrm{}}_2)=E_F`$ in order that the integrand not be large there. We want $`D(\stackrel{}{\mathrm{}}_2)`$ not to grow as $`\stackrel{}{\mathrm{}}_2^{\mathrm{\hspace{0.17em}2}}\mathrm{}`$ in order to satisfy the conditions for the theorem that deforming the contour does not change value of the integral. Since there is no reason to keep any finite contour deformation for large $`\stackrel{}{\mathrm{}}_2^{\mathrm{\hspace{0.17em}2}}`$, we will simply choose to have $`D(\stackrel{}{\mathrm{}}_2)0`$ as $`\stackrel{}{\mathrm{}}_2^{\mathrm{\hspace{0.17em}2}}\mathrm{}`$.

There is another condition that $`D(\stackrel{}{\mathrm{}}_2)`$ should obey: it should vanish at points where $`𝒢`$ has collinear and soft singularities. To see why takes some discussion.

Consider the contributions from three parton cuts, for which there is no virtual loop. For these contributions, we do not want to deform the contours. This is because if any of the loop momenta were complex then at least one of the momenta of the final state particles would be complex. In principle, one could have complex momenta for final state particles as long as the measurement functions $`𝒮_n(\stackrel{}{k}_1,\mathrm{},\stackrel{}{k}_n)`$ are analytic. However, I have in mind applications in which the numerical integration program acts as a subroutine that produces “events” with final state particle momenta $`\{\stackrel{}{k}_1,\mathrm{},\stackrel{}{k}_n\}`$ and weights computed by the subroutine. Then the events could be the input to, for example, a Monte Carlo program that generates parton showers and hadronization. Surely complex momenta for the final state particles are not desirable.

Now recall that there is a cancellation among the contributions $`𝒢_C`$ from different cuts $`C`$ at points where the $`𝒢_C`$ have collinear and soft singularities. Evidently, if we deform the contour for a contribution with a virtual graph but do not deform the contour for the canceling contribution, then the cancellation can be spoiled. We can avoid spoiling the cancellation if we make the contours match at the singularity. That is, $`D(\stackrel{}{\mathrm{}}_2)`$ should vanish at the points where the $`𝒢_C`$ have collinear and soft singularities.

We also need to determine how fast $`D(\stackrel{}{\mathrm{}}_2)`$ needs to approach zero as $`\stackrel{}{\mathrm{}}_2`$ approaches a singularity. Since the integration is in a multidimensional complex space, we need an analysis that makes use of the multidimensional contour deformation theorem. This analysis is given in the Appendix. Here, I present a simpler one dimensional analysis that can serve to clarify the issue.

Consider the following toy integral,

$$I=_0^{x_{\mathrm{max}}}\frac{dx}{x}\left\{\frac{f_V(x)}{x1+iϵ}+f_R(x)\right\}.$$
(69)   
Here the endpoint singularity at $`x=0`$ plays the role of the collinear or soft singularities. The function $`f_V/(x1+iϵ)`$ plays the role of the integrand for the contribution with a virtual subgraph. In this contribution, there is a singularity at $`x=1`$ that comes with an $`iϵ`$ prescription. The function $`f_R`$ plays the role of the integrand for the contribution with no virtual subgraph. We assume that $`f_V(z)`$ and $`f_R(z)`$ are analytic functions. We also assume that $`f_V(0)=f_R(0)`$, so that the apparent singularity at $`x=0`$ cancels.

Now the $`iϵ`$ prescription on the singularity at $`x=1`$ tells us that we can deform the integration contour into the upper half plane, replacing $`x`$ by $`z=x+iy(x)`$ where $`y(0)=y(x_{\mathrm{max}})=0`$. Thus

$$I=_0^{x_{\mathrm{max}}}𝑑x\frac{1+iy^{}(x)}{x+iy(x)}\left\{\frac{f_V\left(x+iy(x)\right)}{x1+iy(x)}+f_R\left(x+iy(x)\right)\right\}.$$
(70)   
Suppose, however, that we want to keep the contour for $`f_R`$ on the real axis. Then we might hope that $`I=\stackrel{~}{I}`$, where

$$\stackrel{~}{I}=\underset{x_{\mathrm{min}}0}{lim}_{x_{\mathrm{min}}}^{x_{\mathrm{max}}}𝑑x\left\{\frac{1+iy^{}(x)}{x+iy(x)}\frac{f_V\left(x+iy(x)\right)}{x1+iy(x)}+\frac{f_R\left(x\right)}{x}\right\}.$$
(71)   
The difference is

$$\stackrel{~}{I}I=\underset{x_{\mathrm{min}}0}{lim}_{x_{\mathrm{min}}}^{x_{\mathrm{max}}}𝑑x\left\{\frac{f_R\left(x\right)}{x}[1+iy^{}(x)]\frac{f_R\left(x+iy(x)\right)}{x+iy(x)}\right\}.$$
(72)   
If we note that $`[f_R(z)f_R(0)]/z`$ is an analytic function even at $`z=0`$ and that the integral of an analytic function around a closed contour vanishes, we have

$$0=\underset{x_{\mathrm{min}}0}{lim}_{x_{\mathrm{min}}}^{x_{\mathrm{max}}}𝑑x\left\{\frac{f_R\left(x\right)f_R\left(0\right)}{x}[1+iy^{}(x)]\frac{f_R\left(x+iy(x)\right)f_R\left(0\right)}{x+iy(x)}\right\}.$$
(73)   
Subtracting these and performing the integral, we have

$$\stackrel{~}{I}I=f_R(0)\underset{x_{\mathrm{min}}0}{lim}_{x_{\mathrm{min}}}^{x_{\mathrm{max}}}𝑑x\left\{\frac{1}{x}\frac{1+iy^{}(x)}{x+iy(x)}\right\}=f_R(0)\underset{x_{\mathrm{min}}0}{lim}\mathrm{log}\left(1+i\frac{y(x_{\mathrm{min}})}{x_{\mathrm{min}}}\right).$$
(74)   
We can draw two conclusions. First, as long as $`y(x)0`$ at least as fast as $`x^1`$ as $`x0`$, we will realize the cancellation of the $`x0`$ singularity and obtain a finite value for $`\stackrel{~}{I}`$. Second, if we choose $`y(x)x`$ as $`x0`$, $`\stackrel{~}{I}`$ will be finite, but it will not be equal to the correct result $`I`$. In order to get a result $`\stackrel{~}{I}`$ that is not only finite but also correct, we need $`y(x)/x0`$ as $`x0`$. A convenient choice is $`y(x)x^2`$ as $`x0`$.

We conclude from the multidimensional extension of this analysis, given in the Appendix, that as $`\stackrel{}{\mathrm{}}_2`$ approaches a singularity, $`D(\stackrel{}{\mathrm{}}_2)`$ should approach zero quadratically with the distance to the singularity.

We now use the qualitative criteria just developed to give a specific choice of deformation. We have chosen

$$\stackrel{}{\mathrm{}}_{2,c}=\stackrel{}{\mathrm{}}_2iD(\stackrel{}{\mathrm{}}_2)\stackrel{}{w},$$
(75)   
where $`\stackrel{}{w}`$, Eq. (66), specifies the direction of deformation. We now specify a deformation function $`D(\stackrel{}{\mathrm{}}_2)`$ that satisfies our criteria. We write $`D`$ in the form

$$D=CG.$$
(76)   
The factor $`C`$ is designed to insure that the deformation vanishes quadratically near the collinear and soft singularities. The factor $`G`$ is designed to turn the deformation off for large $`\stackrel{}{\mathrm{}}_2`$. These factors are explained below and are defined precisely in Eqs. (80) and (83) below.

First, we discuss the factor $`C`$. We want the deformation to vanish at the line $`\stackrel{}{\mathrm{}}_2=x\stackrel{}{\mathrm{}}_4`$ with $`0x1`$, where the amplitude has a collinear singularity. (Since $`\stackrel{}{\mathrm{}}_4=\stackrel{}{\mathrm{}}_1\stackrel{}{\mathrm{}}_2`$, this line is also $`\stackrel{}{\mathrm{}}_1=\lambda \stackrel{}{\mathrm{}}_2`$ with $`0<\lambda <\mathrm{}`$.) Define

$$d_{12}=\frac{\left||\stackrel{}{\mathrm{}}_2|\stackrel{}{\mathrm{}}_1+|\stackrel{}{\mathrm{}}_1|\stackrel{}{\mathrm{}}_2\right|}{|\stackrel{}{\mathrm{}}_1\stackrel{}{\mathrm{}}_2|}=\frac{\left||\stackrel{}{\mathrm{}}_2|\stackrel{}{\mathrm{}}_1+|\stackrel{}{\mathrm{}}_1|\stackrel{}{\mathrm{}}_2\right|}{|\stackrel{}{\mathrm{}}_4|}.$$
(77)   
This function is zero on the line $`\stackrel{}{\mathrm{}}_2=x\stackrel{}{\mathrm{}}_4`$ with $`0x1`$, and furthermore, it vanishes linearly as $`\stackrel{}{\mathrm{}}_2`$ approaches this line. Similarly, we want the deformation to vanish on the line $`\stackrel{}{\mathrm{}}_2=x\stackrel{}{\mathrm{}}_5`$ with $`0x1`$, where the amplitude has its other collinear singularity. The function $`d_{23}`$, where

$$d_{23}=\frac{\left||\stackrel{}{\mathrm{}}_3|\stackrel{}{\mathrm{}}_2+|\stackrel{}{\mathrm{}}_2|\stackrel{}{\mathrm{}}_3\right|}{|\stackrel{}{\mathrm{}}_2\stackrel{}{\mathrm{}}_3|}=\frac{\left||\stackrel{}{\mathrm{}}_3|\stackrel{}{\mathrm{}}_2+|\stackrel{}{\mathrm{}}_2|\stackrel{}{\mathrm{}}_3\right|}{|\stackrel{}{\mathrm{}}_5|},$$
(78)   
vanishes linearly as $`\stackrel{}{\mathrm{}}_2`$ approaches this line. (To see this, use $`\stackrel{}{\mathrm{}}_5=\stackrel{}{\mathrm{}}_2\stackrel{}{\mathrm{}}_3`$.) Let

$$d=\mathrm{min}(d_{12},d_{23}).$$
(79)   
Then $`d`$ vanishes linearly with the distance to either of the collinear singularities. It also vanishes linearly with the distance to the soft singularity at $`\stackrel{}{\mathrm{}}_2=0`$.

Now, we have seen that the deformation should vanish quadratically with the distance to any of the singularities. We can achieve this by letting

$$C(d^2)=\frac{\alpha d^2}{1+4\beta d^2/(|\mathrm{}_4|+|\mathrm{}_5|+|\stackrel{}{q}|)^2},$$
(80)   
where $`\alpha `$ and $`\beta `$ are adjustable dimensionless parameters. Note that, for large $`d`$, $`C(d^2)`$ approaches a constant.

Next, we discuss the factor $`G`$. We want to ensure that the contour deformation vanishes for large $`\stackrel{}{\mathrm{}}_2`$. Let us define

$$a=|\stackrel{}{\mathrm{}}_1|+|\stackrel{}{\mathrm{}}_3||\stackrel{}{q}|$$
(81)   
and

$$A=|\stackrel{}{\mathrm{}}_4|+|\stackrel{}{\mathrm{}}_5||\stackrel{}{q}|.$$
(82)   
Then the singularity that we are avoiding by means of contour deformation is at $`a=A`$. We can turn the deformation off for $`aA`$ by setting

$$G(a)=\frac{1}{A+\gamma a},$$
(83)   
where $`\gamma `$ is an adjustable dimensionless parameter.

There is a subsidiary reason for this choice. At the singularity, $`G=1/[(1+\gamma )A]`$. The factor $`1/A`$ serves to enhance the deformation in the case that $`\stackrel{}{\mathrm{}}_4`$ and $`\stackrel{}{\mathrm{}}_5`$ are nearly collinear, in which case $`d`$ is small on the ellipse $`a=A`$ and the deformation would otherwise be too small.

The reader will note that, while there is a certain uniqueness in defining the direction of the deformation in Eq. (75) to be given by the vector $`\stackrel{}{w}`$, Eq. (66), the normalization $`D=CG`$ with $`C`$ and $`G`$ given in Eqs. (80) and (83) is rather ad hoc. Within the requirements that the deformation should vanish quadratically at the collinear and soft singularities and should vanish for large $`\stackrel{}{\mathrm{}}_2`$, many other choices would be possible. The choice given here is used in the current version of the code . Surely there is some other choice that is better.

## VI The Monte Carlo integration

After the contour deformations, we have an integral of the form

$$=𝑑\mathrm{}\underset{C}{}𝒥(C;\mathrm{})g(C;\mathrm{}+i\kappa (C;\mathrm{})),$$
(84)   
where we use $`\mathrm{}`$ for the loop momenta collectively, $`\mathrm{}=\{\stackrel{}{\mathrm{}}_2,\stackrel{}{\mathrm{}}_4\}`$. The index $`C`$ labels the cut, $`a`$, $`b`$, $`c`$, or $`d`$ in Fig. 3. There is a contour deformation that depends on the cut, as specified by $`\kappa (C;\mathrm{})`$, and there is a corresponding jacobian $`𝒥(C;\mathrm{})`$, Eq. (63). Define

$$f(\mathrm{})=\mathrm{}\left\{\underset{C}{}𝒥(C;\mathrm{})g(C;\mathrm{}+i\kappa (C;\mathrm{}))\right\}.$$
(85)   
We know that $``$ is real, so

$$=𝑑\mathrm{}f(\mathrm{}).$$
(86)   
To perform the integration, we use the Monte Carlo method. We choose points $`\mathrm{}`$ with a density $`\rho (\mathrm{})`$, with

$$𝑑\mathrm{}\rho (\mathrm{})=1.$$
(87)   
After choosing $`N`$ points $`\mathrm{}_1,\mathrm{},\mathrm{}_N`$, we have an estimate for the integral:

$$_N=\frac{1}{N}\underset{i}{}\frac{f(\mathrm{}_i)}{\rho (\mathrm{}_i)}.$$
(88)   
This is an approximation for the integral in the sense that if we repeat the procedure a lot of times the expectation value for $`_N`$ is

$$_N=.$$
(89)   
The expected r.m.s. error is $``$, where

$$^2=\left(_N\right)^2=\frac{1}{N}𝑑\mathrm{}\frac{f(\mathrm{})^2}{\rho (\mathrm{})}\frac{^2}{N}.$$
(90)   
One can rewrite this as

$$^2=\frac{1}{N}𝑑\mathrm{}\rho (\mathrm{})\left(\frac{|f(\mathrm{})|}{\rho (\mathrm{})}\stackrel{~}{}\right)^2+\frac{\stackrel{~}{}^2^2}{N},$$
(91)   
where

$$\stackrel{~}{}=𝑑\mathrm{}|f(\mathrm{})|.$$
(92)   
We see, first of all, that the expected error decreases proportionally to $`1/\sqrt{N}`$. Second, we see that the ideal choice of $`\rho (\mathrm{})`$ would be $`\rho (\mathrm{})=|f(\mathrm{})|/\stackrel{~}{}`$.

Of course, it is not possible to choose $`\rho `$ in this way. But we know that $`|f|`$ has singularities at places where propagator momenta vanish and we know the structure of these singularities. We are not really able to choose $`\rho `$ so that $`|f(\mathrm{})|/\rho (\mathrm{})`$ is a constant, but at least we can choose it so that $`|f(\mathrm{})|/\rho (\mathrm{})`$ is not singular at the singularities of $`|f(\mathrm{})|`$.

Note that it is easy to combine methods for choosing Monte Carlo points. Suppose that we have a recipe for choosing points with a density $`\rho _1`$ that is singular when one propagator momentum vanishes, a recipe for choosing points with a density $`\rho _2`$ that is singular when another propagator momentum vanishes, and in general recipes for choosing points with densities $`\rho _i`$ with several goals in mind. Then we can devote a fraction $`\lambda _i`$ of the points to the choice with density $`\rho _i`$ and obtain a net density

$$\rho (\mathrm{})=\underset{i}{}\lambda _i\rho _i(\mathrm{}).$$
(93)   
### A The density near where a propagator momentum vanishes

Let $`\stackrel{}{\mathrm{}}_J`$ be the momentum of one of the propagators in our graph. We have seen that when particle $`J`$ appears in the final state, there is a factor $`1/|\stackrel{}{\mathrm{}}_J|`$ in the integrand. When propagator $`J`$ is part of a virtual loop, the contribution corresponding to this propagator being put on shell also contains a factor $`1/|\stackrel{}{\mathrm{}}_J|`$. Thus there is a singularity $`1/|\stackrel{}{\mathrm{}}_J|`$ for every propagator in the graph.

The analysis given in the introduction to this section indicates that for each propagator $`J`$ one of the terms $`\rho _i`$ in the density function should have a singularity that is at least as strong as

$$\rho _i(\mathrm{})1/|\stackrel{}{\mathrm{}}_J|$$
(94)   
as $`\stackrel{}{\mathrm{}}_J0`$. It is, of course, easy to choose points with a density proportional to $`1/|\stackrel{}{\mathrm{}}_J|^A`$ as $`\stackrel{}{\mathrm{}}_J0`$ as long as $`A<3`$. (The limitation on $`A`$ arises because for $`A3`$ we would have $`𝑑\stackrel{}{\mathrm{}}_J\rho (\mathrm{})=\mathrm{}`$.) Thus it is easy to arrange that the density of points has the requisite singularities. Specifically, we can choose $`\stackrel{}{\mathrm{}}_J`$ with the density

$$\stackrel{~}{\rho }(\stackrel{}{\mathrm{}}_J)=\frac{1}{2\pi K_0^3}\frac{1}{\left[1+\left(|\stackrel{}{\mathrm{}}_J|/K_0\right)^2\right]^2}\frac{K_0}{|\stackrel{}{\mathrm{}}_J|},$$
(95)   
where $`K_0`$ is a momentum scale determined by the other, previously chosen, loop momenta.

The singularity when $`\stackrel{}{\mathrm{}}_J0`$ can be more severe than $`1/|\stackrel{}{\mathrm{}}_J|`$, depending on the structure of the graph. Consider first the cases $`J=1,3,4,5`$. Here, the singularities for particular cuts, as given in Eq. (28), are $`1/|\stackrel{}{\mathrm{}}_J|^2`$. However, there is a cancellation after one sums over cuts (as for the singularity for $`\stackrel{}{\mathrm{}}_20`$), leaving a singularity $`1/|\stackrel{}{\mathrm{}}_J|`$.

For $`J=2`$ there is a severe singularity of the form $`1/|\stackrel{}{\mathrm{}}_J|^3`$ for particular contributions to Eq. (28). A $`1/|\stackrel{}{\mathrm{}}_2|^3`$ singularity would not be integrable, but, as we have seen in detail, there is a cancellation among the contributions so that only a $`1/|\stackrel{}{\mathrm{}}_2|^2`$ singularity is left. However, it will not do to simply chose $`\rho _i(\mathrm{})1/|\stackrel{}{\mathrm{}}_2|^2`$ because there is also a singularity in the space of the angles of $`\stackrel{}{\mathrm{}}_2`$. It is to this subject that we now turn.

### B The soft parton singularity

When two partons can scatter by exchanging a parton before they enter the final state, there is a severe singularity as the momentum of the exchanged parton goes to zero. For our graph, this happens for $`\stackrel{}{\mathrm{}}_20`$. In this subsection, we consider the behavior of the integrand for small $`\stackrel{}{\mathrm{}}_2`$ as a function of its magnitude $`\stackrel{}{\mathrm{}}_2`$ and of its direction $`\stackrel{}{u}_2=\stackrel{}{\mathrm{}}_2/|\stackrel{}{\mathrm{}}_2|`$.

The singularity for individual cuts, as given in Eq. (28), is of the form $`1/|\stackrel{}{\mathrm{}}_2|^3`$ when we let $`|\stackrel{}{\mathrm{}}_2|0`$ with $`\stackrel{}{u}_2`$ held fixed. This singularity is not integrable. However, as we have seen, the leading term cancels when we sum over cuts, leaving a $`1/|\stackrel{}{\mathrm{}}_2|^2`$ singularity for $`|\stackrel{}{\mathrm{}}_2|0`$ with $`\stackrel{}{u}_2`$ fixed.

Let us now recall from Eq. (53) that, before we deform the integration contour, the contribution for small $`\stackrel{}{\mathrm{}}_2`$ from the cut $`a`$ of Fig. 3 has, in addition to a factor $`1/|\stackrel{}{\mathrm{}}_2|^3`$, a factor $`1/[\stackrel{}{u}_2(\stackrel{}{u}_5\stackrel{}{u}_4)+iϵ]`$. That is, there is a singularity on a surface in the space of $`\stackrel{}{\mathrm{}}_2`$ whose tangent plane is the plane perpendicular to $`\stackrel{}{u}_5\stackrel{}{u}_4`$. We have avoided this singularity by deforming the integration contour. However, the deformation vanishes as $`\stackrel{}{\mathrm{}}_20`$. Thus we must face the question of what happens to the cancellation near the soft parton singularity when the contour deformation is taken into account.

First, let us recall from Eq. (75) that for cut $`a`$ in Fig. 3 the deformation has the form

$$\stackrel{}{\mathrm{}}_{2,c}=\stackrel{}{\mathrm{}}_2iD(\stackrel{}{\mathrm{}}_2)\stackrel{}{w},$$
(96)   
where $`\stackrel{}{w}=\stackrel{}{u}_1+\stackrel{}{u}_3`$. For $`\stackrel{}{\mathrm{}}_20`$,

$$\stackrel{}{w}\stackrel{}{u}_4\stackrel{}{u}_5,$$
(97)   
while $`D`$ has the form

$$D(\stackrel{}{\mathrm{}}_2)\stackrel{}{\mathrm{}}_2^{\mathrm{\hspace{0.17em}2}}\stackrel{~}{D}(\stackrel{}{u}_2).$$
(98)   
Here $`\stackrel{~}{D}`$ vanishes for $`\stackrel{}{u}_2=\stackrel{}{u}_4`$ and for $`\stackrel{}{u}_2=\stackrel{}{u}_5`$ and is positive elsewhere. Thus

$$\stackrel{}{\mathrm{}}_{2,c}=|\stackrel{}{\mathrm{}}_2|\left(\stackrel{}{u}_2+i|\stackrel{}{\mathrm{}}_2|\stackrel{~}{D}(\stackrel{}{u}_2)(\stackrel{}{u}_5\stackrel{}{u}_4)+𝒪(\stackrel{}{\mathrm{}}_2^{\mathrm{\hspace{0.17em}2}})\right).$$
(99)   
Substituting $`\mathrm{}_{2,c}`$ as given above for $`\stackrel{}{\mathrm{}}_2`$ in Eq. (53), we obtain an expression for the contribution from cut $`a`$ to the integrand on the deformed contour near the soft singularity:

$$𝒢_a\frac{𝒮_2}{32|\stackrel{}{\mathrm{}}_4|^2|\stackrel{}{\mathrm{}}_5|^2}\frac{1}{|\stackrel{}{\mathrm{}}_2|^3}\frac{1}{1+\stackrel{}{u}_2\stackrel{}{u}_4}\frac{1}{1\stackrel{}{u}_2\stackrel{}{u}_5}\frac{2\stackrel{}{u}_2(\stackrel{}{u}_5\stackrel{}{u}_4)}{\stackrel{}{u}_2(\stackrel{}{u}_5\stackrel{}{u}_4)+i|\stackrel{}{\mathrm{}}_2|\stackrel{~}{D}(\stackrel{}{u}_2)(\stackrel{}{u}_5\stackrel{}{u}_4)^2}.$$
(100)   
There are two cases to consider. First, when $`|\stackrel{}{\mathrm{}}_2|0`$ with $`\stackrel{}{u}_2`$ fixed, we can drop the second term in the last denominator. Then $`𝒢_ah(\stackrel{}{u}_2)/|\stackrel{}{\mathrm{}}_2|^3`$, where the function $`h(\stackrel{}{u}_2)`$ is the same as on the undeformed contour. As we have seen, the leading $`1/|\stackrel{}{\mathrm{}}_2|^3`$ terms cancel when one sums over cuts. Thus, as noted earlier, the net integrand behaves like

$$𝒢h_{\mathrm{tot}}(\stackrel{}{u}_2)/|\stackrel{}{\mathrm{}}_2|^2$$
(101)   
when $`|\stackrel{}{\mathrm{}}_2|0`$ with $`\stackrel{}{u}_2`$ fixed.

The second case is more interesting. Consider $`|\stackrel{}{\mathrm{}}_2|0`$ and $`\stackrel{}{u}_2(\stackrel{}{u}_5\stackrel{}{u}_4)0`$ with $`\stackrel{}{u}_2(\stackrel{}{u}_5\stackrel{}{u}_4)/|\stackrel{}{\mathrm{}}_2|`$ fixed. Then $`𝒢_a`$ is more singular, $`𝒢_a1/|\stackrel{}{\mathrm{}}_2|^4`$. To see what happens in this region, we analyze the contribution from cut $`b`$ in Fig. 3 in the same fashion. The contour deformation for cut $`b`$ is different from that for cut $`a`$, but the deformations match at leading order as $`|\stackrel{}{\mathrm{}}_2|0`$. (This is an important feature of the choice of contour deformations.) Thus we can use Eq. (99) in Eq. (54) to obtain

$$𝒢_b\frac{𝒮_2}{32|\stackrel{}{\mathrm{}}_4|^2|\stackrel{}{\mathrm{}}_5|^2}\frac{1}{|\stackrel{}{\mathrm{}}_2|^3}\frac{1}{1\stackrel{}{u}_2\stackrel{}{u}_4}\frac{1}{1+\stackrel{}{u}_2\stackrel{}{u}_5}\frac{2+\stackrel{}{u}_2(\stackrel{}{u}_5\stackrel{}{u}_4)}{\stackrel{}{u}_2(\stackrel{}{u}_5\stackrel{}{u}_4)+i|\stackrel{}{\mathrm{}}_2|\stackrel{~}{D}(\stackrel{}{u}_2)(\stackrel{}{u}_5\stackrel{}{u}_4)^2}.$$
(102)   
We see that $`𝒢_b`$ is also proportional to $`1/|\stackrel{}{\mathrm{}}_2|^4`$ in the problematic region. However, since $`u_2u_4u_2u_5`$ in this region, the leading $`1/|\stackrel{}{\mathrm{}}_2|^4`$ behavior cancels when we add $`𝒢_b`$ to $`𝒢_a`$. We are left with the next term, proportional to $`1/|\stackrel{}{\mathrm{}}_2|^3`$.

For the two remaining cuts there is no contour deformation. The contributions from these cuts are each proportional to $`1/|\stackrel{}{\mathrm{}}_2|^3`$. Calculation shows that there is no further cancellation. Thus the net behavior of the integrand is

$$𝒢1/|\stackrel{}{\mathrm{}}_2|^3$$
(103)   
when $`|\stackrel{}{\mathrm{}}_2|0`$ and $`\stackrel{}{u}_2(\stackrel{}{u}_5\stackrel{}{u}_4)0`$ with $`\stackrel{}{u}_2(\stackrel{}{u}_5\stackrel{}{u}_4)/|\stackrel{}{\mathrm{}}_2|`$ fixed.

### C Density near a soft parton singularity

According to the analysis at the beginning of this section, we should choose a density of integration points that has a singularity that is at least as strong as that of $`|𝒢|`$ near the soft singularity at $`\stackrel{}{\mathrm{}}_20`$. Thus we should choose one of the $`\rho _i`$ so that

$`\rho _i(\mathrm{})`$ $``$ $`{\displaystyle \frac{1}{|\stackrel{}{\mathrm{}}_2|^p}},|\stackrel{}{\mathrm{}}_2|0,\stackrel{}{u}_2\mathrm{fixed},`$ (104)   
$`\rho _i(\mathrm{})`$ $``$ $`{\displaystyle \frac{1}{|\stackrel{}{\mathrm{}}_2|^{p+1}}},|\stackrel{}{\mathrm{}}_2|0,{\displaystyle \frac{\stackrel{}{u}_2(\stackrel{}{u}_5\stackrel{}{u}_4)}{|\stackrel{}{\mathrm{}}_2|}}\mathrm{fixed},`$ (105)   
with $`p2`$.

Specifically, having chosen $`\stackrel{}{\mathrm{}}_4`$ we can choose the remaining loop momentum $`\stackrel{}{\mathrm{}}_2`$ with the density

$$\stackrel{~}{\rho }(\stackrel{}{\mathrm{}}_2)=\frac{1}{2\pi K_0^3}\frac{1}{\left[1+\left(|\stackrel{}{\mathrm{}}_2|/K_0\right)^{(3p)}\right]^{(5p)/(3p)}}\left(\frac{K_0}{|\stackrel{}{\mathrm{}}_2|}\right)^p\frac{1}{\mathrm{\Gamma }\sqrt{\mathrm{cos}^2(\theta )+\stackrel{}{\mathrm{}}_2^{\mathrm{\hspace{0.17em}2}}/K_0^2}}.$$
(106)   
Here $`K_0`$ is a momentum scale determined by $`\stackrel{}{\mathrm{}}_4`$, $`\theta `$ is the angle between $`\stackrel{}{\mathrm{}}_2`$ and $`(\stackrel{}{u}_5\stackrel{}{u}_4)`$, and

$$\mathrm{sinh}(\mathrm{\Gamma })=K_0/|\stackrel{}{\mathrm{}}_2|.$$
(107)   
It is easy to choose points with this density by first choosing $`|\stackrel{}{\mathrm{}}_2|`$, then choosing $`\mathrm{cos}(\theta )`$, and finally choosing the corresponding azimuthal angle $`\varphi `$ with a uniform density. Accounting for the fact that $`\mathrm{\Gamma }\mathrm{log}(|\stackrel{}{\mathrm{}}_2|)`$ for $`\stackrel{}{\mathrm{}}_20`$, we see that $`\stackrel{~}{\rho }`$ will have a singularity stronger than that of $`𝒢`$ provided that $`p>2`$. We will see how this works in a numerical example in the next section.

## VII Numerical Example

In this section, I illustrate the principles developed above by means of a particular example. We consider the integral in Eq. (19). We hold $`\stackrel{}{\mathrm{}}_4`$ fixed and consider the integrand as a function of $`\stackrel{}{\mathrm{}}_2`$. In order to simplify the labelling, I define

$$\stackrel{}{\mathrm{}}_2\stackrel{}{\mathrm{}}.$$
(108)   
There is a contribution for each cut $`C`$, with $`C=a,b,c,\mathrm{or}d`$. For each contribution from a cut $`C`$ in which there is a virtual loop, we want to deform the integration contour as discussed in Sec. V. Thus $`\stackrel{}{\mathrm{}}`$ gets replaced by a complex vector $`\stackrel{}{\mathrm{}}_c=\stackrel{}{\mathrm{}}+i\stackrel{}{\kappa }_C`$ and we need to supply a jacobian $`𝒥_C(\stackrel{}{\mathrm{}})`$, Eq. (63). Then the integration over $`\stackrel{}{\mathrm{}}`$ has the form

$$𝑑\stackrel{}{\mathrm{}}\underset{C}{}𝒥_C(\stackrel{}{\mathrm{}})𝒢_C(\stackrel{}{\mathrm{}}).$$
(109)   
The functions $`𝒢_C`$ are the analytic continuations to the deformed contours of the functions given in Eq. (28). As discussed in Sec. VI, the quantity that is relevant for the convergence of the Monte Carlo integration is the integrand divided by the density of points chosen for the integration. In this section, I consider only the integration over $`\stackrel{}{\mathrm{}}`$, so I discuss a choice for the density of integration points $`\rho (\stackrel{}{\mathrm{}})`$ at a fixed $`\stackrel{}{\mathrm{}}_4`$ and display plots of the functions

$$F_C(\stackrel{}{\mathrm{}})\frac{1}{\rho (\stackrel{}{\mathrm{}})}𝒥_C(\stackrel{}{\mathrm{}})𝒢_C(\stackrel{}{\mathrm{}})$$
(110)   
and $`F(\stackrel{}{\mathrm{}})=_CF_C(\stackrel{}{\mathrm{}})`$, as well as plots of the deformation and the density.

For the numerical examples, I choose

$$\stackrel{}{q}=(3,0.5,0)$$
(111)   
and then take $`\stackrel{}{\mathrm{}}_4`$ at the point

$$\stackrel{}{\mathrm{}}_4=(2,1,0).$$
(112)   
Since $`\stackrel{}{\mathrm{}}_5=\stackrel{}{q}\stackrel{}{\mathrm{}}_4`$ we have

$$\stackrel{}{\mathrm{}}_5=(1,0.5,0).$$
(113)   
The singularities of the functions $`𝒢_C(\stackrel{}{\mathrm{}})`$ lie in the plane of $`\stackrel{}{\mathrm{}}_4`$ and $`\stackrel{}{\mathrm{}}_5`$, that is the $`\mathrm{}_z=0`$ plane. In the plots, I choose $`\mathrm{}_z=0`$, so that we see the effect of the singularities. I plot $`|\stackrel{}{\kappa }_a|`$, $`|\stackrel{}{\kappa }_b|`$, $`\rho `$, $`F_a`$, $`F_b`$, $`F_c+F_d`$ and $`F`$ as functions of $`\mathrm{}_x`$ and $`\mathrm{}_y`$ in the domain $`2.5<\mathrm{}_x<1.0`$ and $`1.0<\mathrm{}_y<2.0`$.

Consider first the contour deformation for cut $`a`$, $`\stackrel{}{\mathrm{}}\stackrel{}{\mathrm{}}_c=\stackrel{}{\mathrm{}}+i\stackrel{}{\kappa }_a`$. I take $`\stackrel{}{\kappa }_a=D\stackrel{}{w}`$ as given in Eqs. (75-83) with $`\alpha =\beta =\gamma =1`$. In Fig. 7, I show a graph of $`|\stackrel{}{\kappa }_a|`$ versus $`\mathrm{}_x`$ and $`\mathrm{}_y`$. We see that the deformation is not small. I also display in the figure the lines $`\stackrel{}{\mathrm{}}=x\stackrel{}{\mathrm{}}_4`$ with $`0<x<1`$ and $`\stackrel{}{\mathrm{}}=x\stackrel{}{\mathrm{}}_5`$ with $`0<x<1`$, where the collinear singularities for cut $`a`$ are located. We see that, as desired, the deformation vanishes quadratically as $`\stackrel{}{\mathrm{}}`$ approaches these lines.

There is a different contour deformation for cut $`b`$. The same formulas apply as for cut $`a`$ with the replacements $`\stackrel{}{\mathrm{}}_4\stackrel{}{\mathrm{}}_1`$, $`\stackrel{}{\mathrm{}}_5\stackrel{}{\mathrm{}}_3`$, $`\stackrel{}{\mathrm{}}\stackrel{}{\mathrm{}}`$ and with the sign of $`\stackrel{}{\kappa }`$ reversed. I show a graph of $`|\stackrel{}{\kappa }_b|`$ versus $`\mathrm{}_x`$ and $`\mathrm{}_y`$ in Fig. 8. (This figure does not look like Fig. 7 because $`\stackrel{}{\mathrm{}}`$ varies with $`\stackrel{}{\mathrm{}}_4`$ held fixed, not with $`\stackrel{}{\mathrm{}}_1`$ held fixed as would be needed if we applied the replacement $`\stackrel{}{\mathrm{}}_4\stackrel{}{\mathrm{}}_1`$ to Fig. 7.) I also display in the figure the lines $`\stackrel{}{\mathrm{}}=\lambda \stackrel{}{\mathrm{}}_4`$ with $`0<\lambda `$ and $`\stackrel{}{\mathrm{}}=\lambda \stackrel{}{\mathrm{}}_5`$ with $`0<\lambda `$, where the collinear singularities for cut $`b`$ are located. The deformation vanishes quadratically as $`\stackrel{}{\mathrm{}}`$ approaches these lines.

The jacobian functions $`𝒥_a(\stackrel{}{\mathrm{}})`$ and $`𝒥_b(\stackrel{}{\mathrm{}})`$ associated with the contour deformations are quite unremarkable, so I omit showing them.

Consider next the density of integration points. I choose

$$\rho (\stackrel{}{\mathrm{}})=0.2\rho _1(\stackrel{}{\mathrm{}}_1)+0.6\rho _2(\stackrel{}{\mathrm{}})+0.2\rho _3(\stackrel{}{\mathrm{}}_3),$$
(114)   
as shown in Fig. 9. The function $`\rho _1`$ has a mild singularity as $`\stackrel{}{\mathrm{}}_10`$ and is given by Eq. (95) with $`\stackrel{}{\mathrm{}}_1=\stackrel{}{\mathrm{}}_4\stackrel{}{\mathrm{}}`$ and with $`K_0`$ set equal to 2. The function $`\rho _3`$ has a mild singularity as $`\stackrel{}{\mathrm{}}_30`$; I use the same functional form with $`\stackrel{}{\mathrm{}}_3=\stackrel{}{\mathrm{}}\stackrel{}{\mathrm{}}_5`$. For $`\rho _2`$, I use the function given in Eq. (106) with $`K_0=2`$ and with the power $`p`$ taken as $`p=2.2`$. Then $`\rho _2`$ has a strong $`1/[|\stackrel{}{\mathrm{}}|^{2.2}\mathrm{log}(|\stackrel{}{\mathrm{}}|)]`$ singularity as we approach the $`\stackrel{}{\mathrm{}}=0`$. Furthermore, the density of points is largest near the plane $`\mathrm{}_y=0`$, the plane that is tangent at $`\stackrel{}{\mathrm{}}=0`$ to the ellipsoidal surface that (if we turn off the deformation) contains the scattering singularity. In order to display the dependence of $`\rho `$ on angle near $`\stackrel{}{\mathrm{}}=0`$, I plot in Fig. 10 the angle dependent factor in $`\rho _2`$, namely the factor

$$\frac{|\stackrel{}{\mathrm{}}|/K_0}{\sqrt{\mathrm{cos}^2(\theta )+\stackrel{}{\mathrm{}}^{\mathrm{\hspace{0.17em}2}}/K_0^2}}$$
(115)   
in Eq. (106), in a region near $`\stackrel{}{\mathrm{}}=0`$. Here $`\mathrm{cos}(\theta )=\mathrm{}_y/|\stackrel{}{\mathrm{}}|`$. We see that the density of integration points is heavily concentrated very near the plane $`\mathrm{}_y=0`$ when $`|\stackrel{}{\mathrm{}}|`$ is small.

We are now ready to look at the contribution $`F_a=𝒥_a𝒢_a/\rho `$ to $`F`$ from cut $`a`$. This function is displayed in Fig. 11 with a small rectangle near $`\stackrel{}{\mathrm{}}=0`$ removed from the graph. We see the two collinear singularities, at $`\stackrel{}{\mathrm{}}=x\stackrel{}{\mathrm{}}_4`$ and at $`\stackrel{}{\mathrm{}}=x\stackrel{}{\mathrm{}}_5`$ with $`0<x<1`$. As $`\stackrel{}{\mathrm{}}`$ approaches one of these singularities, $`F_a`$ approaches $`\mathrm{}`$.

In the standard method for calculating $``$, we would perform the integration over $`\stackrel{}{\mathrm{}}`$ analytically for the contribution from cut $`a`$. Because of the singularities, the integration is divergent. However, we can get a finite answer if we regulate the integral by working in $`32ϵ`$ spatial dimensions. Then the result contains terms proportional to $`1/ϵ^2`$ and $`1/ϵ`$ as well as a remainder that is finite as $`ϵ0`$.

What about the contribution to $`F`$ from cut $`b`$, the other cut for which there is a virtual subgraph? This function is displayed in Fig. 12 with the same small rectangle near $`\stackrel{}{\mathrm{}}=0`$ removed from the graph. We see the two collinear singularities, at $`\stackrel{}{\mathrm{}}=\lambda \stackrel{}{\mathrm{}}_4`$ and at $`\stackrel{}{\mathrm{}}=\lambda \stackrel{}{\mathrm{}}_5`$ with $`0<\lambda <\mathrm{}`$. As with $`F_a`$, as $`\stackrel{}{\mathrm{}}`$ approaches one of these singularities, $`F_b`$ approaches $`\mathrm{}`$.

There are two cuts, $`c`$ and $`d`$, for which there are no virtual subgraphs. In Fig. 13 I show the contribution $`F_c+F_d`$ from these cuts. We see that $`F_c+F_d`$ approaches $`+\mathrm{}`$ at just the singularities where $`F_a`$ and $`F_b`$ approach $`\mathrm{}`$.

In the standard method for QCD calculations, we would perform the integration over $`\stackrel{}{\mathrm{}}`$ partially numerically for the contribution from cuts $`c`$ and $`d`$. Of course, we would have to do something about the collinear and soft singularities, since otherwise we would obtain an infinite result. For instance, if we were to use the phase-space slicing method, we would slice away a small part of the integration domain near the singularities and calculate its contribution analytically in $`32ϵ`$ spatial dimensions in the limit that the region sliced away is small. Then we would be left with a numerical integration of $`𝒢_c+𝒢_d`$ in the remaining region (in exactly 3 spatial dimensions). Evidently, the density of points used in the present method would not do for this purpose; we would need to expend more points on the region near the collinear singularities.

We see that the standard method for performing the integrations, in which some parts of the integrations are performed analytically and some are performed numerically, is, of necessity, rather complicated. In the numerical method, we simply combine $`𝒢_a`$, $`𝒢_b`$, $`𝒢_c`$, and $`𝒢_d`$ and integrate numerically. The argument in the preceding sections showed that the contributions from the various cuts cancel as one approaches the collinear singularities. This is illustrated in Fig. 14, where I plot $`F_a+F_b+F_c+F_d`$ versus $`\mathrm{}_x`$ and $`\mathrm{}_y`$. We see, first of all, that the singular behaviors at the collinear singularities cancel, just as the calculation of Sec. IV showed. There is also a cancellation at the soft singularity at $`\stackrel{}{\mathrm{}}=0`$. There is still a singularity in the integrand at $`\stackrel{}{\mathrm{}}=0`$, but it is integrable and is removed from $`F`$ by choosing a suitable density of points $`\rho `$. Thus $`F`$ remains less than about 20 everywhere.

We can see the remnants of the scattering singularity, which is located on an ellipsoidal surface that intersects the plane $`\mathrm{}_z=0`$. If it were not for the contour deformation, $`F`$ would be singular on this surface, approaching $`+\mathrm{}`$ as one approached the surface from one side and approaching $`\mathrm{}`$ as one approached the surface from the other side. Since the deformed contour avoids the singularity, the singular behavior is removed and we are left with a ridge and valley near the ellipsoid. This structure is illustrated in Fig. 15, in which we see a slice through the ridge and valley at $`\mathrm{}_x=0.3`$. Since the amount of deformation vanishes as one approaches $`\stackrel{}{\mathrm{}}=0`$, the width of the ridge and valley structure becomes more and more narrow as $`|\stackrel{}{\mathrm{}}|0`$. Recall that the density $`\rho `$ of integration points is designed to match this increasing narrowness as $`|\stackrel{}{\mathrm{}}|0`$, so that the integration points are concentrated where the structure is.

## VIII Other issues

In the preceding sections, we have seen the most important features of the method of numerical integration for one loop QCD calculations. There are other important issues that are outside of the scope of this paper. I mention these briefly here.

First, full QCD has a much more complicated structure than $`\varphi ^3`$ theory, which was the example for this paper. However, the complications of full QCD are in the numerators of the expressions representing Feynman diagrams, while the cancellations and the analytic structure related to the contour deformation have to do with the denominator structure. Thus one can simply generate the numerator structure with computer algebra and carry it along.

Second, the denominator structure in the example used in this paper is not the only denominator structure that one needs to treat. In the QCD calculation for three-jet-like quantities in electron-positron annihilation at order $`\alpha _s^2`$, there are five possibilities for how a virtual subgraph can occur inside an amplitude. The possibilities are indicated in Fig. 16. For each possibility there is an entering line representing the virtual photon or $`Z`$ boson, which we take to have zero three momentum, and there are three on-shell lines entering the final state. There are graphs of two types, (a) and (b), containing two point virtual subgraphs. There are graphs of two types, (c) and (d), containing three point virtual subgraphs. There is one type, (e), of graph with a four point virtual subgraph. In structure (c), a line with non-zero three momentum enters the three point virtual subgraph and two on-shell lines leave. This is the case that we analyzed in the example of this paper. The structure of the graph led to the singularity structure depicted in Fig. 6. Amplitudes of types (d) and (e) have different singularity structures from that studied here. Case (d) is simpler than the case we have studied, while case (e) is somewhat more complicated. However, the essential features are those that we have already studied.

This leaves virtual self-energy subgraphs. In case (a), there is a self energy subgraph on a propagator that enters the final state. This case requires a treatment different from that discussed in the previous sections. This is evident because there is a nominal $`1/k^2`$ where $`k^2=0`$. The treatment required is to represent the virtual self-energy via a dispersion relation . In this representation, the subgraph is expressed as an integral over the three-momentum in the virtual loop with an integrand that is closely related to the integrand for the corresponding cut self-energy graph. The point-by-point cancellation between real and virtual graphs is then manifest. It is convenient also to use the dispersive representation for the much easier case (b).

Third, we have to do something about ultraviolet divergences in virtual subgraphs. These are easily removed  by subtracting an integrand that, in the region of large loop momenta, matches the integrand of the divergent subdiagram. The integrand of the subtraction term should depend on a mass parameter $`\mu `$ that serves to make the subtraction term well behaved in the infrared. Then, with the aid of a small analytical calculation for each of the one loop divergent subdiagrams that occur in QCD, one can arrange the definition so that this ad hoc subtraction has exactly the same effect as $`\overline{\mathrm{MS}}`$ subtraction with scale parameter $`\mu `$.

Fourth, I use Feynman gauge for the gluon field, but then self-energy corrections on a gluon propagator require special attention . The one loop gluon self-energy subgraph, $`\pi ^{\mu \nu }(k)`$, contains a term proportional to $`k^\mu k^\nu `$ that contributes quadratic infrared divergences , while the cancellation mechanism that we have studied in this paper takes care of logarithmic divergences. This problem can be solved by replacing $`\pi ^{\mu \nu }`$ by $`P_\alpha ^\mu \pi ^{\alpha \beta }P_\beta ^\nu `$, where $`P_\alpha ^\mu =g_\alpha ^\mu k^\mu \stackrel{~}{k}_\alpha /\stackrel{~}{k}^2`$, with $`\stackrel{~}{k}=(0,\stackrel{}{k})`$. The terms added to $`\pi ^{\mu \nu }`$ are proportional to either $`k^\mu `$ or $`k^\nu `$ and thus vanish when one sums over different ways of inserting the dressed gluon propagator into the remaining subgraph. Since $`P_\alpha ^\mu k^\alpha =0`$, the problematic $`k^\mu k^\nu `$ term is eliminated. Effectively, this is a change of gauge for dressed gluon propagators from Feynman gauge to Coulomb gauge.

Most of these issues are discussed briefly in . A quite detailed, but not very pedagogical, treatment can be found in . Further analysis of these issues is left for a future paper.

I have also not given a complete presentation of algorithms for choosing integration points. As discussed in Sec. VI C, the crucial issue is to have the right singularities in the density of points near a soft parton singularity of the Feynman diagram. This is not the only issue that needs to be addressed in a complete algorithm. Of course, the demonstration program  has a complete algorithm. However, this algorithm is quite a hodge-podge of methods and it seems that a detailed exposition should be reserved for a better and more systematic method, which remains to be developed.

## IX Results and conclusions

In the preceding sections, we have seen some of the techniques needed for the numerical integration method for QCD calculations. Of course, since not all of the techniques have been explained, the explanation does not constitute a very convincing argument that such a calculation is feasible. A truly exhaustive explanation would help, but an actual computer program that demonstrates the techniques is better. Results from such a program were presented in . Since then, I have found and corrected one bug that resulted in errors a little bit bigger than 1% and have made some other improvements in the code. The resulting code and documentation are available at .

The program is a parton-level event generator. The user is to supply a subroutine that calculates how an event with three or four partons in the final state contributes to the observable to be calculated. The program supplies events, each consisting of a set of parton momenta $`\{\stackrel{}{k}_1,\stackrel{}{k}_2,\stackrel{}{k}_3\}`$ or $`\{\stackrel{}{k}_1,\stackrel{}{k}_2,\stackrel{}{k}_3,\stackrel{}{k}_4\}`$, together with weights $`w`$ for the events. Then the user routine calculates $``$ according to

$$\frac{1}{N}\underset{i=1}{\overset{N}{}}w_i𝒮(k_i).$$
(116)   
The weights used are the real parts of complex weights; the imaginary parts can be dropped since we always know in advance that $``$ is real. Thus the weights are both positive and negative. It would, of course, be more convenient to have only positive weights, but one can hardly have quantum interference without having negative numbers along with positive numbers.

The first general purpose program for QCD calculation of three-jet-like quantities in $`e^+e^{}\mathrm{ℎ𝑎𝑑𝑟𝑜𝑛𝑠}`$ at order $`\alpha _s^2`$ was that of Kunszt and Nason . This program uses the numerical/analytical method of Ellis, Ross, and Terrano. In Table I, I compare the results of the Kunszt-Nason program to those obtained with the numerical method for the $`\alpha _s^2`$ contributions to moments of the thrust distribution,

$$_n=\frac{1}{\sigma _0(\alpha _s/\pi )^2}_0^1𝑑t(1t)^n\frac{d\sigma ^{[2]}}{dt}.$$
(117)   
In Table II, I compare the results of the two methods for moments of the $`y_{\mathrm{cut}}`$ distribution for the three jet cross section. To define these quantities, let $`f_3(y_{\mathrm{cut}})`$ be the cross section to produce three jets according to the Durham algorithm  with resolution parameter $`y_{\mathrm{cut}}`$. Let $`g_3(y_{\mathrm{cut}})`$ be the negative of its derivative,

$$g_3(y_{\mathrm{cut}})=\frac{f_3(y_{\mathrm{cut}})}{dy_{\mathrm{cut}}}.$$
(118)   
Then we calculate moments of the $`\alpha _s^2`$ contribution to this quantity,

$$_n=\frac{1}{\sigma _0(\alpha _s/\pi )^2}_0^1𝑑y_{\mathrm{cut}}(y_{\mathrm{cut}})^ng_3^{[2]}.$$
(119)   
In each table, the results for the numerical method are shown with their statistical and systematic errors. (The systematic error is estimated by changing the cutoffs that remove small regions near the singularities where roundoff errors start to become a problem.) The corresponding results for the numerical/analytical method are shown with the statistical errors as reported by the program. Inspection of the tables shows that there is good agreement between the two methods.

We have explored some of the most important techniques necessary for a QCD calculation for three-jet-like quantities in electron-positron annihilation at order $`\alpha _s^2`$ using numerical integration throughout the calculation. For the techniques covered, this explanation expands on the brief presentation in . We have also seen that the method works. The older and very successful numerical/analytical method for QCD calculations has its complications. The numerical method has its own complications, but they are different complications. Thus one may expect that the classes of problems for which each of the methods is well adapted may be different. There may be some classes of problems for which the natural flexibility of the numerical method makes it the more useful method. It remains for the future to explore the possibilities.

###### Acknowledgements.

I thank M. Seymour and Z. Kunszt for providing advice and results from their programs to help with debugging the program described here. I thank P. Nason for providing the Kunszt-Nason program that I used in preparing the Tables I and II. I thank M. Krämer for criticisms of the manuscript. Finally, I am most grateful to the TH division of CERN for its hospitality during a sabbatical year in which much of the writing of this paper was accomplished.    

## Contour deformation in many dimensions

The calculational method described in this paper makes use of Cauchy’s theorem in a multi-dimensional complex space. Since this theorem is not proved in most textbooks on complex analysis, I provide a proof here, including the special case, needed for our application, in which there is a singularity on the integration contour.

Let $`f(z)`$ be a function of $`N`$ complex variables $`z^\mu =x^\mu +iy^\mu `$, with $`\mu =1,\mathrm{},N`$, where $`x^\mu `$ and $`y^\mu `$ are real variables. Consider a family of integration contours $`𝒞(t)`$ labeled by a parameter $`t`$ with $`0t1`$ and specified by

$$z^\mu (x^1,\mathrm{},x^N;t)=x^\mu +iy^\mu (x^1,\mathrm{},x^N;t),\mu =1,\mathrm{},N.$$
(120)   
Let $`(t)`$ be the integral of $`f`$ over the contour $`𝒞(t)`$,

$$(t)=_{𝒞(t)}𝑑zf(z)=𝑑xdet\left(\frac{z(x;t)}{x}\right)f(z(x;t)).$$
(121)   
Suppose that $`f(z)`$ is analytic in a region that contains the contours. Then we have Cauchy’s theorem:

$$(1)=(0).$$
(122)   
To prove this theorem, we simply prove that $`d(t)/dt=0`$. Define

$$A_\nu ^\mu =\frac{z^\mu }{x^\nu }=\delta _\nu ^\mu +i\frac{y^\mu }{x^\nu }.$$
(123)   
Let $`B_\nu ^\mu `$ be the inverse matrix to $`A_\nu ^\mu `$. Then

$$B_\nu ^\mu detA=\frac{1}{(N1)!}ϵ^{\mu \mu _2\mathrm{}\mu _N}ϵ_{\nu \nu _2\mathrm{}\nu _N}\frac{z^{\nu _2}}{x^{\mu _2}}\mathrm{}\frac{z^{\nu _N}}{x^{\mu _N}},$$
(124)   
where $`ϵ^{\mu _1\mathrm{}\mu _N}`$ is the completely antisymmetric tensor with $`N`$ indices, normalized to $`ϵ^{\mathrm{1\hspace{0.17em}2}\mathrm{}N}=1`$, and $`ϵ_{\mu _1\mathrm{}\mu _N}`$ is the same tensor. This has the immediate consequence that

$$\frac{}{x^\mu }\left(B_\nu ^\mu detA\right)=0.$$
(125)   
Also,

$$detA=\frac{1}{N!}ϵ^{\mu _1\mathrm{}\mu _N}ϵ_{\nu _1\mathrm{}\nu _N}\frac{z^{\nu _1}}{x^{\mu _1}}\mathrm{}\frac{z^{\nu _N}}{x^{\mu _N}},$$
(126)   
so

$$\frac{}{t}detA=\frac{A_\mu ^\nu }{t}B_\nu ^\mu detA=i\frac{y^\nu }{x^\mu t}B_\nu ^\mu detA.$$
(127)   
We need one more result:

$$\frac{f}{x^\mu }=\frac{z^\nu }{x^\mu }\frac{f}{z^\nu },$$
(128)   
so

$$\frac{f}{z^\nu }=B_\nu ^\mu \frac{f}{x^\mu }.$$
(129)   
Then, using the results (125), (127), and (129) and an integration by parts, we find

$`{\displaystyle \frac{d}{dt}}(t)`$ $`=`$ $`{\displaystyle 𝑑x\frac{d}{dt}\left[detAf\right]}`$ (130)   
$`=`$ $`{\displaystyle 𝑑xdetA\left\{i\frac{y^\nu }{x^\mu t}B_\nu ^\mu f+\frac{f}{z^\nu }i\frac{y^\nu }{t}\right\}}`$ (131)   
$`=`$ $`{\displaystyle 𝑑xdetA\left\{i\frac{y^\nu }{x^\mu t}B_\nu ^\mu f+iB_\nu ^\mu \frac{f}{x^\mu }\frac{y^\nu }{t}\right\}}`$ (132)   
$`=`$ $`i{\displaystyle 𝑑xB_\nu ^\mu detA\frac{}{x^\mu }\left\{\frac{y^\nu }{t}f\right\}}`$ (133)   
$`=`$ $`i{\displaystyle 𝑑x\frac{y^\nu }{t}f\frac{}{x^\mu }\left\{B_\nu ^\mu detA\right\}}`$ (134)   
$`=`$ $`0.`$ (135)   
This proves the theorem.

Consider now a more complicated problem. Suppose that we have an integral of the form

$$=𝑑x[f(x)+g(x)].$$
(136)   
where $`f`$ and $`g`$ are both singular on a surface $`𝒫`$ in the space of the real variables $`x`$. Suppose that the strength of the singularities are such that the integral of either function would be logarithmically divergent. Suppose further that there is a cancellation in the sum such that the integral of $`f+g`$ is convergent. Let $`d(x)`$ be the distance from any point $`x`$ to the surface $`𝒫`$. Let us cut out a region of radius $`R`$ around $`𝒫`$ and write

$$=\underset{R0}{lim}\left[_{d>R}𝑑xf(x)+_{d>R}𝑑xg(x)\right].$$
(137)   
Now we wish to explore the consequences of deforming the integration contour for the integral of $`f`$. Thus we investigate (with the same notation as above)

$$(t)=\underset{R0}{lim}\left[_{d>R}𝑑xdet\left(\frac{z(x;t)}{x}\right)f(z(x;t))+_{d>R}𝑑xg(x)\right].$$
(138)   
Following the previous proof we find that there is a surface term in the integration by parts

$`{\displaystyle \frac{d}{dt}}(t)`$ $`=`$ $`\underset{R0}{lim}\left[i{\displaystyle _{d<R}}𝑑x{\displaystyle \frac{}{x^\mu }}\left\{B_\nu ^\mu detA{\displaystyle \frac{y^\nu }{t}}f\right\}\right]`$ (139)   
$`=`$ $`\underset{R0}{lim}\left[i{\displaystyle 𝑑S_\mu \left\{B_\nu ^\mu detA\frac{y^\nu }{t}f\right\}}\right],`$ (140)   
where the integration is over the surface $`d=R`$ and $`dS_\mu `$ is the surface area differential normal to the surface.

We want to arrange the deformation specified by $`y^\mu (x;t)`$ so that $`d(t)/dt=0`$. For this to happen, it is clear that $`y`$ will have to approach 0 as $`x`$ approaches the surface $`𝒫`$. Then $`B_\nu ^\mu \delta _\nu ^\mu `$ and $`detA1`$ as $`x`$ approaches $`𝒫`$. Let the dimensionality of the singular surface $`𝒫`$ be $`Na`$. If the function $`f`$ was such that the original integral was logarithmically divergent, then $`fR^a`$ for $`R0`$. The integration over the surface gives a factor $`dS^\mu R^{a1}`$ for $`R0`$. Suppose that the deformation vanishes proportionally to $`R^b`$. Then

$$\frac{d}{dt}(t)\underset{R0}{lim}\left[R^{a1}R^aR^b\right].$$
(141)   
Then $`d(t)/dt=0`$ if $`b>1`$. The choice made in the main text of the paper is $`b=2`$.