Dipole Splitting Algorithm -- A practical algorithm to use the dipole subtraction procedure --

The Catani--Seymour dipole subtraction is a general and powerful procedure to calculate the QCD next-to-leading order corrections for collider observables. We clearly define a practical algorithm to use the dipole subtraction. The algorithm is called the Dipole splitting algorithm (DSA). The DSA is applied to an arbitrary process by following well defined steps. The subtraction terms created by the DSA can be summarized in a compact form by tables. We present a template for the summary tables. One advantage of the DSA is to allow a straightforward algorithm to prove the consistency relation of all the subtraction terms. The proof algorithm is presented in the subsequent article. We demonstrate the DSA in two collider processes, $pp \to \mu^{-}\mu^{+}$ and $2\,jets$. Further as a confirmation of the DSA it is shown that the analytical results obtained by the DSA at the Drell-Yan process exactly agree with the well known results obtained by the traditional method.


Introduction
The operating CERN Large Hadron Collider (LHC) discovered a new boson whose mass is around 125 GeV. The new boson is identified as the higgs boson in the standard model. In order to specify the property of the field and the interaction precisely we need more results from the LHC experiments and the various examinations of the results must be carried by comparing to the theoretical predictions. In this article we treat with the subject of the theoretical prediction for arbitrary process at hadron collider like LHC. The calculation for the prediction consists of two ingredients, the parton distribution function (PDF) and the partonic cross section. The PDF is a process independent quantity and determined as a numerical function from the experimental data. The recent reviews about the PDF are seen in [1-4]. The partonic cross section is calculated by the perturbative expansion of the strong coupling constants of QCD. The prediction which includes only the leading order (LO) has a large dependence on the renormalization scale µ R in the strong coupling constants, α s (µ R ), and the factorization scale µ F in the PDF, f (x, µ F ). The large dependences on µ R and µ F lead to the large uncertainty on the prediction. The QCD Next-to-Leading Order (NLO) correction reduces the µ R and µ F dependences and makes the prediction more precise. In this reason we need the QCD NLO corrections for the precise prediction.
There are two main prediction schemes: matrix element prediction and showered prediction, which are mentioned at [5]. The matrix element prediction is the following: the partonic cross section includes the matrix element which represents the transition amplitude from free initial partons to free final partons. The jet algorithms like [6][7][8][9][10] define the jet observables. A jet algorithm is directly applied to the partons in the final state of a matrix element and the distributions about jet observables are compared to the experimental results. On the other hand the showered prediction is the following: The final state partons of a matrix element are showered by a shower algorithm. Then a jet algorithm is applied to the partons after showering and finally the distributions about jet observables are compared to the experimental results. The hadronization effect may be also included for the better simulation. The merit of the matrix element prediction is to be less involved and simpler than the showered prediction. The merit of the showered prediction is to simulate the phenomena happening after the collision better than the matrix element prediction. In this article we deal with the QCD NLO corrections only in the matrix element prediction scheme.
The tradition of the QCD NLO correction at the hadron colliders probably started from the Drell-Yan process, pp/pp → µ − µ + , at the pioneer works, [11][12][13][14][15][16][17][18][19], which may be the simplest one among the hadron collider processes. At the pioneer works the analytical expression of the NLO correction is obtained by using the dimensional regularization with d = 4 − 2ǫ, throughout. All the ultraviolet (UV) and the infrared, more precisely soft and collinear, divergences are regulated as the poles of 1/ǫ and 1/ǫ 2 . The NLO correction to cross section is generally written as (1.1) where the symbols, σ R , σ V and σ C , represent the real emission correction, the virtual 1-loop correction, and the collinear subtraction term, respectively. The 1-loop matrix element in the virtual correction includes the UV divergences, which are subtracted by the renormalization program. Each of three terms, σ R , σ V and σ C , includes the integration over the phase space (PS) in d-dimension. The complete cancellation of the soft and collinear divergences from three terms can be realized only after the analytical integration of the d-dimensional PS. Then we obtain the finite physical cross section. For more complicated multi-parton processes this method encounters mainly with two difficulties: the evaluation of the 1-loop matrix element and the analytical integration of the d-dimensional PS.
Because the matrix element for multi-parton leg processes becomes complex and long expressions. The difficulty about the evaluation of the 1-loop matrix element has being solved by the technical developments. The reviews about the recent developments are seen in [20][21][22].
The new methods are invented to overcome the difficulty about the analytical integration of the d-dimensional PS: the phase space slicing method [23][24][25][26] and the subtraction method [27][28][29]. Among them the Catani-Seymour dipole subtraction procedure [27,28] in the subtraction method is quite successful and has been widely used. At the dipole subtraction procedure the subtraction terms are introduced and the NLO correction is reconstructed as σ NLO = (σ R − σ D ) + (σ V + σ I ) + σ P + σ K , (1.2) where the symbols, σ D , σ I , σ P , and σ K , represent the the dipole-(D-), I-, P-, and K-terms, respectively. The σ D and σ I subtract all the soft/collinear singularities from σ R and σ V , respectively, at the level of the square of the matrix elements. Then the PS integration of the subtracted cross sections, (σ R − σ D ) and (σ V + σ I ), can be executed in 4-dimension to be finite. The σ P and σ K are also separately finite under the 4-dimensional integration.
Now that the dipole subtraction has been applied for many processes, we can have some drawbacks. Among them we would like to point out the following difficulties: the person who has obtained the results for the NLO corrections by the dipole subtraction, sometimes has the difficulties to confirm the results, because many subtraction terms are involved and the large amount of the calculation is executed as numerical evaluation for the Monte Carlo integration. For the other person who does not calculate by himself the confirmation of the shown results is more difficult. All the soft and collinear singularities from the D and I terms must exactly cancel the singularities from the real correction σ R and the virtual correction σ V , respectively. If we use any wrong collection or any wrong expression for the D and I terms, the cancellation is spoiled. The successful cancellation provides one check on the singular parts of the D and I terms. On the other hand the P and K terms are separately finite by themselves and the check by the cancellation is impossible. In this sense the uses of the P and K terms are the place where we can easily make mistakes. In many papers it is not clarified which and how many subtraction terms are used. Of course it is too long and unreasonable that the explicit expressions of all subtraction terms are written on paper. But it is possible that the minimal information to specify the subtraction terms is shown. At least the total number of the used subtraction terms should be clearly mentioned. At many works the dipole subtraction is applied by using the automated programs. In that case the similar kinds of the criticisms should be thrown on the implementation, especially the algorithm to create the subtraction terms. At the packages, [95] and [97], the creation algorithms are not clearly defined. We can not know which and how many D, I, and P/K terms are created under a given input process. At the execution of the packages the information of the created terms is not printed out.
At the package [93] the creation algorithm is presented on the paper. At the output of the Mathematica program the information of all the created subtraction terms and the total number are separately printed out. The output codes of the D, I, and P/K terms are separately saved in the corresponding folders and can be easily identified. At all three packages any algorithm to check the consistency of all the created subtraction terms is not provided.
In order to solve one part of the difficulties and the criticisms, we really wish the clear definition of a practical algorithm to use the dipole subtraction. Although the general algorithm is given in the original papers [27,28], we wish more practical algorithm which can be directly used step by step. Also in order to automatize the procedure into a computer program we need such a practical algorithm to be applied for arbitrary process.
In consequence we wish an algorithm which provides the clear definitions of the following items: 1. Input, output, and creation order, 2. All necessary formulae on document, 3. Necessary information to specify each subtraction term, 4. Summary table of all created subtraction terms, 5. Proof algorithm of the consistency.
The proof algorithm at the last entry in the wishlist means the followings. Eqs. (1.1) and (1.2) are combined and the consistency of the subtraction terms is represented as (1. 3) The proof algorithm of the consistency means the algorithm to prove the Eq. (1.3) for arbitrary process.
The purpose of this article is to present a practical algorithm which satisfies all the requirements in the above wishlist. We actually present such an algorithm, called as the dipole splitting algorithm (DSA). At the DSA the input is all real emission processes, which contribute to an observable process like, pp → n jets. The input of one subpartonic real process like, uū → uūg, creates the output of the D, I, P/K terms, all of which have the same initial states with the input real process. In this sense the subtraction terms are sorted by the initial state partons. At the DSA all the subtraction terms are sorted also by the kinds of the splitting which each subtraction term possesses as a part. The sort by splittings is equivalent with the sort by the reduced Born process which each subtraction term possesses as a remaining part when the splittings part is removed. In order to specify each subtraction term uniquely we introduce a bijection mapping, called field mapping.
Each field mapping is made for each subtraction term. The filed mapping exactly specify the connection between the legs of the Born process reduced from an input real process and the fixed legs of the same reduced Born process. By using the field mapping we can specify each term in a compact form without confusion. Thanks to the well-defined compact form for the specification, we can summarize the created subtraction terms in tables. We will give a template for the summary tables. We intend that the person who does not create the subtraction terms by himself can specify and write down them only by reading the tables on a document without any direct communication to the author of the tables. About the last entry in the wishlist we could find a straightforward algorithm to prove the consistency of the created subtraction terms. The cancellations in Eq. (1.3) are realized among the terms which have both of the same initial states and the same reduced Born processes. According to the feature of the DSA that the subtraction terms are sorted by both of the initial states and the reduced Born processes, the groups of the subtraction terms which are canceled each other are systematically identified. In this reason the DSA makes possible the construction of the straightforward proof algorithm.
We here mention about the relation between the DSA and the algorithm implemented in the AutoDipole package [93]. The creation algorithm of the D and I terms in the DSA is essentially same with the algorithm in the AutoDipole. In this sense the DSA originates in the algorithm of the AutoDipole. In the DSA the concrete expressions are more clearly defined by using the above mentioned field mapping and the necessary information to specify each subtraction is given in a compact form. The creation algorithm for the P and K terms in the DSA is different from the algorithm in the AutoDipole. As explained in the DSA the initial states of all created subtraction terms under one input real process are same with the initial states of the real process. On the other hand in the algorithm of the AutoDipole the initial states of the P and K terms with the non-diagonal splittings are different from the initial states of the input real process.
This article is organized as : in Sec. 2 the DSA is defined. In the subsections the master formulae and the algorithm steps, the creation of the D, I, and P/K terms, and the advantages of the DSA, are explained. In the sections, 3 and 4, the DSA is demonstrated at the example processes, pp → µ − µ + and pp → 2 jets, respectively. In Sec. 5 we have a confirmation that the analytical results obtained by the DSA at the Drell-Yan process coincide with the results by the traditional method. Sec. 6 is devoted to a summary. In Appendix. A the formulae for the subtraction terms are collected. In Appendix. B the summary tables for the dijet process are shown.

