{'fraction_non_alphanumeric': 0.03699319325244155, 'fraction_numerical': 0.009132033991459856, 'mean_word_length': 4.029449287688709, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Next-to-leading order QCD calculations with parton showers II: soft singularities

Programs that calculate observables in quantum chromodynamics at next-to-leading order typically generate events that consist of partons rather than hadrons -- and just a few partons at that. These programs would be much more useful if the few partons were turned into parton showers, which could be given to one of the Monte Carlo event generators to produce hadron showers. In a previous paper, we have seen how to generate parton showers related to the final state collinear singularities of the perturbative calculation for the example of e+ + e- -->3 jets. This paper discusses the treatment of the soft singularities.


I. INTRODUCTION
This is the second of two papers concerning how to modify a next-to-leading order (NLO) calculation in quantum chromodynamics by adding parton showers in such a way that when the generated events are used to calculate an infrared safe observable, the observable is calculated correctly at next-to-leading order. In the companion paper [1], Krämer and the present author explain the aspects of this problem that relate to the divergences that appear in perturbation theory when two massless partons become collinear. This part of the method is quite simple and easy to understand. The treatment of the divergences that appear when a gluon becomes soft is a bit more technical and will be given in this paper.
We find that it it is conceptually simplest to match showers to a NLO calculation if the NLO calculation is performed in the Coulomb gauge. That is because the collinear divergences of the theory in a physical gauge like the Coulomb gauge are isolated in the cut and virtual self-energy diagrams -just the diagrams that you would draw to illustrate a parton splitting to two partons. Accordingly, the present research program was initiated in Ref. [2], which shows how to do NLO calculations in the Coulomb gauge. I should point out that, although there is a certain conceptual simplicity that emerges through the use of a physical gauge, one can use any gauge, as the recent work of Frixione, Nason, and Webber on the same subject demonstrates [3,4].
I begin with a brief restatement of the problem. Consider the calculation of an experimental observable that describes a hard scattering process involving strongly interacting particles and suppose that the chosen observable is infrared safe, which makes a perturbative calculation for the observable sensible. The simplest sort of calculation for this purpose consists of calculating the hard process cross section at the lowest order in α s , call it α B s , at which it occurs. However, one often performs a next-to-leading order calculation, including terms proportional to α B s and α B+1 s . Then the estimated error from yet higher order terms (which are usually unknown) is typically smaller than with a leading order calculation.
The NLO calculations are typically in the form of computer programs that act as Monte Carlo generators of partonic events. They produce simulated events with final states f consisting of a few partons with specified momenta. For example, for three-jet production in electron-positron annihilation, there are three or four partons in the final state. Each event comes with a weight w n . Then the predicted value for an observable described by a measurement function S(f ) is w n S(f n ). (1) A serious weakness of most current NLO Monte Carlo generators is that they produce simulated final states that are not close to physical final states. This contrasts sharply with the highly successful leading order Monte Carlo event generators such as Pythia, Herwig, and Ariadne [5][6][7], which do generate realistic final states. What could be done to generate realistic final states in an NLO calculation? One might think of using the three and four parton final states in the NLO partonic generators to generate a parton shower from each outgoing parton. This gives a final state with lots of partons. Then one could use one of the leading order Monte Carlo event generators to turn the many partons into hadrons. There should not be a serious problem with the hadronization stage, since this is modelled as a long distance process that leaves infrared safe observables largely unchanged. The problem lies with the parton showers. Here a high energy parton splits into two daughter partons, which each split into two more partons. As this process continues, the virtualities of successive pairs of daughter partons gets smaller and smaller, representing splittings that happen at larger and larger distance scales. The late splittings leave infrared safe observables largely unchanged. However, the first splittings in a parton shower can involve large virtualities. They represent a mixture of long distance and short distance physics. Thus one must be careful that the showering does not reproduce some piece of short distance physics that was already included in the NLO calculation. To be a little more precise, if one expands the prediction of the calculation including showers in a power series in α s (Q), where Q is the hard process scale, in the form then our goal is to insure that the coefficients C 0 and C 1 are exactly the same in the calculation with showers as in the purely perturbative NLO calculation. The algorithms that we propose apply to showers from final state massless partons. They are illustrated using the process e + + e − → 3 jets in quantum chromodynamics. 1 While our research program was underway, Frixione and Webber succeeded in matching showers to a NLO calculation for a problem with massless partons in the initial state (but no observed colored partons in the final state) [3]. Subsequently, Frixione, Webber, and Nason have extended this method to a problem with massive colored partons in the final state [4]. Their method is similar to the method presented in Ref. [1], or more precisely to a variant of that method in which the first level of parton splitting, based on Coulomb gauge splitting functions in Ref. [1], is instead based on the splitting functions of Herwig. The soft gluon effects treated in this paper are not part of the algorithms of [3,4]. Instead, these references employ what amounts to a soft gluon cutoff. 2 There has been other work on this problem [8]. That of Collins [9] is particularly instructive.

II. THE PROBLEM OF THE SOFT SINGULARITIES
In this section, I describe the particular problem left over from Ref. [1] that remains to be addressed in this paper, recalling just enough of the argument of Ref. [1] to set up the issue.
Consider one of the Born graphs for e + + e − → 3 jets. Three partons emerge into the final state. The idea presented in Ref. [1] is to let each of these partons split into a pair of partons, as illustrated in Fig. 1. Two features of the splitting are important. First, when the virtualityq 2 of the pair is small, the splitting probability is proportional to (P (x)/q 2 )dq 2 dx where x is the share of the momentum carried by one of the daughter partons and P (x) is the appropriate Altarelli-Parisi splitting function. Second, the collinear singularity atq 2 → 0 is damped by a Sudakov factor with the behavior exp(−α s c log 2 (q 2 )) forq 2 → 0. Thus each of the three partons makes a jet, but in the limit of small α s , each jet is usually very narrow and appears in an infrared safe measurement like a single massless parton. There is a qualifying adverb "usually" here. A fraction α s of the time, one of the splittings has a 1 I discuss e + + e − → 3 jets but not e + + e − → 2 jets. Thus the present calculation is to be applied with a measurement function that give zero for events that are close to the two-jet configuration. 2 The parameter C 1 above is then modified by an amount proportional to a positive power of the cutoff parameter β, so as long as β is chosen small enough, the cutoff effects in [3,4] are negligible.
FIG. 1: Primary parton splitting from Ref. [1]. The filled circles represent graphs for the Born amplitude and complex conjugate amplitude. Each of the partons emerging from the Born amplitude splits into two partons with a vertex, represented by the squares, that includes a Sudakov suppression factor. Each of the six daughter partons undergoes further, secondary, splittings and enters the final state as a complete shower. The secondary splittings are represented by the diamonds. substantial virtuality and we get a four-jet final state. Thus there is an order α B+1 s effect in which some probability is removed from the three-jet final state and given to a four-jet final state. In a purely NLO calculation, this effect is included in the α B+1 s graphs. Thus, to keep the calculation correct to NLO, we subtract these probabilities from the α B+1 s perturbative graphs.
The subtractions are required for NLO consistency of the calculation, but they have another important effect. Since the splitting functions have just the structure found in the order α B+1 s graphs in the collinear limit, the subtractions cancel the collinear singularities (and the collinear × soft singularities) of the order α B+1 s graphs. The real and virtual singularities are separately cancelled, whereas without the subtractions the real and virtual singularities are left to cancel against each other (which, of course, they do 3 ). If there were no other singularities, we could let the three and four parton final states of the order α B+1 s graphs develop into parton showers, which would affect the result of the calculation only at order α B+2 s . The calculation would be numerically stable because the integrals for the three and four parton states would be separately finite.
But, of course, there are other singularities, namely the singularities from emission of soft gluons at wide angles from the three hard partons, as well as the corresponding singularities in the α B+1 s virtual graphs. These singularities cancel between the order α B+1 s graphs with 3 We employ a calculation [10][11][12] in which the integrals over the three-momenta in virtual loops are performed numerically. Then this cancellation happens automatically point by point in the integration. In other methods, the virtual integrals are performed entirely analytically. In this case, the cancellation mechanism is less simple. three and four final state partons, but we cannot let these two kinds of final states develop separately into showers without losing this cancellation and thus encountering numerical instability. This is not just a technical problem to be somehow sidestepped. The emission of a soft gluon at order α B+1 s is a real physical effect. It is an effect that does not fit a picture of separate evolution of the three jets because the soft gluon is sensitive to the whole final state. That is, the probability to emit a soft gluon does not break up into the sum of three independent probabilities, one for emission from each jet.
There is a natural way to deal with this. We can consider the three jets to form an antenna that radiates a soft gluon. The radiation probability is proportional to (f (θ, φ)/E)dEdΩ, where E is the energy of the gluon, θ, φ are its angles, and f (θ, φ) is a function that is determined by the directions of the three jets. We can account for the normalization of the radiation probability provided by the virtual graphs by providing a suppression factor with the behavior exp(−α s c | log(E)|) for E → 0. This is depicted in Fig. 2. Thus a soft gluon is always emitted, but when α s is small the gluon is usually too soft to matter to an infrared safe measurement. However, a fraction α s of the time the gluon can have a substantial energy and we get a four-jet final state. Thus there is an order α B+1 s effect in which some probability is removed from the three-jet final state and given to a four-jet final state. In the purely NLO calculation, this effect is included in the α B+1 s graphs. Thus, to keep the calculation correct to NLO, we subtract these probabilities from the α B+1 s perturbative graphs. 4 Of course, the required subtraction terms will be just the soft singularity subtractions that we needed to stabilize the α B+1 s contributions if we add showers separately to α B+1 s graphs with three and four-jet final states.
Does the necessity to add soft radiation from a three-jet antenna mean that the standard parton shower Monte Carlo programs such as Herwig, Pythia, and Ariadne need to be rewritten in order to be compatible with an NLO calculation? Not at all. We do not need soft radiation with its correct angular dependence at every stage of a parton shower. All that we need is to radiate one soft gluon from the hard event with its correct angular dependence. This radiation, along with the initial splittings of the three hard partons, can be considered to be part of the NLO calculation, so that the Monte Carlo programs do not need any modification at all. I will add more comments on this subject in the conclusions section of this paper.

III. SOFT GLUON EMISSION
In Ref. [1], we saw that one can turn "Born × [1 + real − |virtual |]" into "Born × real exp(−|virtual |)" with respect to self-energy graphs, which contain the collinear and collinear × soft singularities in the Coulomb gauge. In this section, we begin a detailed description of the singularities generated by the emission of soft gluons at wide angles to the outgoing partons. Up to now in our description, these singularities cancel between real and virtual emissions in the form "Born ×[1+real −|virtual |]." In the following sections, we turn this into the form "Born × real exp(−|virtual |)" together with some left over nonsingular terms.  Fig. 1 are also used here. The extra gluon coming from the circles on the right and left represents the soft gluon radiated from the three jets. As with the other partons, the soft gluon undergoes secondary splittings and generates a shower, represented by the diamonds.
Consider a cut Born graph in which the final state has a quark (momentum q 1 ), a gluon (momentum q 2 ) and an antiquark (momentum q 3 ). These partons will make collinear showers, but we ignore that for now. What concerns us here is that these outgoing partons act as an antenna that radiates soft gluons. Let a soft gluon with momentum l be emitted by parton i and, in the complex conjugate amplitude, absorbed by parton j.
We need to specify the kinematics in a definite way. We choose to route the momenta for the {i, j} term so that the momenta { p 3 } carried by, respectively, the quark, the gluon, and the antiquark in the final state are given by Thus for a graph with i = j, the momentum k before the soft gluon emission is k = q and afterwards it is p = q − l. For a graph with i = j, the momentum k before the soft gluon emission is k = q + 1 2 l and afterwards it is p = q − 1 2 l.

A. The soft gluon approximation
We now specify the approximation to be applied for the soft gluon emission [13]. Let a parton with momentum k emit a soft gluon with momentum l and enter the final state with momentum p. Then l µ = k µ − p µ . The momenta p and k are approximately in the direction of a light-like vector u. We take u = (1,q) whereq = q/| q |, and q is 1 2 ( k + p) in the case of an i = j graph or k in the case of an i = j graph, as specified in Eq. (3).
The approximation is based on the fact that the soft gluon sees a current J µ that is approximately proportional to u µ . Thus if we denote the polarization vector of the gluon by ǫ µ , we can replace This replacement becomes exact in the limit in which J µ approaches λu µ for some constant λ. It is important the the momentum l not point along u. That is, the gluon momentum should be small so that both k and p and thus J are approximately proportional to u, but the gluon momentum should not be collinear with the momentum of the emitting parton. For the graph in which the soft gluon goes from parton i to parton j, the factors u · ǫ/u · l, together with appropriate color matrices, give a factor (in the Coulomb gauge) Here D µν (l) is the Coulomb gauge numerator factor wherel is (0, l 1 , l 2 , l 3 ). The constants C ij are the color factors for this process, What about the factor J · l? In the case of emission from a quark, we have (since so we simply include a factor / p for the final state quark. The final state phase space associated with the original graph has a factor 1/(2| p |). We retain this factor without change.
For emission from an antiquark we have [noting that the momenta p and q are directed against the fermion number flow and that we should extract a factor (−1) to associate with the color factor for the graph] so we simply include a factor −/ p for the final state antiquark. The final state phase space associated with the original graph has a factor 1/(2| p |). We retain this factor without change.
The case of emission from a gluon is a bit more complicated. We denote the tensor associated with the vertex function by Then the three gluon vertex function with momenta The vertex dotted into (k−p) times the adjoining propagator and the propagator numerator for the final state gluon, with the −ig and the color matrix −if abc removed, is Here in the second to last step we have noted that D(p) µ αp α = 0. Then in the last step we drop the term proportional to l so as to simplify our soft gluon approximation. Thus the net result of the approximation is to include, in addition to the factor u · ǫ/u · l, a factor D(p) µν for the final state gluon. The final state phase space associated with the original graph has a factor 1/(2| p |). We retain this factor without change.
It is important that the factors in Eqs. (8), (9) and (11) revert to the factors in the Born diagram when l → 0.

B. The added terms
We can now apply these ideas. We write the contribution to the cross section from a given cut Born level graph as Here S 3 is the final state measurement function [see Eq. (4) of Ref. [1]] and f 3 contains the factors from the Feynman rules corresponding to the three final state particles, which, suppressing the Dirac and Lorentz indices, are The function G 3 contains everything else, including an amputated Green function times a hermitian conjugate amputated Green function. With this notation, we write the terms to be added for real soft gluon emission as where p (ij) k are defined in Eq. (3) and S 4 is the measurement function for four final state particles.
The factor θ( l 2 < M 2 soft ) is inserted to suppress contributions from momenta l that are outside the range of validity of the soft gluon approximation. We take where √ s 0 and t 0 are, respectively, the c.m. energy and the thrust of the final state with parton momenta { q 1 , q 2 , q 3 } and where we take the parameter λ soft to be 1/3. The functions F ij are, for i = j, as specified in Eq. (5). For i = j we begin with a factor We use Eq. (7) to break C ii up into Then in the term proportional to C ia , we impose the condition that the angle between l and q i should be greater than the angle between q a and q i . That is, our soft gluon for i = j should be soft but not collinear, and we impose a cut to enforce that condition. Thus we take It will be useful to have at hand an abbreviated notation based on the definitions and With this notation, the Born and soft-real contributions are and The indicated integration overl is an integration over the unit sphere.
We define a corresponding virtual soft gluon contribution that will cancel the real soft gluon contribution for small l as long as the measurement function is infrared safe, so that R ij → R 0 for l → 0: Conveniently, the integral over the angles of l in I[soft, virtual] is simple: The integral is largest whenq 1 ·q 3 = 1, which impliesq 1 ·q 2 =q 3 ·q 2 = −1. Then the integral is 2N C = 2C A . The integral is smallest whenq 1 ·q 2 = 1, or elseq 3 ·q 2 = 1. Then the integral is N C − 1/N C = 2C F . We will, in due course, see how to exponentiate the added contributions I[soft, real] + I[soft, virtual]. First, however, we must subtract these terms from I so that the cross section remains unchanged.

IV. ADJUSTING THE COLLINEAR TERMS
The terms in Eqs. (23) with i = j represent soft gluon approximations to the cut selfenergy diagrams for outgoing partons i. Thus we should subtract them from the exact cut self-energy diagrams. Consider, for example, i = j = 1, corresponding to a self-energy insertion on the outgoing quark propagator. This contribution has the form given in Eq. (7) of Ref. [1], where q = q 1 and the notation, including the variables (q 2 , x, φ), is explained in Sec. III of Ref. [1]. We can translate the {1, 1} term in Eq. (23) into the notation used in this equation by using the change of variable Then Now we can use to rewrite this as Thus we can subtract I 11 [soft, real] by replacing M g/q (q 2 , x, φ) by where The effect of this is to remove the part of M that is singular when | l | → 0 at a fixed angle greater than the angle between the quark and the other two outgoing partons. Note that | l | → 0 at a fixed angle corresponds to ∆ → 0 and x → 0 with ∆/x fixed.
The term in Eq. (23) for emission from the antiquark line, I 33 [soft, real] is given by Eq. (31) with (/ q − / l) replaced by (−/ q + / l) and the index 1 replaced by 3. Thus we can subtract I 33 [soft, real] by replacing M g/q (q 2 , x, φ) in the analogue of Eq. (26) by with the function P soft 3 (q 2 , x) given in Eq. (33) with i = 3. Finally, the term in Eq. (23) for emission from the gluon line, I 22 [soft, real] is given by with the function P soft 2 (q 2 , x) given in Eq. (33) with i = 2. Thus we can subtract I 22 [soft, real] by replacing M µν g/g (q 2 , x, φ) in the analogue of Eq. (26) by (There is also a g → qq splitting function M µν q/g , but there is no soft singularity for g → qq so this function does not get modified.) The terms in Eq. (24) with i = j represent soft gluon approximations to the virtual self-energy diagrams for outgoing partons i. Thus we should subtract them from the exact virtual self-energy diagrams. Consider the virtual self-energy insertion on the outgoing quark propagator. This contribution has the form given in Eq. (12) of Ref. [1], The functions P a/b are given in the Appendix of Ref. [1]. When translated into the notation used in this equation, the {1, 1} term in Eq. (24) is where P soft 1 is defined in Eq. (33). Thus we can subtract I 11 [soft, virtual] by replacing P g/q (q 2 , x) by Similarly, we can subtract I 33 [soft, virtual] by replacing P g/q (q 2 , x) = P g/q (q 2 , x) in the analogue of Eq. (37) by and we can subtract I 22 [soft, virtual] by replacing P g/g (q 2 , x) in the analogue of Eq. (37) by After the collinear splitting terms have been adjusted to remove their wide angle soft parts, we turn "Born × [1 + real − |virtual |]" for the collinear splittings into "Born × real exp(−|virtual |)" as explained in Ref. [1]. The only difference is that the meaning of "real " and "virtual " here have been slightly changed.

V. ADJUSTING THE NLO GRAPHS
It remains to subtract the terms I ij [soft, real] and I ij [soft, virtual] for i = j from I. Recall that we include in the calculation the cut order α B+1 s graphs that do not have a cut self-energy subdiagram or a virtual self-energy subdiagram with the immediately adjacent propagator cut. (We also include diagrams with cut self-energy subdiagrams, but with a theta function that requires the virtuality in the self-energy subdiagram to be large.) The included graphs have soft gluon divergences, which cancel between real and virtual soft gluon emissions. We will subtract the approximate soft gluon contributions from these exact contributions, so that the soft gluon divergences are cancelled separately for the real and virtual contributions.

A. Real soft gluons
This is simple for the real emission terms, I[soft, real]. Consider a cut order α B+1 s graph with four parton lines cut. The cut lines are a quark and an antiquark plus either two gluons or another quark and antiquark. In the case of a quark and an antiquark plus another quark and antiquark, we do nothing. In the case of a quark and an antiquark plus two gluons, we designate one of the gluons as potentially soft. Then this cut graph corresponds to the contribution from one of the possible Born graphs to one of the terms in I[soft, real], which has two gluons in the final state, one of them "hard" and one soft. We simply subtract this contribution to I[soft, real] from the cut graph in question. Now we designate the other final state gluon as potentially soft, find the corresponding contribution to I[soft, real] and subtract it.
This procedure has two effects. First, we use up all of the terms in I[soft, real] that we needed to subtract from I. Second, the cut next-to-leading order graph with these counterterms has only integrable soft-gluon singularities since the singularities of the counterterms match the singularities of the cut order α B+1 s graph when either of the gluons become soft. Furthermore, the cancelling terms have the same final states: we are cancelling graphs with four final state partons against counterterms with the same four final state partons.

B. Virtual soft gluons
Subtracting the virtual terms, I[soft, virtual] requires a bit more analysis. We want to subtract these terms from the cut next-to-leading order graphs with three partons in the final state -a quark, a gluon, and an antiquark. Such a cut graph has a virtual loop either to the right of the final state cut or to the left. We want the contributions from −I[soft, virtual] to cancel the soft gluon singularities from these virtual loops. However, the terms in I[soft, virtual], as written, do not have the right structure to effect the cancellation. Thus we need to rearrange the terms.
In the {i, j} term in I[soft, virtual] we have the function In the {j, i} term we choose to redefine the integration variable so thatl → −l. Then we have a contribution with We use the identity where and The function F L ij gives the singular factor that arises from the exchange of a soft virtual gluon in a loop graph to the left of the final state cut. Similarly, F R ij gives the soft gluon singularity for a loop graph to the right of the final state cut. The essential part of a demonstration of this results from considering the integral Here we have used the notation of Eq. (5) for the emission and absorption of a real gluon and have applied the same approximations. Upon performing the integration, we find In F L ij (l), the two parts of the first term represent respectively the process in which parton j emits a (transversely polarized) soft gluon that is absorbed by parton i and the process in which parton i emits a soft gluon that is absorbed by parton j. The last term is the phase produced by the Coulomb force between the outgoing partons.
With these functions at hand, we consider a cut next-to-leading order graph with three parton lines cut. Then there is a virtual loop. In the virtual loop there may be one or two gluons that connect two lines that go to the final state. Designate one of these gluons as potentially soft. This graph then corresponds to one of the contributions to I[soft, virtual]. We subtract a term of the form with J = L or J = R depending on whether the virtual loop was to the left of the cut or to the right. We do this for each potentially soft gluon in the virtual loop. This procedure takes care of the terms in I[soft, virtual] that we needed to subtract from I. Additionally, the singularity of the counterterms matches the singularity of the graph when either of the gluons become soft (see [11]) so that the graph with these counterterms has only integrable soft-gluon singularities.
There are some technical issues to attend to. First, one should match the integration variables in the loop graph to the integration variables in the counterterms, so that { q 1 , q 2 , q 3 } are the momenta of the final state partons and l is the momentum of the potentially soft gluon in each case. Second, one needs to deform the integration contour for the counterterms according to the iǫ prescriptions indicated. This deformation should closely match that of the next-to-leading order graph in the limit of small l, as discussed in the Appendix.

VI. SECONDARY SHOWERING FROM NLO GRAPHS
We have started with the order α B+1 s graphs simply as specified by the Feynman rules and then, because there are order α B+1 s effects incorporated in the splittings of the partons emerging from the α B s graphs and in the soft gluon radiation from the α B s graphs, we have subtracted certain quantities from the α B+1 s graphs. The result is that the subtracted α B+1 s graphs with and without virtual loops are separately free of infrared divergences.
One can attach a shower to each final state parton emerging from a subtracted order α B+1 s graph. Each shower affects the expectation value of the observable being calculated. However, according to the argument in Ref. [1], as long as the observable is infrared safe, the effect is to add terms proportional to α B+2 s and higher powers of α s . This holds independently of the exact functions used to generate the showers, as long as these functions have the required generic properties. For our present purposes, we can use the simple algorithm described in Sec. VIII of Ref. [1].
The only thing that we need to specify is the choice for width parameter κ for the first splitting in the shower. We will take where q is the momentum of the parton that is to split and t is the thrust of the partonic final state in the order α B+1 s graph. The factor (1 − t) is important. In the case that we generate an event with thrust very near to 1, we do not want the showering to turn this two-jet event into a three-jet event, an event with thrust substantially less than 1. The factor (1 − t) prevents this. [For thrust greater than a value t max close to 1, the code [14] provides a stronger cutoff: Here the three outgoing partons radiate a soft gluon according to the radiation pattern specified by the F ij . The gluon radiation at small l is suppressed by the exponential function, which gives the probability that the gluon was not already radiated with a higher energy than | l |. As we found in the companion paper [1] on collinear emissions, there is a subtraction to multiplication theorem that says that The proof of this is simple. We have where and Now I[radiate, 1] can be expanded in powers of α s to give where Here C is the angular integral in the exponent, given in Eq.
This completes the proof of the theorem.

VIII. THE RESULT
We are now able to generate a useful expression for the observable I. We begin with a cut graph that contributes to I[Born], Eq. (12). For each of the final state partons, we introduce an additional integration, representing the splitting of that parton, as explained in Ref. [1]. We also introduce an integration, representing the radiation of a soft gluon from the final state partons. The function representing the Born matrix element then depends on all of the integration variables. This gives a probability to create a total of seven partons. Each of these seven partons then serves as the source of a secondary shower. 5 This is illustrated in Fig. 2.
In addition, there are contributions from order α B+1 s cut Feynman graphs plus their soft gluon counterterms. Included here are cut self-energy graphs, but only for virtuality above the cutoff λ V q 2 i , so that there are no collinear divergences. Each of the three of four partons emerging from such a graph then serves as the source of a secondary shower. These order α B+1 s remnant terms are complicated, but they have only integrable singularities and they are just what is needed to leave the calculation correct to next-to-leading order.

IX. A NUMERICAL TEST
In this section I present some numerical results based on the algorithm described above. I have computed the thrust distribution dσ/dt for thrust t = 0.86, in the middle of the three-jet region. We will compare this distribution calculated with parton showers to the same distribution calculated with a straightforward NLO computation with no showers. The difference should be of order α B+2 s . We divide the difference by the NLO result, forming The ratio R should have a perturbative expansion that begins at order α 2 s . We can test this by plotting the ratio against α 2 s . We expect to see a curve that approximates a straight line through zero for small α 2 s . For comparison, I exhibit also the ratio with the order α B s graphs with showers but with the α B+1 s corrections omitted, The ratio R LO should have a perturbative expansion that begins at order α 1 s . Thus we expect to see a curve proportional to the square root function [α 2 s ] 1/2 for small α 2 s . The size of the coefficient of [α 2 s ] 1/2 in R LO does not have any great significance, since it is quite sensitive to the choice of renormalization scale.
The results of this test are shown in Fig. 3 for √ S = M Z with a renormalization scale µ = √ S/6. The values of α s (M Z ) range from (1/2) × 0.118 to 2 × 0.118, where 0.118 represents something close to the physical value for α s (M Z ). One should note that 2 × 0.118 amounts to quite a large α s in this calculation since α s (M Z ) = 0.236 gives α s (µ) = 0.586.
We see the expected shape of the R LO curve. We then examine whether the R curve approaches a straight line through the origin as α 2 s → 0. Within the errors, it does. However, the slope of the straight line is quite small. Presumably with other choices of parameters in the program, the absolute value of the slope would be bigger.
In Fig. 4, I show the same comparison, but this time for t = 0.71. This is near the value t = 2/3 that marks the far end of the three-jet region, with the three partons in what is sometimes called the Mercedes configuration. The results are similar to the results for t = 0.86.
In Fig. 5, I show the same comparison, but this time for t = 0.95. This is near the two-jet limit at t = 1. For α s = 0.118, one would normally not use a calculation that did not include a summation of logs of 1 − t for t this close to 1, since log(0.05) 2 ≈ 9. Thus I would not recommend using the code discussed in this paper for a comparison to data this near to the two-jet limit. Nevertheless, we can still test for the absence of an α 1 s term in R. Looking at the graph, we see that, within the errors, there is no evidence for a nonzero α 1 s term in R. We expect an α 2 s term, but it appears that the coefficient of α 2 s is quite small. On the other hand, it appears that some of the yet higher order terms are quite substantial.

X. AREAS WHERE MORE WORK IS NEEDED
The present work suffers from a number of deficiencies. First, the code is designed to calculate three-jet cross sections in electron-positron annihilation correctly to NLO, but it ignores the two-jet cross section. It should be rather simple to calculate the two-jet cross section at next-to-leading order with showers. Then one should merge the two calculations so that three-jet observables are calculated correctly at order α 2 s while the two-jet cross section is calculated correctly at order α 1 s . Such a merged calculation would have benefits even for the calculation of three-jet observables. Specifically, there would be an advantage when calculating three-jet observables close to the two-jet region. For example, consider the thrust distribution near to t = 1. The two-jet calculation with showers, before any merging, would approximately sum terms proportional to α N s log J (1 − t), generating these terms by showering from the two-jet hard process. For dσ/dt, the merged calculation should contain the α 1 s and α 2 s contributions correctly, including of course the term proportional to α 2 s log(1 − t). It should also include the terms α n s ×logs for n > 2 generated by showering from a two-jet configuration. The three-jet algorithm presented in this paper does not include merging with a two-jet calculation and thus does not generate these logarithms. Thus the present program, when applied to the thrust distribution, should not be used for small values of 1 − t.
Second, it is desirable to add not only parton showers but hadronization to NLO perturbative calculations. For hadronization, one can use an existing model [5][6][7]. However, most hadronization models require color information on the final state partons. Since the current code does not assign a color structure to these partons, it will be necessary to add this capability.
Third, in principle, hadronization models should be infrared safe in the sense that final partonic states that differ by a parton splitting into two almost collinear partons (or one parton and another with almost zero momentum) will produce almost the same hadronic final state. Since hadronization models currently in use do not necessarily have this property, the NLO calculation should provide it by combining parton pairs with virtuality less than some cutoff parameter matched to the hadronization model to be used.
Finally, the current paper includes hardly any numerical investigations to check how well the program works for various observables and how the results depend on parameters such as the renormalization scale µ and the soft gluon limiting energy M soft . Addressing these needs remains a topic for future research.

XI. CONCLUSIONS
In this paper and [1], we have seen how to add parton showers to a next-to-leading order calculation in QCD in such a way that the result obtained for an infrared safe observable remains correct to next-to-leading order. The main feature of this procedure is to turn the collinear singularities of the order α B+1 s graphs into the primary parton splittings, that is the first splittings of the partons that emerge from an α B s graph. This was accomplished in [1] in a particularly simple fashion by virtue of the fact that Ref. [1] is based on a NLO calculation in the Coulomb gauge. In this gauge the collinear and collinear × soft singularities are entirely contained in the cut self-energy graphs. Thus we simply eliminate the small-virtuality part of cut self-energy subgraphs in α B+1 s graphs and instead incorporate these contributions into the primary parton splittings from the corresponding α B s graph. One could subtract a different quantity from the α B+1 s graphs as long as it has the same collinear and collinear × soft singularities, 6 as in the work of Refs. [3,4].
The next step, and the main subject of this paper, is to treat the soft, wide angle singularities of the theory. Again, we subtract terms with these singularities from the α B+1 s graphs and then incorporate them as radiation from the α B s graphs. Here, the hard partons emerging from an α B s scattering form an antenna that radiates the soft gluon. The gluon is radiated coherently from the antenna and is not part of any jet.
This gives an algorithm that can be viewed as a simulation of the production of three hard partons (for the three-jet process considered in this paper) and one soft gluon, with each of the three hard partons splitting into two daughter partons. The six daughter partons and one soft gluon then serve as the progenitors of parton showers according to the simple algorithm used in Ref. [1] or any standard algorithm. There remains the remnants of the order α B+1 s graphs, with three and four parton final states. Each of these final state partons serves as the progenitor of whatever sort of parton shower is desired.
The results of this paper suggest a strategy for building future NLO calculations matched to parton showers. The matching concerns the primary parton splittings, that is the first splittings of the partons emerging from a Born graph. The matching also needs the radiation of a soft gluon from the antenna formed by partons emerging from the Born graph. This suggests that the author of an NLO calculation might build in the primary parton splittings and the soft gluon radiation. There are also secondary splittings. These are the splittings from the soft gluon produced from a Born graph, the splittings of the daughters produced by the primary splittings, the splittings of the partons from an order α B+1 s graph, and, finally, the further splittings of all of these partons as a complete parton shower is formed. The secondary splittings do not need to be matched to the NLO calculation. That is, the soft gluon and the partons after the primary splittings and the partons produced by the α B+1 s graphs can be fed to any shower and hadronization algorithm as long as it is infrared safe. The result is that for a given process such as e + + e − → 3 jets we can have N NLO next-to-leading order calculations and N MC shower and hadronization calculations and we do not need N NLO × N MC matching schemes as long as each NLO calculation includes its primary parton splittings and treatment of soft gluon radiation.
This does not preclude debate over how best to do the primary parton splitting and soft gluon radiation. This paper and Ref. [1] have, of necessity, chosen one scheme, but that scheme may well be less than optimal.
The code described in this paper is available at [14].
R to remain unchanged. This is in contrast to our treatment of showering from the Born graphs and is required for the NLO graphs in order to preserve the soft gluon cancellations. Next, in the integral of the full graph f , the integration is deformed 7 into complex loopmomentum space, The imaginary part of the deformed loop momentum is a definite function κ( l) of the real part. The Jacobian J is the determinant of the matrix (δ IJ + i∂κ I /∂l J ). By deforming the contour, we avoid the singularity of the form in f , the {i, j} scattering singularity. It is, however, not allowed to deform the contour away from the soft singularity at l ij = 0. For that reason κ( l ) is proportional to l 2 ij as l 2 ij → 0. Hence we come very near to the scattering singularity when l ij is small.
What about the subtraction term? Here there is a singularity of the form We may call this the approximate scattering singularity. We note that for l ij → 0, so the singularity (A.4) is indeed an approximation to Eq. (A.3) if we identify l ′ with l ij . As in the f integral, we need to deform the contour to avoid the singularity (A.4).
As in the f integral, the deformation will have to vanish as we approach the soft singularity at l ′ = 0. In order for the cancellation between the full graph and the soft gluon subtraction to work, we need to be careful about the matching of the singularities in the two integrands.
We are now in a position to state the matching problem more precisely. We define complex vectors s ij , one for each potentially soft gluon line {i, j} in Eq. (A.1). Each s ij is a function of the corresponding l ij that is asymptotically equal to l ij when l ij is small. Then we write the integral in Eq. (A.1) in the form The Jacobian J s is the determinant of the matrix ∂s I ij /∂l J ij . Our goal will be to specify the functions s ij so that we avoid the approximate scattering singularity in the second term of Eq. (A.6) and so that the l ij → 0 singularities in the sum of the two terms cancel sufficiently well that the integral is convergent.
We will approach this problem in stages. First, we concentrate on the l ij → 0 region.
We immediately recognize a small problem. The soft singularity (with κ = 0) lies on an ellipse, while the approximate soft singularity lies on a plane. To make them match at the required accuracy, we define s ij = l ij + (1 − l 2 ij /M 2 soft ) ξ( l ij ) + i ζ( l ij ) (to be improved), (A.7) where ξ and ζ are certain functions of l ij to be defined below and M 2 soft is the maximum value of l 2 in the soft gluon subtraction, Eq. (14).
The idea is to have the soft subtraction cancel the leading singular behavior of f point by point in l. To do this, we need to match the singularity There are three parts in the relation (A.7) between s ij and l ij . There is a real displacement ξ, a contour deformation ζ in the imaginary direction, and a factor (1− l 2 ij /M 2 soft ). The factor (1 − l 2 ij /M 2 soft ) serves to set the deformation to zero at the edge of the integration region and needs no further discussion. We choose the real displacement to be This choice is motivated by carrying out the expansion (A.4) to one more order. Then one finds that ( l ij + ξ( l ij )) · (q j −q i ) ∼ | q i | + | q j | − | q i + l ij | − | q j − l ij | + O(l 3 ij ). (A.11) The imaginary displacement is ζ( l ) = C ζ l 2 min (1 +l ·q i ), (1 −l ·q j ) (q j −q i ). (A.12) Here C ζ is a parameter with dimensions of inverse momentum. This function matches the function κ( l ) in the limit of small l as long as one chooses C ζ to be the same parameter as is used in κ( l ). There are other singularities that have to match. First, both f and the soft subtraction terms have 1/| l ij | 3 soft singularities that evidently match since s ij = l ij to lowest order. Second, both terms have collinear singularities along the lines 1+l ij ·q i = 0 and 1−l ij ·q j = 0. These are pinch singularities: it is not allowed to deform the contours to avoid them. For this reason ζ( l) vanishes when l is in the same direction as − q i or q j . Furthermore, the part of ξ that is not proportional to q i vanishes when l is in the same direction as − q i and the part of ξ that is not proportional to q j vanishes when l is in the same direction as q j .
We are now ready for the second stage of our construction. In the first stage, we concentrated on the l ij → 0 region. The function s ij given in Eq. (A.7) works fine for small l ij . However, our integration region ends at | l ij | = M soft , so we supplied a factor that turns the contour deformation off at the boundary of the integration region. 9 Unfortunately, the singularity at l ij · (q j −q i ) = 0 intersects the surface | l ij | = M soft . Thus we should retain some contour deformation even at the boundary | l ij | = M soft of the integration region. This is simple to do.
Let us consider the problem anew, this time ignoring the special problems at | l | → 0 that were addressed in the analysis above. Let n = (q j −q i )/|q j −q i | and write l ij in spherical polar coordinates (l, θ, φ) based on n as the z axis. Then we have an integral where x = cos θ, that is x = l ij · n/| l ij |. (A.15) The integrand has a singularity of the form 1/(x + iǫ). All that we need to do is deform the contour for the x-integral away from the singularity. There are collinear singularities at l ij ·q j = 1 andl ij ·q i = −1, that is at φ = 0, π and x = x 0 , where The factor l 2 ij /M 2 soft is included here simply to turn this deformation off for small l 2 ij , since we already have a method for dealing with the small l 2 ij region. In a coordinate frame independent notation, this is l ij → s ij with s ij = l ij + η( l ij ) (to be improved), (A.20) Now we have two methods for moving the integration contour, each of which works in its region and turns off in the region of the other method. All that we have to do is to add the two deformations,