Dipole splitting algorithm
In this section the DSA is defined. At Sec. 2.1 all the steps of the DSA and the master formulae are presented. At the sections 2.2, 2.3, and 2.4, the creation algorithms of the D, I, and P/K terms are explained, respectively. The advantaged of the DSA is claridied in Sec. 2.5. The formulae for all the subtraction terms are collected in Appendix A.

Definition
The DSA consists of the following steps, Step 1. List of real correction processes: (2.1) We assume that we wish to make a prediction for observables on a collider process like pp → 2 jets. Once a collider process is selected, the contributing partonic real emission processes are specified. The Step 1 is to specify all the real processes denoted as R i and to write down the list. The number of the real processes is written as n real and the number of the fields in the final states is denoted as (n + 1). For example the collider process, pp → 2 jets, has the real correction processes as The Steps 2, 3, and 4, are explained in the following sections, 2.2, 2.3, and 2.4, respectively. The last step of the DSA, Step 5, is to obtain the NLO correction as where each σ(R i ) is defined as where the f F(x a/b ) (x 1/2 ) is the PDF and the subscript, F(x a/b ), represents the field species of the initial state parton at the leg a/b, which is again defined at the next section. The symbols,σ(R i ), with the subscripts, R, D, V, I, P, and K, represent the partonic cross sections of the real correction, the dipole (D) term, the virtual correction, the I term, the P term, and the K term, respectively. Each partonic cross section is defined aŝ (2.10) In Eqs. (2.5) and (2.6) the symbol, S R i , is the symmetric factor of the real process R i and The energy is denoted as E I = p µ=0 I for I = a, b, 1, · · · , n+1. The flux factor is denoted as is the square of the matrix element of the real process R i after the average over spin and color in 4-dimension. The n s (a/b) represents the spin degree of freedom of the leg a/b at R i in 4-dimension. In Eqs. (2.7) and (2.8), the process, B1, is the abbreviation of the B1(R i ) which is made by removing one gluon from the final states of the R i as B1(R i ) = R i − g. The S B1 is the symmetric factor of the process B1, and the symbol, Φ(B1) d , represents the d-dimensional n-body PS with the flux factor as The |M virt (B1)| 2 d is the abbreviation of the quantity, (M LO (B1) M 1-loop (B1) * + M LO (B1) * M 1-loop (B1)), after the average over spin and color in d-dimension, where the M LO (B1) and the M 1-loop (B1) are the matrix elements of the LO and the 1-loop correction of the process B1, respectively. In Eqs.(2.9) and (2.10) the B j is the abbreviation of B j (R i ) and it is a reduced Born process belonging to the R i , which is precisely defined at the next section 2.2. The S B j is the symmetric factor of the B j and the Φ a (R i : B j , x) 4 is the 4-dimensional n-body PS with the scaled momentum, xp a , and with the flux factor F (xp a , In Eqs. (2.6), (2.8), (2.9), and (2.10) the concrete expressions of the subtraction terms, D(R i ), I(R i ), and P/K(R i : B j , xp a ), are presented in the following sections, 2.2, 2.3, and 2.4, respectively.
The PS integration in Eq.(2.4) is finite in 4-dimension and we here see the finite parts separately. The real correction,σ R (R i ), has the soft and collinear singularities, which are subtracted by the dipole terms,σ D (R i ). The subtracted cross section ( (2.14) The virtual correction,σ V (B1), includes the poles of the soft and collinear divergences, 1/ǫ and 1/ǫ 2 , after the subtraction of the UV divergences, 1/ǫ UV , by the renormalization program. The I-term,σ I (R i ), cancels all the soft and collinear poles in d-dimension aŝ After the cancellation of the poles, we can integrate the PS in 4-dimension aŝ (2. 16) The cross section of the P and K terms,σ P (R i ) andσ K (R i ), are finite themselves separately in 4-dimension. In this way the PS integrations in Eq.(2.4) are executed in 4-dimension.
If the selected collider process includes any jet in final state, the jet functions which are supplied by an infrared safe jet algorithm must be multiplied.
To complete the master formulae we add the LO contribution as where the sub-partonic LO processes which contribute to the selected collider process are denoted as L i and the number of the LO processes is denoted as n LO . The Φ(L i ) 4 is the 4-dimensional n-body PS including the flux factor. Then the prediction at the NLO accuracy is written as

D term creation
In this section the 'Step 2. D(R i )' is explained. The input and output of this step are written as The input process is each real process R i among all real processes, R 1 , · · · , R n real , which are specified at the 'Step 1'. The input trivially defines the set {x} with the field species F(x I ) and the momenta p I for the values, I = a, b, 1, · · · , n + 1, as The output, D(R i ), is the summation of all the created dipole terms. The creation is repeated over all the input processes, R i with i = 1, · · · , n real , and the outputs are the corresponding dipole terms, D(R i ) with i = 1, · · · , n real .
In the original paper by the Catani and Seymour [27] each dipole term is specified by three legs (I, J, K), where the pair of legs (I, J), more precisely the combied legs, IJ, is called the emitter (pair) and the leg K is the spectator. All possible combinations of (I, J, K) are chosen without duplicate from all (n+3)-legs of R i . In the DSA the creation algorithm of the dipole terms is given as 1. Choose all the possible emitter pairs (x I , x J ) from the set {x} in the order of the splittings from (1) to (7)  where the Dipole2-u, means the species of the fields, x i and x j , are the up and anti-up quarks as (F(x i ), F(x j )) = (u,ū), and the Dipole3-ū, means (F(x a ), F(x i )) = (ū,ū), for instance. The spectator x K can be a quark or gluon in either of the final or initial state.
If the spectator is in the final/initial state, we denote the case as the subcategory -1/2. Then the category of the dipole terms are further divided into the subcategories as The summation of the dipole terms which belong to the same category is written respectively as where the subcategories about the quark flavors and about the spectator in final or initial state are all summed. The dipole terms belong to Dipole 1 are summed as 25) and the summation of all the dipole terms as output is written as The concrete expression of each dipole is given at the original paper [27] in the form, where the s IJ is defined as s IJ = 2p I · p J , and the x IJK is specified in the part 'Concrete expressions'. The Bj is a Born process which is reduced from the input process R i by removing the splitting part. Then the dipole terms, D(R i : dipj), in the category Dipole j have the the reduced Born process Bj, which is made from the R i on the rules as where the symbols, g f /i and u f /i , represent a gluon and a quark in the final/initial state.
The operation, ±g f , means to add/remove a gluon to/from the final state. The other operations like, ±u f , are similarly defined. About the B2u, B3u, and B4u, the other subcategories with the other quark flavors also exist. The symbol, (T IJ · T K )/T 2 IJ , represents the operators of the color factor insertions and the V IJ,K is the dipole splitting function with the helicity correlation. The actions of the color-and the helicity-operators onto the amplitude of the reduced Born process are clearly defined by using the field mapping explained at next part.

Field mapping
Each dipole term includes the square of a reduced Born amplitude shown in Eq. (2.27).
The original (n+3)-legs of the input process R i are connected to the (n+2)-legs of the reduced Born process. In order to specify the way of the connection clearly in a compact form we introduce a bijection mapping for each dipole, called the field mapping. Per one choice of a combination, (x I , x J , x K ), we can make a new set {x} by the unification of the elements, (x I , x J ) → x IJ , and the replacement, x K → x K . To explain the definition of the set {x} precisely, we separate the dipole terms into four categories as (IJ, K) = (ij, k), (ij, a), (ai, k), and (ai, b), where the indices, i, j, k, represent the legs in the final state and the indices, a and b, represent the legs in the initial state. The four categories are sometimes called as the final-final, final-initial, initial-final, and initial-initial, dipoles, respectively The relation between the four categories and the categories defined at the previous part 'Creation order' is shown as (2.29) • Final-final : (ij, k) In the category the set {x} is defined as , P(x b ); P(x 1 ), · · · , P(x n−2 ), P(x n−1 ), P(x n )} = {P(x a ), P(x b ); P(x 1 ), · · · , P(x n+1 ), P(x ij ), P(x k )}, (2. 32) where the symbols, F(x α ) and P(x α ), represent the field species and the momenta of the elementsx α with the indices, α = a, b, 1, · · · , n. The field species are determined as F(x n−1 ) = F(x ij ), (2.33) The following other choice of the set {x} is also possible: (2.47) The field species and the momenta in this case are similarly determined. This note about the freedom of the order of the elements in the {x} stands for the following other cases, (ij, a), (ai, k), and (ai, b) as well.
• Final-initial : (ij, a) For the category of the final-initial dipoles (ij, a), the set {x} is made as The field specifies are defined in the similar way with the previous case in Eqs. F(x α ) = F(x L ) for α = b, 1, · · · , (n − 1) , (2. 53) where the F(x ij ) is same in Eq. (2.33). About Eq.(2.53) the elements,x α , with α = b, 1, · · · , n − 1, are identified with the elements, x L , asx α = x L in Eq.(2.48) for L = b, 1, · · · , n + 1, skipping the indices i and j. The momenta are also similarly defined as where the reduced momenta,p ij andp a , are defined in Eqs. (A.36) and (A.37).
where the reduced momenta,p ai andp k , are defined in Eqs. (A.41) and (A.42).
• Initial-initial : (ai, b) For the initial-initial dipoles (ai, b), the set {x} is made as The definition of the field specifies At Eq.(2.69) the definition of F(x ai ) is same with the previous case of (ai, k). Thẽ x a represent the other elements with α = 1, · · · , n. The momenta are defined as It is noted that in this case the momentum of the spectator is not changed as shown in Eq.(2.73) and the momenta of all the other elementsx α with α = 1, · · · , n are changed intok L with L = 1, 2, · · · , (n + 1), skipping the index i as shown in Eq.(2.74).
Next we treat with the definition of the set {y} which is made from a reduced Born process. Under the input R i the reduced Born process Bj which belongs to Dipole j is made by the rules shown in Eq. (2.28). One Born process, Bj, determines the set {y} with the field species and the momenta as {y} = {y a , y b ; y 1 , · · · , y n }, (2.75) , P(y b ); P(y 1 ), · · · , P(y n )}. (2.77) The number of the elements of the set {y} is (n+2) which is same with the set {x}. We can always find a bijection (one-to-one onto) mapping from the set {x} to the set {y}, , with the indices α, β = a, b, 1, · · · , n, which satisfies with two conditions, -The argumentx α and the image y β are both in final, or both in initial state.
The two conditions mean that the mapping, f , connects the elements whose species are identical, and it does not mix the elements in the initial and final states. The inverse mapping is denoted asx α = f −1 (y β ). After the construction of a mapping, the element, x α , is identified with the image, y β . Using the inverse mapping the identification of the elements are generally written as For convenience we introduce the notations, y emi and y spe , defined as y emi = f (x IJ ), and y spe = f (x K ). One example of the field mapping is shown in Fig.2, where the number of the final states of an input process R i is assumed to be, (n + 1) = 4, and the dipole term with the combination, (IJ, K) = (x 1 x 3 , x 4 ), is picked up. Using the field mapping the concrete form of each dipole term can be clearly expressed as (2.80) Using the notation, y emi and y spe , the concrete form of each dipole terms is slightly simplified as The subscript indices I, J, and K, at the quantities, s IJ , x IJ,K , and V y emi IJ,K , refer to the elements of the set {x} and the momenta, {p a , p b , p 1 , · · · , p n+1 }, in Eq. (2.23). On the other hand the indices, y emi and y spe , at the operators, T y emi , T yspe , and V y emi IJ,K , refer to the elements of the set {y}. The Casimir operator T 2 F(y emi ) is defined as the constant, C F = 4/3 in the case F(y emi ) =quark, and C A = 3 in the case F(y emi ) =gluon. The momenta inputted into the reduced Born amplitude Bj are the momenta, P(y β ), defined in Eq.(2.79). The legs of the Born process, on which the color and the helicity operators act, are clearly specified on the basis of the set {y}. The action of the color and the helicity operators at the square of the reduced Born process is illustrated in Fig. 3. In the DSA Figure 3: The structure of the square of the reduced Born process with the helicity-and color-correlation, Bj|T y emi · T yspe V y emi | Bj , is shown. The inner product of two color operators is denoted as T · T = c T c T c . the field species F({y}) of the reduced Born process Bj are fixed in one category Dipolej, namely, the identical set of the field species, is used for all the dipole terms which belong to the category, Dipolej. On the other hand the momenta of the Bj, P({y}) = {P(y a ), P(y b ); P(y 1 ), · · · , P(y n )}, are generally different functions of the original momenta, {p a , p b , p 1 , · · · , p n+1 }, associated with the set {x}, depending on the choice, (IJ, K).

Concrete formulae
In this part we specify the full expressions of the dipole terms. The expressions are separated into the cases that the field species of the emitter, F(y emi ), is a quark, or gluon.
In the case, F(y emi ) =quark, the expression in Eq.(2.81) is written as where the helicity correlation of the dipole splitting function, V y emi IJ,K , disappears and the function is fully factorized to the reduced Born amplitude. In the case, F(y emi ) =gluon, the expression is written as The formulae for the dipole terms in all categories in Eq.(2.24) are collected in Appendix A.1. Here we pick up two examples. The first example is the category, Dipole 1-(1)-1. In the category the dipole terms is where the dipole splitting function is given in Eq. (A.3) as The reduced momenta are given in Eqs. (A.30) and (A.31). The scalars, y ij,k and z i , appearing in the above formulae are defined in Eqs. (A.33) and (A.34). The second example is the category, Dipole 3-(6)-2, which is given in Eq. (A.24) as where the dipole splitting function is given in Eq. (A.25) as (2.87) The reduced momenta,p ai andk L , are defined in Eqs. (A.45) and (A.46). The scalar x i,ab is defined in Eq. (A.47). Sometimes it is convenient that the gluon polarization vector in the matrix element is taken in the basis of the helicity eigenstate as the circular polarization vector. The formulae for the dipole terms with the helicity correlation in the helicity basis are given in [93,98].

Examples
In order to demonstrate the creation algorithm we take the dijet event, pp → 2 jets, and pick up one input real process as, The input defines the set {x} with the field species and the momenta as Following the creation order in Eq.(2.24) and in Fig. 1, we create the all dipole terms from the input R 1 by choosing the three legs (x I x J , x K ) as where the notation (IJ, K) is the abbreviation of (x I x J , x K ). We pick up four dipole terms and write down the concrete expressions.
• Example 1: Dipole1 (1)-1: (13,2) The set {x} is defined with the field species and the momenta as A field mapping, y β = f (x α ), is found as The field mapping is unique in this case. The field mapping is interpreted as the identification between the elements in the sets, {x} and {y}, as Since the field species of the set {y} are fixed in the category Dipole 1, the expression in Eq. (2.100) can be abbreviated without confusion as (a, b ; 13, 2). (2.101) The momenta P({y}) are determined as Then the dipole term is given in Eq. (A.2) and is written down as where the dipole splitting function, V 13,2 , is given in Eq. (A.3) and the reduced momenta are determined in Eq.(2.102). The legs on which two color insertion operators act are clearly specified by referring to the elements y 1 and y 2 of the set {y}. It is again noted that in the same category Dipolej, in the present case, Dipole1, the same field species, F({y}), is used. Under the agreements we can abbreviate the expression in Eq.(2.103) as with the definition 1, 2 = B1 = {y} |T y 1 · T y 2 | B1 . While the indices, 1, 2, and 3, at the quantities, s 13 and V 13,2 , refer to the momenta of the original legs of R 1 , {p a , p b ; p 1 , p 2 , p 3 }, the arguments '1,2' inside the bracket, , refer to the elements of the set {y}.
• Example 2: Dipole1 (1)-1: (23,1) The set {x} is defined as The definition of the set {y} is same with Eqs. (2.96), (2.97), and (2.98). The field mapping is found as The mapping determines the identification of the elements as which can be abbreviated as (a, b; 1, 23). The momenta of the set {y} is determined as The dipole term is written as The set {x} is defined with the field species and the momenta as The reduced Born process is fixed as B2 = uū → gg, which determines the set {y} with the field species and the momenta as In this case the field mapping has two possibilities due to two identical fields in the final state, two gluons, as Both field mappings are equally allowed. To write down the concrete expression either of two possibilities must be chosen. Here we choose the first case in Eq. (2.118).
The identification of the elements are written as which is abbreviated as ( a, b ; 12, 3). The momenta is determined as (2.120) Referring to Eq. (A.20) the dipole term is written as Keeping in mind the reduced Born process, B2u, the expression is abbreviated as (2.122) • Example 4: Dipole3u (6)-2: (a1,b) The set {x} is defined as Next we fix the reduced Born process as The field mapping is uniquely found as which is interpreted as the identification as The expression is abbreviated as ( a1, b ; 2, 3). The momenta are also determined as The dipole term is written as which is abbreviated as For this dipole term we have another possibility to fix the reduced Born process as In this case the field mapping is found as This is abbreviated as ( b, a1 ; 2, 3). The momenta are determined as which is abbreviated keeping B3u ′ in mind as At the first case of the B3u, the element, x a1 , is identified with the element, y a , and at the second case of B3u ′ , the, x a1 , is identified with the y b . The first case appears to be a simpler expression in the sense that the leg a of the input process R 1 is connected to the leg a of the reduced Born process, B3u. In this reason the first case may be favored than the second case.

Summary
The hadronic cross section of a real correction subtracted by the dipole terms is written which is a part in Eq.(2.4). The partonic cross section is written aŝ It is sufficient that the real correction, |M(R i )| 2 4 , and the dipole term, D(R i ), are obtained in 4-dimension. The real correction |M(R i )| 2 4 is summed and averaged over both of the spin and color. The dipole term, D(R i ), is the summation of all the dipole terms under the input R i and is separated into the subcategories as (2.146) In each category, dipolej, the reduced Born process, Bj, is fixed, which defines the set {y} with the field species as (2.147) Once the reduced Born process and the set {y} are fixed, the necessary information to specify each dipole is the three elements of the set {x} and the field mapping as The form to specify the information can be abbreviated as 2. (a, b; 1, · · · , IJ, · · · , K, · · · , n + 1) . (2.150) Using the notation, y emi = f (x IJ ) and y spe = f (x K ), each dipole term is simply written down as which is abbreviated as The square of the reduced Born amplitude with the color-and spin-correlations, Bj |T · T V| Bj , is summed and averaged over the color degree of freedom. It is also summed over the spin configurations, but not averaged. Instead the dipole terms are divided by the spin average factor of the input real process as shown in Eq. (2.145). It is also noted that the symmetric factor by which the dipole terms are divided is not the symmetric factor of the reduced Born processes, S B j , but the symmetric factor of the input process, S R i , as also shown in Eq. (2.145).

I term creation
In this section 'Step 3. I (R i )' is explained. The input and the output of this step are: The B1(R i ) is the reduced Born process of the category Dipole 1, which is made from the input process R i on the rule, B1 = R i − g f , shown in Eq. (2.28). The creation algorithm given in the original paper [27] is to choose all combinations of two elements (y I , y K ) without duplicate from the set {y} = {y a , y b ; y 1 , · · · , y n }, of the process B1 (R i ), where the elements to be chosen are quark or gluon as F(y I/K ) = quark or gluon. In the DSA the creation algorithm of the I-terms is given as  Figure 4: The order to choose the first element y I is shown.

Creation order
In the DSA the order to choose the first element y I is determined as ( where the f f in/ini represents a quark in the final/initial state and the g f in/ini represents a gluon in the final/initial state as shown in Fig. 4. Each choice of the y I is followed by the choice of the second element y K . The choice of the y K in the final state is first and the choice in the initial state is second, which are denoted as the subcategories, 1 and 2, respectively. Then the creation order is written as Each pair (y I , y K ) specifies each I-term which is denoted as I (R i ) I, K . The summation of all created I-terms is the output, I (R i ), which is written as When all (n+2)-legs of B1, equally all (n+2)-elements of the set {y}, are quark or gluon, the indices, I and K, run over I, K = a, b, 1, · · · , n, and the total number of the I-term is (n + 2)(n + 1).

Concrete formulae
The concrete expression of each I term is given in universal form as where the common factor A d is defined as with the free parameter µ introduced in the dimensional regularization at d = 4 − 2ǫ.
The definition of the Casimir operator T 2 F(y I ) is same in Eq. (2.81). The universal singular function V F(y I ) is given as (2.158) where the n-body PS integration is defined in Eq. (2.12) as (2.162) The square of the color correlated Born amplitude, B1|T y I · T y K | B1 , is obtained in ddimension, and summed and averaged over the spin and color including the spin average factor, which is contrast to the case of the dipole terms in Eq. (2.6). Again in contrast to 'Step 2. D(R i )', in the present 'Step 3 ' the only one set {y} and the only one process B1 appear. Then we can drop the subscript y and the specification of the B1 from Eq.(2.156) without confusion as (2.163) where we introduce the notation for convenience as (2.164)

Complete set
Each I-term, I I, K , includes the square of a reduced Born amplitude with color correlations as T I · T K . The summation of I-terms, I = (I,K) I I, K , includes the square of a Born amplitude with all combination of the pairs (I, K). We call the set of the elements, T I · T K , with all the combinations of (I, K) as 'complete set of the square of the color correlated Born amplitude B1'. The name is sometimes abbreviated as 'complete set of the B1', which is explicitly written down as a, b , a, 1 , a, 2 , · · · , a, n , b, a , b, 1 , b, 2 , · · · , b, n , 1, a , 1, b , 1, 2 , · · · , 1, n , 2, a , 2, b , 2, 1 , · · · , 2, n , · · · , n, a , n, b , n, 1 , · · · , n, n − 1 } , where all the legs of the reduced Born process B1 are assumed to be quark or gluon. The number of the elements is (n + 2)(n + 1), which is same with the number of the I-terms.
The complete set of B1 is always included in the dipole terms in the category Dipole 1, because the dipole terms, D(R i : dip1), include the square of the reduced Born amplitude B1 with all the combination of the pair (y emi , y spe ) shown in Eq. (2.81) as B1 |T y emi · T yspe | B1 . Once we obtain the analytical or numerical expressions of the complete set of B1 as a function of the arbitrary input momenta, {P(y a ), P(y b ); P(y 1 ), · · · , P(y n )}, for the calculation of the dipole terms, the expressions can be reused for the calculation of the I-terms as well. Such a reuse of the expressions can save a certain amount of work to construct the subtraction terms.

Examples
We show some examples at the same process, R 1 = uū → uūg, in Eq. (2.88). The input for this 'Step 3. The concrete expression for the I-term, I 1, 2 , for example, is written as · V F(y 1 ) · s −ǫ 12 · B1 = {y a , y b ; y 1 , · · · , y n } |T y 1 · T y 2 | B1 , which is abbreviated as (2.171) The input momenta into the Born amplitude B1, {P(y a ), P(y b ); P(y 1 ), P(y 2 )}, are the momenta at the 2-body phase space in Eq.(2.162). The output, I(R 1 ), is simply written as I(R 1 ) = (I,K) (2.172) The I(R 1 ) is calculated in term of the complete set of the B1 = uū → uū, { B1|T I · T K | B1 } comp . The number of the elements of the complete set is twelve, which exactly correspond to the twelve ones, B1 |T y emi · T yspe | B1 , included in the dipole terms in the category Dipole 1 in Eq. (2.92).

Summary
The contributions of the virtual correction and the I-term to the hadronic cross section are written as which is in a part of Eq. (2.4). The partonic cross section is written aŝ where the virtual correction, |M virt (B1)| 2 d , is obtained in d-dimension, and summed and averaged over the spin and color. The I-term, I (R i ), is also obtained in d-dimension, and summed and averaged over the spin and color including the spin average factor. The output, I (R i ), is summation of all the created I-terms as (2.175) Once we determine the reduced Born process, B1, and the associated set {y}, each I term, I (R i ) I, K , is specified by the information of the pair, (I, K) . (2.176) The concrete expression of each I term is written in universal form as After the cancellation of the poles the PS integration is carried in the 4-dimension to be finite asσ (2.178)

P and K terms creation
In this section ' Step 4. P(R i ) and K(R i ) ' is explained. The input and the output are: Input: R i and Bj (R i ) , Output: P(R i ) and K(R i ) .
The symbol, Bj(R i ), represents the reduced Born process of the category Dipole j, which is made from the input process R i on the rules in Eq. 3. If F(x ai ) = F(y a ), create the pairs (y a , y K ) with K = 0, 1, · · · , n, b.
4. Write down the concrete expressions of all the P and K terms.
The steps 1, 2, and 3 are explained in the next part 'Creation order'. At the step 4 the concrete formulae for the P-term and the K-term are shown in the following parts, 'Concrete formula for the P-term' and 'Concrete formula for the K-term', respectively. Some examples are shown at 'Examples'. Finally we will have 'Summary'. The formulae for the P and K terms are collected in Appendix. A.3.
t Figure 5: The order to choose the pair (x a , x i ) is shown.

Creation order
We take the set {x} of the process R i in Eq. (2.21). The creation order is divided into the cases with the leg-a (x a ) and with the leg-b (x b ). We start with the case with lega. We choose the possible pairs (x a , x i ) from the set {x} in the order of the splittings, If F(x ai ) = F(y a ) → Create pairs : (y a , y K ) for K = 0, 1, 2, · · · , n, b , (2.179) where each pair creates each P and K terms as P(R i , x a : Bj, y a , y K ) for K = 1, 2, · · · , n, b , K(R i , x a : Bj, y a , y K ) for K = 0, 1, 2, · · · , n, b , (2.180) or If F(x ai ) = F(y b ) → Create pairs : (y b , y K ) for K = 0, 1, 2, · · · , n, a , (2.181) where each pair creates the P and K terms as P(R i , x a : Bj, y b , y K ) for K = 1, 2, · · · , n, a , K(R i , x a : Bj, y b , y K ) for K = 0, 1, 2, · · · , n, a . (2.182) Among the elements y K with K = 1, 2, · · · , n, a, and b, only the colored fields are taken as F(y K ) = quark or gluon. For convenience we categorize the P and K terms by the kinds of the second leg, y K . The P and K terms with y K = y 0 , y k for k = 1, · · · , n, and y a/b , are respectively categorized into 0,1, and 2. It is noted that when the final state of the R i includes the identical fields, only one pair of (x a , x i ) is taken and the others must be discarded. For example we take the process, R i = uū → ggg, and the set {x} = {x a , x b ; x 1 , x 2 , x 3 }. From this input we can find three pairs of the splitting (3) as (x a , x 1 ), (x a , x 2 ), and (x a , x 3 ). (2.183) Among three pairs we can take only one pair, for instance, (x a , x 1 ), and must discard the other two pairs, (x a , x 2 ) and (x a , x 3 ). The discard rule is contrast to the 'Step 2', where all three pairs must be taken for the creation of the dipole terms. The creation order with the leg-b is completely analog to the case with the leg-a.

Concrete formula for P term
The concrete formula for the P term with the leg-a (x a ), in the case of F(x ai ) = F(y a ) is given in universal form as . The Lorentz scalar, s xay K , is defined as s xay K = 2 p a · P(y K ) with p a in Eq. (2.23). The square of the reduced Born amplitude with color correlation, Bj|T ya · T y K | Bj , is obtained in 4-dimension, which is the same function of the momenta, P({y}), in the dipole term in Eq. (2.82), except for the spin average factor.
The squared amplitude, T · T , in the P term is summed and averaged over the spin and color including the spin average factor. The input momenta into the Born amplitude are written as P({y}) = {P(y a ), P(y b ); P(y 1 ), · · · , P(y n )}, (2.185) which are defined in the contribution to the cross section aŝ σ P (R i , x a : Bj, y a , y K ) = 1 0 dx 1 S B j Φ a (P(y a ), P(y b ) → P(y 1 ), · · · , P(y n )) 4 · P(R i , x a : Bj, y a , y K ). (2.186) The n-body phase space including the flux factor is defined as Φ a (P(y a ), P(y b ) → P(y 1 ), · · · , P(y n )) 4 = 1 F (P(y a ), P(y b )) with the initial momenta, (P(y a ), P(y b )) = (xp a , p b ) and the energy E y i = P(y i ) µ=0 . The phase space is identical with one in Eq. (2.13) as with the identification, p i = P(y i ) for i = 1, 2, · · · , and n. The expression for the P term in the case F(x ai ) = F(y b ) is similarly given as The phase space in this case, Φ a (P(y a ), P(y b ) → P(y 1 ), · · · , P(y n )) 4 , is the same expression in Eq. (2.187) and the initial momenta is defined as (P(y a ), P(y b )) = (p b , xp a ). The concrete formula for the P term with the leg-b, P(R i , x b : Bj, y b/a , y K ), and the phase space, Φ b (P(y a ), P(y b ) → P(y 1 ), · · · , P(y n )) 4 , are completely analog to the case with the leg-a. We define the output P(R i ) as the set which consists of all the created P-terms, where the element, P(R i , x a ), is defined as the subset, P(R i , x a ) = { P(R i , x a : B1), P(R i , x a : B3), P(R i , x a : B4)} . (2.191) The elements, P(R i , x a : Bj), are the summation over the P terms including Bj as P(R i , x a : Bj) = n K=1 P(R i , x a : Bj, y b/a , y K ) + P(R i , x a : Bj, y b/a , y a/b ) . (2.192) In the case with the leg-b, the set P(R i , x b ) and the summation, P(R i , x b : Bj), are similarly defined.

Concrete formula for K term
The concrete formula for the K term with the leg-a (x a ), in the case of F(x ai ) = F(y a ), are separated into three categories, 0, 1, and 2, introduced above. The formulae are written as K(R i , x a : Bj, y a , y 0 ) = α s 2πK F(xa)F(ya) (x) · Bj |Bj , (2.193) K(R i , x a : B1, y a , y k ) = α s 2π (2.195) The symbol, Bj |Bj , is the abbreviation of Bj = {y a , y b ; y 1 , · · · , y n } | Bj , which is the usual squared amplitudes at the LO process. The symbol, Bj|T ya · T y K | Bj , is the abbreviations of Bj = {y a , y b ; y 1 , · · · , y n } |T ya · T y K | Bj , which is the same quantity in Eq. (2.184). The functions of x,K F(xa)F(ya) (x), h(x), and K F(xa)F(ya) (x), and the symbol, γ F(y k ) , are defined in Appendix A.3. It is noted that the K terms with pair (y a , y k ) with k = 1, 2, · · · , and n, exist only in the case of the diagonal splittings, namely, the case including the process B1 shown in Eq. (2.194).
Same with the P term the input momenta into the Born amplitude, P({y}), are given in the contribution to the cross section aŝ where the n-body PS, Φ a , is same in Eq. (2.187). The formulae in the case of F(x ai ) = F(y b ) are similarly given as (2.199) The formulae for the K term with the leg-b, K(R i , x b : Bj, y b/a , y K ), and the phase space, Φ b , are completely analog to the case with the leg-a.
Again similar to the case of the P-term we define the output K(R i ) as the set which consists of all the created K-terms, where the element, K(R i , x a ), is defined as the subset,

Examples
To demonstrate the creation of the P and K terms we take the same input process for the dipole terms creation in Eq. (2.88) as The sixteen pairs are created and each pair corresponds to a P and a K terms. There are two kinds of exceptions noted above already. The first exception is the pairs of the type, (y a/b , y 0 ), in this example, 1, 5, 9, and 13, which produces only a K-term. The second exception is the following. As noted the K terms with the pair (y a/b , y k ) with k = 1, 2, · · · , and n, exist only for the diagonal splittings. Then the K terms with the pairs with the non-diagonal splittings, in this example, the pairs, 6, 7, 14, and 15, do not exist. We show the concrete expressions specified from three pairs, 1.(y a , y 0 ), 2.(y a , y 1 ) and 16.(y b , y a ), for examples.

Summary
The contributions of the P and K terms to the hadronic cross section are written as (2.215) The partonic cross sections are written aŝ where the PS, Φ a , is given in Eq. (2.187). The PS integrations of the P and K terms are separately finite in 4-dimension. The outputs, P(R i ) and K(R i ), are the sets defined in Eqs. (2.190) and (2.200), respectively. Once we select an input process R i , a leg-a or -b, and a reduced Born process, Bj, each P and K terms, P/K(R i , x a/b : Bj, y a/b , y K ), are specified by the information of the pair, (y a/b , y K ), which is abbreviated as (a/b, K) . (2.217) The concrete formulae for the P and K terms are collected in Appendix A.3.

Advantages of the DSA
This formula shows that the real correction,σ R (R i ), the virtual correction,σ V (B1(R i )), and all the subtraction terms,σ D (R i ),σ I (R i ),σ P (R i ), andσ K (R i ), which are created from an input process, R i , have the same initial parton states, F(x a ) and F(x b ), and they are all multiplied by the same PDFs, . In other words the subtraction terms are sorted by the initial parton states. This is the first feature. The second feature is that the subtraction terms are sorted by the reduced Born processes as well. As defined in the previous sections the creation order of the D, I, P, and K terms are sorted by the kinds of the splittings and the reduced Born processes, Bj, with j = 1, 2, 3, and 4, where the processes, B2, B3, and B4, may have the subcategories about the quark flavors. The third feature is that we introduce the sets, {x},{x} and {y}, and the field mapping y = f (x). Using the sets and the mapping each subtraction term is specified in a well defined compact form.
The above mentioned three features of the DSA lead to the following three advantages of the DSA: 1. Consistency proof of the subtraction terms, 2. Easy construction of the codes for the Monte Carlo integration, 3. The subtraction terms and the summary tables in a compact form.
We start the explanation with the first advantage. By the construction of the dipole subtraction procedure the summation of all the introduced subtraction terms must vanish as written in Eq. (1. 3 (2.219) The first advantage is that the straight forward proof of the consistency of the subtraction terms in Eq. (2.219) is possible. The cancellation among the subtraction terms can be realized among the subtraction terms which have the same initial parton states and the same reduced Born processes. According to the first two features it is clear that the cancellation is realized among which categories of the subtraction terms and then a systematic proof of the consistency becomes possible in the DSA. We can construct a straightforward algorithm to prove Eq. (2.219), which is presented at the accompany paper. The second advantage is the following. In order to construct the computer codes for the Monte Carlo integration, we must collect the subtraction terms which have the same initial parton states to be multiplied by the same PDFs. This collection is realized in the DSA thanks to the first feature that the created subtraction terms are sorted by the initial parton states. The third advantage is the followings. According to third feature we can specify all the subtraction terms in a compact form. For example, the compact forms for the D, I, P, and K terms are shown in Eqs. of the P and K terms is different from the DSA. In order to demonstrate the difference we take the same example process, R 1 = uū → uūg, in Eq. (2.203). The creation algorithm of the P and K terms in the AutoDipole takes the reduced Born process, B1 = R 1 − g f = uū → uū, as the input. Then the P and K terms are created by adding to the process B1 the possible splittings in Fig. 5. The splitting (3) can be added to the y a of the B1 and the elements y K are chosen. The choice creates the P and K terms as P/K(R 1 , x a : B1, y a , y K ), which are same with the case of the DSA. As next choise the splitting (7) can be added and the y K are chosen. The choice creates the P and K terms, which are written in the notation defined in Eq. (2.180) as P/K(R i = ug → uūu, x a : B4u, y a , y K ). (2.220) As the notation shows, these P and K terms are created from the input R i = ug → uūu, when the splitting (7) is applied. In this way the creation place of the non-diagonal splittings (6) and (7), are different between the AutoDipole algorithm and the DSA. The advantage of the AutoDipole algorithm can be that all the P and K terms include only one kind of the reduced Born process B1. In this sense the collected expressions of the P and K -terms, which are created form the input R i , are simpler than the case of the DSA. The disadvantage of the AutoDipole alogorithm may be that the P and K terms with non-diagonal splittings, created by the AutoDipole algorithm, have the different initial parton states from the other subtraction terms. The involvement of the different initial states spoils the first feature of the DSA and hence lose the the first and second advantages of the DSA. Namely in the AutoDipole algorithm the proof of the consistency of the subtraction terms becomes complex and for the Monte Carlo integration the recollection of the non-diagonal P and K terms is required as extra work of the users. The third advantage of the DSA, the expressions and the summary tables in a compact form, is held also for the AutoDipole algorithm, because the subtraction terms are sorted by the reduced Born process in the AutoDipole algorithm as well.

List of R i
At Step 1 we make the list of the contributing real correction processes, {R i }, as shown in Eq. (2.2), There are three independent processes as n real = 3 . The number of the final states is (n + 1) = 3. At Drell-Yan process in order to exhaust all independent partonic processes, it is sufficient to take into account for one quark flavor, u, for instance. Finally we can generalize to the full five massless quark flavors.

D(R 1 ) creation
The input process, R 1 = uū → µ − µ + g, determines the set {x} with the field species and the momenta as We create the dipole terms in the order shown in Fig .1 as Then we specify the field mapping for each dipole term and write down the concrete expression.

(a3,b)
The set {x} is defined with the field species and the momenta as where the reduced momenta,p a3 andk 1/2 , are defined in Eqs. (A.45) and (A.46). Then we construct the field mapping as x 2 ) = (y a , y b ; y 1 , y 2 ), (3.12) which is interpreted as the identification of the elements as (y a , y b ; y 1 , The expression is abbreviated as ( a3, b ; 1, 2). The field mapping determines the momenta as The dipole term is written in Eq. (A.12) as

(b3,a)
The set {x} is defined as The set {y} is fixed in Eqs. (3.6) and the field mapping is found as (y a , y b ; y 1 , y 2 ) = (x a , x b3 ; x 1 , x 2 ) , (3.19) which is abbreviated as ( a, b3 ; 1, 2). The momenta are determined as The output, D(R 1 ), is written as The input process, R 2 = ug → µ − µ + u , determines the set {x} as We create the dipole term as Dipole4 (7) -2 : 1. (b3, a) . which is abbreviated as ( a, b3 ; 1, 2). The momenta are determined as The dipole term is written in Eq. (A.28) as The output, D(R 2 ), is written as Similar to the creation of D(R 2 ), the dipole term is created in the following way. Input process: Field mapping: ( a, b3 ; 1, 2) .

Summary of creation
The created dipole terms from the inputs, {R 1 , R 2 , R 3 }, are summarized at Table. 1.

I term
At Step 3 we create the I term, I (R i ), from the input, B1 (R i ). Among the real correction processes, {R 1 , R 2 , R 3 }, only the R 1 has the reduced Born process B1 as The set {y} is fixed in Eqs. The concrete expressions are given in Eq. (A.50) as The output I (R 1 ), is written as

P and K terms
At Step 4 we create the P and K terms, P(R i ) and K(R i ), from the inputs, R i and B j (R i ).

P/K(R 1 ) creation
The process R 1 defines the set {x} in Eq. The possible reduced Born process is only the B1(R 1 ) as  4.(y b , y a ) . 2.(y a , y b ) P(R 1 , x a : B1, y a , y b ) = α s 2π and with the leg-b as The output for the P term is written in Eq. (2.190) as where the elements, P(R 1 , x a/b ), are written as P(R 1 , x a/b ) = P(R 1 , x a/b : B1) = P(R 1 , x a/b : B1, y a/b , y b/a ) . (3.67) The output for the K term is written in Eq. (2.200) as where the elements, K(R 1 , x a/b ), are written as The concrete expressions are written down as The outputs are written as and (3.75)

P/K(R 3 ) creation
Similar to the P/K(R 2 ) creation, the P and K terms are created in the following way.

Summary of creation
The created P and K terms from the inputs of the real processes, {R 1 , R 2 , R 3 }, and the reduced Born processes are summarized at Table. 3.

At
Step 5 we obtain the NLO correction, σ NLO in Eq. (2.3) as The σ(R 1 ) is concretely written as where the finite combinations of the partonic cross sections are separately written aŝ The subtraction terms, D(R 1 ), I(R 1 ), P(R 1 , x a/b ), and K(R 1 , x a/b ), are given in Eqs.
where the finite cross sections are written aŝ where the contributions are written aŝ

List of R i
At Step 1 we make the list of the real processes {R i } as There are eleven independent processes as n real = 11 . The number of the final states is (n + 1) = 3. We assume the five massless quark flavors, u, d, s, c, and b. The contributing real processes can be exhausted by the independent processes which are produced from the field contents of only two quark flavors and a gluon. We take the u and d quarks as the two flavors. The real processes in the curly bracket in Eq. (4.1) represent the processes which are obtained by the replacements of the quark flavors. For example the process, R 1d , is concretely written as R 1d = dd → ddg, which is obtained by the replacements, u → d andū →d at the process, R 1u = uū → uūg. It is sometimes useful to categorize the partonic processes by the crossing symmetry. The processes, R 1u , R 2u and R 3u , are categorized into the master process, 0 → uuūūg. The processes, R 4u , R 5ud , R 6ud , and R 7u , are categorized into the process, 0 → udūdg. The processes, R 8u , R 9u , and R 10u , are into 0 → uūggg. The process, R 11 , is into 0 → ggggg.

D term
At Step 2 we create the dipole terms, D (R i ), from the inputs, {R i }, in Eq. We pick up the dipole term, 16. D (R 9 : dip1-(4)-1) b2,1 , among twenty seven dipole terms in Table 12. For the process R 9 , the reduced Born process B1 is fixed as B1(R 9 ) = ug → ug , For the dipole term, D b2,1 , the set {x} is defined as mapping is specified as (y a , y b ; y 1 , y 2 ) = (x a , x b2 ; x 1 , x 3 ) , (4.11) which is abbreviated as (a, b2; 1, 3) in Table 12. The field mapping determines the momenta of the set {y} as P({y}) = {p a ,p b2 ;p 1 , p 3 } . (4.12) The concrete expression of the dipole term is given in Eq. (A.14) as where the dipole splitting function, V b2,1 , and the Lorentz scalar, Eq. (2.6) asσ D (R 9 ) = 1 S R 9 Φ(R 9 ) 4 1 n s (u)n s (g) D (R 9 : dip1-(4)-1) b2,1 + · · · . (4.14) It is noted that the symmetric factor of the reduced Born process, B1(R 9 ), is S B 1 = 1, but the contribution of the dipole term is not divided by the factor. The contribution of all the dipole terms involved in D (R 9 ), must be divided by the symmetric factor of the input real process R 9 , S R 9 = 2.

I term
At Step 3 we create the I terms, I (R i ), from the inputs, B1 (R i ). The summary tables   of all the created I terms are shown at Tables 15 ∼ 19 at Appendix B.2. Since the real processes, R 3 and R 7 , do not have the reduced Born process B1, the I terms, I(R 3 ) and I(R 7 ), do not exist. The details of the creation of I (R 1 ) have been explained from Eq. (2.166). Here we see the I (R 9 ) concretely. The input to create I (R 9 ) is the process, B1(R 9 ) in Eq. (4.5) and the associated set {y} is defined in Eqs. (4.6) and (4.7). Twelve I terms are created and they are listed at Table 17. Referring to the formula in Eq. (A.50) the concrete expression of I (R 9 ) is written as where the V f /g are defined in Eqs. (A.52) and (A.53). The contribution to the cross section reads in Eq. (2.8) asσ where the cross section is divided by the symmetric factor of the B1(R 9 ), S B 1 = 1.

P and K terms
At Step 4 we create the P and K terms, P(R i ) and K(R i ). The summary tables of all the created P and K terms are shown in Tables 20 ∼ 30 at Appendix B.3. The details of the creation of P/K (R 1 ) are presented with some examples from Eq. (2.203). In this section we show only one P term in P (R 9 ) and one K term in K (R 9 ). The input R 9 defines the set {x} in Eqs. (4.2), (4.3), and (4.4). The possible reduced Born processes B j (R 9 ) are fixed during the creation of the dipole term D (R 9 ), which are explicitly shown in Table 12  Here we show the concrete expressions for the P and K terms, 10. P/K(R 9 , x b : B1, y b , y 1 ), in (4.22) The contributions of the P and K terms, P/K(R 9 , x b : B1, y b , y 1 ), are divided by the symmetric factor of the reduced Born process B1, S B 1 = 1. For comparison, the contributions from the terms, P/K(R 9 , x b : B 4u ), are divided by the symmetric factor of the reduced Born process B 4u , S B 4u = 2. In this way the P and K terms are divided by the symmetric factor of the reduced Born processes, S B j , not by the symmetric factor of the input real process, S R i . In this sense the use of the symmetric factor for the P and K terms are the inverse manner to the use for the dipole terms, which is explained at the end of the Sec. 4.2.

At
Step 5 we obtain the NLO correction, σ NLO , in Eq. (2.3) as where the summation over the different quark flavors is suppressed. The NLO cross section, σ(R i ), is given in the formula in Eq. (2.4) as (4. 24) In the cases of the NLO cross sections, σ(R 3 ) and σ(R 7 ), the formula is simplified as

Analytical check at Drell-Yan
In this section we have an analytical check at the Drell-Yan process. In Sec. 5.1 we review the well known analytical results obtained by the traditional method. In Sec. 5.2 we obtain the analytical results by the DSA. We will show that both of the results exactly coincide.

Traditional method
In this section we review the well known results, which were obtained for the first time at the pioneer works, [11][12][13][14][15][16][17][18][19]. The method used in the works became the traditional method to calculate the QCD NLO correction at the hadron collider processes since then. In the method both of the real and virtual corrections are calculated in the d-dimension, that is, the matrix elements, the PS integrations, and the spin-average factors, are all defined in the d-dimension. In the method not only the UV, soft, and collinear divergences in the virtual correction, but also the soft and collinear divergences in the real correction are regularized as the poles, 1/ǫ or 1/ǫ 2 . We start with the review of the LO contribution. The general formulae are given in Eqs. (2.17) and (2.18). In this method we redefine the partonic cross section in Eq. (2. 19) in the d-dimension. At the Drell-Yan process we assume one quark flavor, u, and finally generalize to five massless flavors. At LO only one independent process exists as The Feynman diagram is shown in Fig. 6 2 . The contribution is written as where the partonic cross section is defined in the d-dimension aŝ We fix the kinematical valuables as follows. The squared energy of two protons in the initial state is denoted as S = (P a + P b ) 2 = 2P a · P b with the momenta, P a/b . The squared energy of two partons in the initial state is denoted asŝ where the momenta are denoted as p a/b = x 1/2 P a/b . The square of the muon pair invariant mass is denoted as M 2 µ + µ − = q 2 = (p µ + + p µ − ) 2 with the momenta of the anti-muon/muon, p µ + /µ − . We write the total cross section asσ tot (L 1 ) =σ LO (q 2 ) and calculate it aŝ where the color degree of freedom of quark is denoted as N c = 3, and the symbol µ represents a dimensionful free parameter introduced in the d-dimensional space-time. The last factor, P(q 2 ), is defined as with the electric charges, (Q u , Q µ ) = (2/3, −1), and the constants from the Z boson coupling, (v u , v l , a u , a µ ) = (1/4−2s 2 /3, −1/4+s 2 , −1/4, 1/4). The square of the Weinberg angle and the Z boson mass are denoted as s 2 = sin 2 θ W and M z , respectively. We omit the contribution from the decay of the on-shell Z boson for simplicity. For reference the total cross section of the process, L 1d = dd → µ − µ + , is obtained by the replacements, (Q u , v u , a u ) → (Q d , v d , a d ) = (−1/3, −1/4 + s 2 /3, 1/4) . We predict the distribution of the squared muon pair invariant mass, M 2 µ + µ − = q 2 , because the observable may be the simplest and most typical one at the Drell-Yan event. The q 2 -distribution is calculated as where the variable z is defined as z = q 2 /ŝ. The NLO corrections are written for general processes as where the processes, {R i }, are the real correction processes, which are same ones listed at the Step 1 in the DSA. Each cross section, σ trad (R i ), is written as where the real correction,σ R (R i ), is defined in d-dimension, contrasted to the quantity in Eq. (2.5). The virtual correction,σ V (B1(R i )), is same in Eq. (2.7). The symbol,σ C (R i ), represents the so-called collinear subtraction term, which is concretely shown later. Figure 7: The diagram of the real process, R 1 = uū → µ − µ + g .
The Drell-Yan process has three independent real correction processes listed in Eq. (3.1) as Each cross section is written as 11) and the cross section, σ trad (R 3 ), is obtained by the replacement, u →ū, at σ trad (R 2 ). First we show the partonic cross sections included in the σ trad (R 1 ) in Eq. (5.10).
The Feynman diagram of the real correction is shown in Fig. 7. The q 2 -distribution is calculated as where the factor, C(q 2 , ǫ), is defined as Figure 8: The one-loop diagram of the process, B1( The a R is written as (5.15) The expression includes the so-called '+'-distribution, that is defined for arbitrary functions, g(z) and h(z), as with a value, 0 ≤ τ < 1.
•σ V (B1(R 1 )) The virtual correction,σ V (B1(R 1 )), is defined in Eq. (2.7) aŝ The Feynman diagram of the virtual one-loop correction is shown in Fig. 8. The q 2 -distribution is calculated as where the factor, C(q 2 , ǫ), is same in Eq. (5.14) and after the subtraction of the UV pole, 1/ǫ UV , by the renormalization program, the a V is obtained as The q 2 -distribution of the collinear subtraction term is written as where the q 2 -distribution of the LO process, L 1 , with the rescaled initial energy, xŝ, is defined as The symbol, µ F , represents the mass factorization scale. The four-dimensional Altarelli-Parisi splitting function, P uu (w) = P f f (w), is given in Eq. (A.58). There are two contributions to the collinear subtraction term. One is the contribution of the gluon radiation from the u-quark leg in the initial state at the process, R 1 = uū → µ − µ + g. The other is the contribution of the radiation fromū-quark leg at R 1 . Both contributions are written as the identical expressions and the expression in Eq. (5.20) includes the both contributions. The q 2 -distribution is calculated as where the factor, C(q 2 , ǫ), is same in Eq. (5.14) and the a C is written as At this stage with the finite results in 4-dimension we can return the common factor, C(q 2 , ǫ), also into 4-dimension as The real correction,σ R (R 2 ), is defined aŝ The Feynman diagram of the real correction is shown in Fig. 9. The q 2 -distribution is calculated as where the splitting function, P gf (z), is given in Eq. (A.61).
The q 2 -distribution of the collinear subtraction term is written as where the q 2 -distribution of the LO process is same in Eq. (5.21) and the splitting function, P gū (z) = P gf (z), is given in Eq. (A.61). The distribution is written in the form, The summation of the two contributions in Eqs. (5.29) and (5.32), is free from the divergences as Taking account of the exchanged initial states and all five massless flavors, we obtain the prediction of the q 2 -distribution at NLO accuracy as The contributions from the distributions, dσ(L 1 )/dq 2 and dσ trad (R 1 )/dq 2 , are written as where the combination of the PDFs are denoted as H qq (x 1 , x 2 ) = f q (x 1 )fq(x 2 )+fq(x 1 )f q (x 2 ). The valuable, x 2 , is expressed as x 2 = τ /(x 1 z) with τ = q 2 /S. The partonic distributions, dσ(L 1 )/dq 2 and dσ trad (R 1 )/dq 2 , are given in Eqs. (5.6) and (5.24), respectively. The contributions from the distributions, dσ trad (R 2 )/dq 2 and dσ trad (R 3 )/dq 2 , are written as with the PDFs, H qg (x 1 , x 2 ) = f q (x 1 )f g (x 2 ) + f g (x 1 )f q (x 2 ). The partonic distribution, dσ trad (R 2 )/dq 2 , is given in Eq. (5.34). The expression of the partonic distribution, dσ trad (R 3 ) /dq 2 , is identical to the distribution, dσ trad (R 2 )/dq 2 .

DSA
In this section we use the dipole subtraction procedure through the DSA to obtain the analytical results of the q 2 -distribution at the Drell-Yan process. The DSA has been already applied to the Drell-Yan process at Sec. 3. The results are summarized in Sec. 3.5. We start with the explicit calculation of the σ(R 1 ) in Eq. (3.83). The subtracted real cross section, (σ R (R 1 )−σ D (R 1 )), in Eq. (3.84) is defined in 4-dimension. But the analytical integration of the PS in 4-dimension seems not to be easy. Instead we redefine the cross section in d-dimension and regularize the soft and collinear singularities as the poles of 1/ǫ and 1/ǫ 2 , which are produced by the d-dimensional PS integration. The poles from theσ R (R 1 ) and theσ D (R 1 ) are canceled each other, and the finite analytical expression at 4-dimension is obtained. The distribution, dσ R (R 1 )/dq 2 , is calculated in Eqs. (5.13) and (5.15). Then we proceed to the calculation ofσ D (R 1 ), which is now defined in d-dimension asσ There are two dipole terms for the process R 1 , shown in Eq. (3.22). The first dipole term, D(dip1-(3)-2) a3,b , is written in Eq. (3.15), which is now defined in d-dimension. The dipole splitting function, V a3,b , in d-dimension reads in [27] as 40) and the square of the reduced Born process, B1 |T·T| B1 , is also defined in d-dimension.
The contribution to the partonic cross section of the dipole terms can be analytically integrated over the PS region including the soft and collinear singularities. The analytical integration of the dipole terms for arbitrary process is one core part for the construction of the dipole subtraction procedure [27]. The contribution of the first dipole term to the cross section in Eq. (5.39) is written asσ D (R 1 : D a3,b ), and the integration formula is applied to the cross section aŝ where the Lorenz scalar is denoted as s ab = 2p a · p b , and the function, V u,u (x; ǫ), reads in [27] as The correlation of two color insertion operators in the square of the Born process, B1, is fully factorized as B1 |T ya · T y b | B1 = −C F B1|B1 . The q 2 -distribution is written as where the distribution with the scaled initial energy, dσ(L 1 , xŝ)/dq 2 , is defined in Eq. (5.21). The contribution from the second dipole in Eq. (3.21) is written in the identical expression.
Then the contribution of all dipole terms in D(R 1 ) is written in the form, where the common factor C(q 2 , ǫ) is same in Eq. (5.14), and the a D is written as (5.45) which is expanded by the 1/ǫ as Recalling the expression of the a R in Eq. (5.15) the q 2 -distribution of the subtracted real correction is written as Next we calculate the subtracted virtual correction, (σ V (B1(R 1 )) +σ I (R 1 )), in Eq. (3.85). The concrete expression of the virtual correction,σ V (B1(R 1 )), is given in Eqs. (5.18) and (5.19). The expression of the I term, I(R 1 ), is given in Eq. (3.52) and the contribution to the partonic cross section is written aŝ The q 2 -distribution is calculated in the form, with Then we obtain the q 2 -distribution of the subtracted virtual correction as The cross section of the P and K terms is written in Eq. (3.86). The P and K terms are obtained in Eqs. (3.66) and (3.68), and the cross section is explicitly written down aŝ The q 2 -distribution is calculated as Then the q 2 -distribution, dσ(R 1 )/dq 2 , is written as where the finite quantities, (a R − a D ), (a V + a I ), and a PK , are given Eqs. (5.48), (5.54), and (5.58). The summation of three contributions is calculated as The results exactly coincide with the results obtained by the traditional methods in Eq. (5.25). Next we calculate the q 2 -distribution of σ(R 2 ) in Eq. (3.87). The subtracted real correction is written in Eq. (3.88). Similarly to the case of σ(R 1 ), we redefine the cross section in d-dimension. The distribution, dσ R (R 2 )/dq 2 , is obtained in Eqs. (5.29) and (5.30). Then we proceed to the calculation ofσ D (R 2 ), which is defined in d-dimension aŝ For the process R 2 , only one dipole exists in Eq. (3.35) and the expression is given in Eq. (3.34). The cross section is analytically integrated over the soft and collinear region, and the q 2 -distribution is calculated as where the function, V g,f (z; ǫ), reads in [27] as The result is written in the form, The difference between the a R,ug in Eq. (5.30) and a D,ug in Eq. (5.65) is calculated as The cross sections of the P and K terms are written in Eq. (3.89). The P and K terms, P(R 2 ) and K(R 2 ), are created in Eqs. (3.74) and (3.75). The cross section is explicitly written aŝ Then the q 2 -distribution is calculated in the form, Finally we obtain the q 2 -distribution, dσ(R 2 )/dq 2 , as The results exactly coincide with the results in Eq. (5.35). The calculation of the distribution, dσ(R 3 )/dq 2 , is completely analog to the dσ(R 2 )/dq 2 , and the results are the identical expression. In this way it is shown that the analytical results by the DSA exactly agrees with the well known results by the traditional methods as for i = 1, 2, and 3 .

Summary
We treat with the subject of the QCD NLO correction at hadron collider like LHC. At the simple processes the analytical results are obtained by the traditional method which may originate in the pioneer works for the Drell-Yan process, [11][12][13][14][15][16][17][18][19]. The traditional method is reviewed at Sec. 5.1. At the complex processes, namely, the multi-parton leg processes, which are typically the processes, pp → Standard model particles + n jets , it is almost impossible to obtain the analytical results for the NLO corrections. The dipole subtraction procedure overcomes some difficulties of the calculation and made it possible to obtain the NLO correction at the multi-parton leg processes. The price to employ the dipole subtraction is mainly two things: one is that many subtraction terms must be created and the expressions are not so simple. The other one is that the large amount of the calculation is executed as the numerical evaluation for the Monte Carlo integration.
As the consequence the person who has obtained the results of the NLO corrections, has the difficulty to confirm whether the obtained results are true or false. For the other person who does not calculate by himself, the confirmation is more difficult. In order to solve one part of the problem we wish a practical algorithm to use the dipole subtraction which satisfies the following requirements, -Well defined steps with the documented formulae, -Summary tables to specify all the created subtraction terms, -Associated algorithm to prove the consistency of the created subtraction terms.
In this article we present such a algorithm which satisfies all the above requirements. We name it the dipole splitting algorithm (DSA) and define it in Sec. 2. The master formulae and all the steps are defined in Sec. 2.1 and the formulae for the subtraction terms are given in Appendix. A. The creation algorithm and the concrete expressions for the D, I, and P/K terms are explained in the sections, 2.2, 2.3, and 2.4, respectively. The advantage of the DSA, including the realization of the above requirements, is clarified in Sec. 2.5.
We demonstrate the DSA at the Drell-Yan process in Sec. 3 and at the dijet process in Sec. 4. The summary tables for the Drell-Yan process are shown in Sec. 3. The summary tables for the dijet process are shown at the tables in Appendix. B. Those tables can be a template for the summary table to specify all the subtraction terms created by the DSA.
Regarding the use of the summary tables we intend that the other person who does not execute the DSA by himself can specify the created subtraction terms and write down the concrete expressions, only by reading the tables on a document without any direct communication with the author of the tables. The proof algorithm of the consistency of the subtraction terms will be presented in the accompany paper, [100]. As one most reliable confirmation of the DSA, we have the analytical check against the results by the traditional method at the Drell-Yan process in Sec. 5.2.
We plan the following subjects in future. We will apply the DSA to some processes at LHC and make the predictions at the NLO accuracy. The DSA can be applied to the processes at the e − e + -and e − p -colliders, as well. The application is easy, because the DSA becomes simplified for the processes than the case of the hadron collider presented in this article. The DSA is presented for arbitrary process including only the massless quarks in this article. We will present the DSA for the processes including massive quarks. The extension should be straightforward, because the construction algorithm of the subtraction terms for the massive quarks given in the article [28] are same with the case of the massless quarks given in the article [27]. Regarding the automation of the DSA as a computer package, the AutoDipole [93] is a good candidate to be implemented, because the creation algorithm of the dipole and I terms implemented in AutoDipole is essentially same with the DSA and only the creation algorithm of the P and K terms are different from the DSA. Then it is sufficient that only the part to create the P and K terms in the whole codes is modified. It is hoped that the DSA will serve to obtain the reliable predictions at the NLO accuracy.
A Formulae for DSA A.1 D term Bj |T y emi · T yspe V y emi IJ,K | Bj . (A.1)