ME-PS matching in the simulation of multi-jet production in hadron collisions using a subtraction method

The subtraction method for the matching between the matrix element (ME) and parton shower (PS), that has been developed for combining 0-jet and 1-jet production processes in association with electroweak-boson production in hadron collisions, is extended to multi-jet production. In order to include multi-jet MEs, we have to address the soft-gluon divergence together with the collinear divergence. We introduce an approximation which simultaneously reproduces both divergences in a form suitable for application to our subtraction method. The alteration in the subtraction can be compensated by applying an appropriate correction to corresponding non-radiative events. We demonstrate that $W$ + 0, 1, and 2 jet production processes can be consistently combined using the developed matching method.


Introduction
As a result of higher-order quantum chromodynamics (QCD) interactions, hard interactions producing large transverse momentum particles or heavy particles, e.g., weak bosons, are frequently associated with hadron jets in high-energy hadron collisions such as the proton-proton collisions provided by CERN LHC. A precise understanding of the jet production is necessary not only to achieve precise measurements of hard interactions but also to reliably estimate the background of unknown new phenomena. Since the jet production and other soft QCD interactions complicate the produced events, simulations based on Monte Carlo (MC) event generators are indispensable tools for measurements at hadron colliders. Hence, it is important that the jet production is precisely reproduced in MC event generators. The hadron jets are considered to originate from the production of energetic partons (light quarks and gluons) that can be evaluated in the framework of perturbative QCD (pQCD). Multi-jet production can be simulated by parton showers (PSs) in MC event generators. However, since the PS is based on the collinear approximation of pQCD, we cannot expect satisfactory precision for the production of isolated energetic jets that mimic the signature of various new phenomena. Simulations based on the matrix element (ME) including the jet (energetic parton) production are necessary for reproducing such phenomena.
We also have to apply PS simulations to events generated according to jet-including MEs in order to provide realistic simulations. A problem arises when we simultaneously employ ME and PS for the jet production. There may be an overlap (double-count) between the two simulations. The problem PTEP 2015, 053B04 S. Odaka et al. subsequently demonstrate that the developed matching method provides a reasonable simulation of W + jet production up to two jets.
The rest of this article is organized as follows. We introduce a subtraction method that simultaneously subtracts the collinear divergence and the soft-gluon divergence in Sect. 2. The correction to PS-applied events that compensates for the alteration in the subtraction is described in Sect. 3. The performance of the ME-PS matching method employing the developed techniques is tested for W production processes in Sect. 4, and the discussions are concluded in Sect. 5.

Combined subtraction
The collinear approximation that is subtracted from the squared MEs of radiative processes in the LLL subtraction can be expressed as [15] T coll = a T coll a , T coll where the sum is taken over all particles in the initial and final states (external legs) that can emit the considered radiation. We call the radiating external leg a the emitter. The parameter α s is the QCD coupling, and Q 2 a is the squared momentum transfer of the considered branch for the initial-state radiation (ISR) and the squared invariant mass of the branch products for the final-state radiation (FSR). The factor P a (z a ) represents the corresponding Altarelli-Parisi splitting function [25]. Although the exact definition depends on the kinematics model [15], the variable z a can be approximately considered as the energy fraction carried by the child of the branch, i.e., the branch product other than the radiation. Hence, the radiation carries the energy fraction of approximately 1 − z a . The factor T 0,a represents the squared ME of the non-radiative event that is reconstructed by removing the considered radiation from the emitter a, by exactly reversing the kinematics model of PS branches. Finally, the approximation is multiplied by the ratio between the squared center-of-mass (CM) energies of the radiative and non-radiative events,ŝ/ŝ 0,a [17], in order to correct for the change of the flux factor in the cross-section calculation; usually,ŝ/ŝ 0,a = 1/z a for ISR.
We have demonstrated in previous studies that the cross sections can be made finite by subtracting the approximation in Eq. (1) from 1-jet MEs [12,15]. However, the subtraction is not successful for 2-jet processes based on 2-jet production MEs. Figure 1 shows the p T distribution of the smaller p T parton in the qg → W q g reaction after the subtraction, where p T is measured with respect to the beam direction. The quark combination, qq , represents all possible pairs of up-type and down-type quarks up to the b quark which couple to the W boson, including CKM non-diagonal combinations. The simulation has been carried out for the design LHC condition, proton-proton collisions at the CM energy of 14 TeV, using a leading-order parton distribution function (PDF) CTEQ6L1 [26]. The renormalization scale (μ R ) and the factorization scale (μ F ) are set to the W -boson mass (m W = 80.419 GeV). The cut, p T ≥ 5 GeV and R j j ≥ 0.2, is applied to the two final-state partons in order to cut off the divergence, where R j j is the separation between the two partons defined as the quadratic sum of the separations in the pseudorapidity and azimuthal angle, In order to deal with 2-jet production MEs, we have to subtract the final-state divergence that appears at the limit where the two final-state partons are produced collinearly, together with the initial-state collinear divergence that has been subtracted from 1-jet production MEs. We import the subtraction technique that we have developed and tested for QED photon radiation in the final state [17]. The definition of the approximation in Eq. (1) is also valid for the final-state divergence; the difference lies only in the definition of the splitting function.   Fig. 1. p T distribution of the smaller p T parton in the qg → W q g reaction after the LLL subtraction. The simulation is carried out for the 14 TeV LHC condition. The cut, p T ≥ 5 GeV and R j j ≥ 0.2, is applied to cut off the remaining divergence. The distribution of positive-weight (+1) events is illustrated with the dashed histogram and that of negative-weight (−1) events with the dotted histogram. The unweighted result obtained from the difference between the two distributions is shown with the solid histogram.
We evaluate the approximation in Eq. (1) if the parton that is taken as the radiation is judged to be soft. The radiation is defined as soft when the Q value for the assumed radiation-emitter combination is smaller than a given energy scale. We define this energy scale to be equal to m W for both ISR and FSR. The subtraction component is determined by adding the approximations for all possible assignments of the radiation.
In general, the non-radiative event that is defined by removing the soft radiation may again include a soft radiation when we treat multi-jet events. In order to make such doubly soft contributions finite, we need to subtract next-to-leading-logarithmic (NLL) components together with LL ones. The evaluation of NLL terms is beyond the scope of the present study. Hence, we simply ignore their contributions. We terminate the evaluation and set the ME value to zero once such a doubly soft combination is found.
Since the subtraction is unphysical, the remaining cross section may become negative in some phase-space regions. GR@PPA generates events having weights of +1 and −1 according to the sign of the cross section. The spectra of the positive-and negative-weight events are separately shown with dashed and dotted histograms, respectively, in Fig. 1. The final spectrum shown with the solid histogram is obtained as the difference between them. The obtained spectrum indicates a negative divergence at p T → 0. Namely, the collinear approximation in Eq. (1) overestimates the divergence in the soft region.
The properties of another divergence, the soft-gluon divergence, have been discussed by Catani and Seymour [27]. They have shown that the asymptotic form of this divergence, the soft-gluon approximation, can be expressed in a form similar to the collinear approximation in Eq. (1) as 2 where q denotes the four-momentum of the considered soft gluon. In this approximation, we have to consider the spectator b together with the emitter a because this divergence originates from the interference between the radiations from different external legs. The variables p a and p b represent their four-momenta, and T a · T b is the corresponding color charge. This divergence emerges in gluon radiations, while no such divergence appears in quark radiations. Starting from this expression, Catani and Seymour introduced the dipole subtraction method for next-to-leading-order calculations, in which the proposed subtraction function converges to the collinear approximation at the collinear limit and to the soft-gluon approximation at the soft limit of the gluon radiation. In the following, we introduce another form of the subtraction function that is suitable for application in our event generator.
Here, we have to note that the color factor is assumed to factorize in the expression in Eq. (2). However, the factorization is not trivial. The amplitude of a diagram m for jet production processes can be subdivided into a color part M m C and a kinetic part M m p . If we do not assume the color factorization, the soft-gluon approximation can be written as where M m 0 ⊗R a C represents the color part for the radiative diagram that is transformed to the nonradiative diagram m 0 by the removal of a gluon attached to the external leg a, and the kinetic part for the non-radiative diagram m 0 is written as M m 0 0, p . If the color part factorizes, i.e., the ratio of the color parts, is independent of the diagram combination (m 0 , n 0 ), the ratio can be considered as the color charge, Y m 0 n 0 ab = T a · T b , and the expression in Eq. (3) can be reduced to Eq. (2). In general, the color part factorizes only when the number of colored external legs included in the non-radiative event is two or three [27], whereas it does not factorize for more complicated events. Actually, we confirmed that the ratio in Eq. (4) depends on the diagram combination in W + 3 jet production processes by explicitly calculating the ratio. However, we confine ourselves to 2-jet production processes in this article. The color part factorizes in this case. Hence, in the following, we assume that the soft-gluon divergence can be approximated by the function in Eq. (2).
First of all, we compare the collinear-limit behavior of the soft-gluon approximation with the softlimit behavior of the collinear approximation. We assume hereafter that all partons are massless. Hence, the a = b term vanishes in Eq. (2). As the energy fraction of the radiation can be taken as 1 − z a , the four-momentum of the radiation can be approximated as q → (1 − z a ) p a for ISR and q → (1 − z a ) p a /z a for FSR at the collinear limit. Since Q 2 a = −( p a − q) 2 for ISR and Q 2 a = ( p a + q) 2 for FSR, the denominator, p a q, in Eq. (2) is equal to Q 2 a /2 both in ISR and FSR. Therefore, since the soft-gluon approximation in Eq. (2) can be further approximated as at the collinear limit. In order to examine the soft-limit behavior of the collinear approximation, we define the function P a (z) as The color factors for gluon radiations are T 2 q→qg = C F = 4/3 and T 2 g→gg = C A = 3. The definitions for ISR and FSR are slightly different from each other. The explicit definitions are given asP Note that all theP a (z) functions approach unity at the soft limit, z → 1. Sinceŝ/ŝ 0,a → 1 at the soft limit, we obtain the soft-limit behavior of the collinear approximation in Eq. (1) as Namely, the collinear limit of the soft-gluon approximation is identical to the soft limit of the collinear approximation.
Although the above conclusion may sound trivial, we gain a crucial insight into divergence approximations from this discussion. Let us consider a modified soft-gluon approximation defined as SinceP a (z a ) → 1 andŝ/ŝ 0,a → 1 at z a → 1, this approximation approaches the soft-gluon approximation at the soft limit. Simultaneously, because of the property in Eq. (6), the approximation in Eq. (10) approaches the collinear approximation in Eq. (1) at the collinear limit. Therefore, we must be able to subtract all divergences in multi-jet MEs by using the approximation in Eq. (10) instead of the collinear approximation in Eq. (1). However, direct use of the approximation in Eq. (10) is problematic; the cross-section integration is hard to converge. We apply the subtraction not only in regions close to the soft limit, but also significantly away from the limit. The approximation in Eq. (2) exhibits unexpected behavior in the latter. It is better to eliminate the contribution of the soft-gluon approximation away from the soft limit. As a solution, we introduce a combined approximation that is defined as where r soft is an arbitrary dumping coefficient that should approach unity at the soft limit and reduce to zero far away from the limit. It is obvious that this approximation has the same properties as that in Eq. (10) at the soft limit and the collinear limit. We have tested several definitions of r soft and found that the simplest definition, provides good convergence properties in all cases that we have tested. Hence, we use this definition through the following discussions. It should be noted here that, if there are only two colored particles in the non-radiative event, the soft-gluon approximation in Eq. (2) coincides with the soft limit of the collinear approximation, Eq. (9), without taking the collinear limit if it is evaluated in a frame in which the two colored  The simulation condition and the histogram notation are the same as those of Fig. 1, except that the cut for the produced partons is relaxed to p T ≥ 1 GeV and R j j ≥ 0.01.
particles are aligned back-to-back. Thus, it is not necessary to account for the soft-gluon divergence separately. This condition applies to the 1-jet production processes in hadron collisions that we have studied previously. Other cases where this condition holds include the e + e − → qqg process in the frame in which the qq pair is produced back-to-back, and the deep inelastic scattering, q → qg, in the Breit frame. Figure 2 shows the p T spectrum of the smaller p T parton in the qg → W q g reaction after the combined subtraction, in which the combined approximation defined in Eqs. (11) and (12) is used for the subtraction when the final-state gluon is considered as the radiation, while the ordinary LLL subtraction is applied when q is assumed to be the radiation. The gluon is always treated as the radiation when the FSR component is evaluated. For the calculation of the soft-gluon approximation, the color charge can be assigned as T q · T q = −(C F − C A /2) and T q · T g = T q · T g = −C A /2 from the sum rule in Eq. (5). The total subtraction factor is evaluated by adding all these approximations. The cut is relaxed to the usual one, p T ≥ 1 GeV and R j j ≥ 0.01, that is applied only for numerical stability. The other conditions are the same as those for the simulation shown in Fig. 1. We can see that the negative divergence at p T → 0 has disappeared, although a finite negative offset remains. This offset might be caused by the mismatch in the kinematics model because we are using different models for ISR and FSR. In any case, even if it is due to the mismatch, the existence of such a small offset is harmless since the spectrum at small p T is overwhelmed by PS radiations from non-radiative events.
A similar result has also been obtained for another W + 2 jet production process, qq → W gg. The LLL subtraction provides finite cross sections for the other W + 2 jet processes that include no gluon in the final state: gg → W qq , qq → W q q , qq → W q q , etc. The p T distribution of the smaller p T parton for the combined W + 2 jet production process is shown in Fig. 3. The spectrum is similar to that of the qg → W q g subprocess in Fig. 2 because this subprocess dominates the combined result in the LHC condition.
Since the color factorization is assumed, the above discussions are restricted to 2-jet production processes. However, the factorization is assumed only for proving the identity between the collinear limit of the soft approximation and the soft limit of the collinear approximation. If the identity is   Fig. 3. p T distribution of the smaller p T parton after the subtraction. All W + 2 jet production processes are included. The combined subtraction is applied to the qg → W q g and qq → W gg reactions, while the LLL subtraction is applied to the other processes in which no gluon emerges in the final state. The simulation condition and the histogram notation are the same as those of Fig. 2. generally proved, the above subtraction method can also be applied to larger jet-multiplicity processes. Although it is yet to be formally proved, numerical studies indicate that the identity holds in 3-jet production processes.

Soft-gluon correction to PS-applied events
The subtracted components are restored by PS radiations from lower jet-multiplicity events in our matching method. Since the subtraction is altered by accounting for the soft-gluon divergence, the PS simulation should also be altered accordingly. As we can see in Eq. (2), the soft-gluon approximation does not fully factorize; i.e., the kinetic factor depends on the momenta of particles composing the non-radiative event. Therefore, it is impossible to implement the corresponding correction as an alteration to PS branches. Instead, we adopt a method to alter the production frequency of lower jet-multiplicity events to which the PS simulation is applied.
The divergent components subtracted from W + 2 jet MEs have to be restored by PS radiations from W + 1 jet events. Hence, we discuss the correction to PS-applied W + 1 jet (W + 1 jet ⊗ PS) events in the following. In this PS simulation, the PS has to be applied to the final-state parton as well as the initial-state partons. Since the kinematics model of the PS branches plays a key role in the ME-PS matching [12], the performance of the initial-state (spacelike) PS model has been intensively examined in previous studies on leading-jet matching [22][23][24]. On the other hand, the performance of the final-state (timelike) PS has not been explicitly examined in these studies. The kinematics model of the final-state PS that is implemented in GR@PPA is described in a previous article [15]. The performance of this model has been studied, and necessary corrections to achieve precise matching have been established in the study of diphoton production [17] by using the QED photon radiation from quarks as the probe. We adopt this model and the established corrections in the present study. A similar timelike PS is also applied to the partons radiated in the initial-state PS [15].
The soft-gluon correction to W + 1 jet ⊗ PS events is evaluated as follows. First of all, we determine the hardest PS branch in each event. The search is performed over both initial-state and final-state branches. The hardness is defined by the p T value of each branch, where p T is defined 8 as p 2 T = (1 − z)Q 2 for ISR and p 2 T = z(1 − z)Q 2 for FSR. Subsequently, we boost the included hard-interaction event to its CM frame, and rotate the hardest PS branch so that it aligns along the emitter parton in the hard-interaction event in order to remove the effects of softer branches. The rotation is performed so as to minimize the angular change in order to preserve the angular information of the radiation as strictly as possible. From the boosted hard-interaction event and the hardest branch, we construct a radiative event following the kinematics model that we adopt in PS branches. We evaluate the collinear and combined approximations for this radiative event, and take the ratio, T comb /T coll , as the event weight. The squared ME value is multiplied by this event weight before it is passed to the MC integration and event generation utility BASES/SPRING [28,29] included in GR@PPA using the LabCut framework [15,17]. The PS simulation is allowed before passing the ME value to BASES/SPRING in this framework. In this way, the evaluated event weight is automatically accounted for in the cross-section integration and event generation.
There are some remarks on the hardest-branch search. Starting from the hard interaction, the search is performed in decreasing order of Q 2 for both ISR and FSR. In ISR, the search is done for radiations from spacelike partons only. In FSR, when PS branches attached to a quark in the hard interaction are investigated, only those radiations from this quark line are examined, while all g → gg branches are examined when PS branches attached to a gluon are investigated. The search is terminated if we encounter a flavor-singlet branch, g → qq, or reach the end of the branch chain. We ignore the other branches in order to ignore irreducible higher-order effects included in the PS simulation. For example, when PS branches are attached to a qq → W g event, one of the gluon radiations from the quark that is produced by a g → q q branch of the final-state gluon may become the hardest. Since the property of the gluon radiation of a quark is different from that of a gluon, the observed hardest gluon radiation cannot be considered as a radiation from the qq → W g event; instead, it should be attributed to a radiation from a qq → W q q event. We ignore such higher-order effects by terminating the search at the g → q q branch. If the flavor-singlet branch is the hardest, we return the unit event weight because such an event corresponds to a quark-radiation event for which we do not need to consider the soft divergence.
In Fig. 4, the distribution of the event weight evaluated for qg → W q ⊗ PS events, in which a gluon is selected as the hardest radiation in PS, is compared with the distribution of the T comb /T coll ratio of qg → W q g events generated according to the LLL approximation. The event generation has been carried out in the 14 TeV LHC condition. Although the event weight of qg → W q ⊗ PS events is evaluated, the ME values are not corrected for it in the event generation. The compared distribution of qg → W q g events shows the alteration in the subtraction. Therefore, the comparison shows us how precisely the alteration is compensated by the correction to qg → W q ⊗ PS events.
The cuts are carefully arranged to ensure that the two distributions are comparable. The Q 2 value is required to be larger than m 2 W for q with respect to the initial-state gluon in both qg → W q events in qg → W q ⊗ PS and non-radiative qg → W q events reconstructed in the LLL approximation to qg → W q g. A p T condition, p T ≥ 5 GeV, is required for the hardest branch in qg → W q ⊗ PS. Accordingly, the same condition is required for the gluon in the qg → W q g events with respect to the initial-state partons and to the combined momentum of the two final-state partons.
The two distributions compared in Fig. 4 are in good agreement in their shapes, although there is a significant difference in the normalization. The difference in the normalization is due to the difference in the gluon energy spectrum. The spectrum of the hardest PS branch is suppressed at low energy as an effect of multiple radiation. The comparison is also made in Fig. 4 in three regions of the gluon energy in the CM frame, E cm (g). We can see that the difference in the normalization is  large at low energy and becomes smaller at high energy, while the agreement in the shape holds in all regions. Hence, we expect that the thus-evaluated correction to the qg → W q ⊗ PS events will reasonably compensate for the alteration in the subtraction from qg → W q g MEs in visible highenergy regions. The possible difference in the normalization, particularly in low-energy regions, can be attributed to a higher-order effect that we do not need to worry about. In any case, the soft-gluon correction is at most 10% on average even for very soft gluon radiations. Therefore, even though the soft-gluon correction plays a significant role in the subtraction of divergent components from W + 2 jet MEs, its impact on the finite W + 1 jet cross section is limited.

Matched event generation
The performance of the matching method described in previous sections is tested for W -boson production in the LHC condition in the following. We combine W + 0, 1, and 2 jet production events generated according to W + 0, 1, and 2 jet MEs, respectively, using the GR@PPA event generator. The W bosons are assumed to decay to an electron and neutrino pair and the PS simulation is fully applied down to Q = 5.0 GeV. The combined subtraction described in Sect. 2 is applied in the generation of W + 2 jet events, while the LLL subtraction and the soft-gluon correction described in 10 Sect. 3 are applied in the generation of W + 1 jet events. The non-divergent contribution of doubly soft events is ignored in the generation of W + 2 jet events as described in Sect. 2. The soft-gluon correction is applied not only to qg → W q ⊗ PS events but also to qq → W g ⊗ PS events in order to compensate for the alteration in the subtraction from qq → W gg MEs. On the other hand, no correction is applied to W + 0 jet events. The other generation conditions are the same as those of previous simulations in this article, except that the collision CM energy is decreased to 7 TeV in order to compare the result with experimental data [30]. As described in Sect. 2, the energy scales, renormalization scale (μ R ), factorization scale (μ F ), and scales for ISR (μ ISR ) and FSR (μ FSR ), are all set to m W as the default. The scales μ ISR and μ FSR define the maximum Q 2 of the PS and subtraction in the initial and final states, respectively. Although this choice may not be appropriate for highp T jet production, fixed energy scales are convenient for testing the matching properties since the relation to visible quantities such as the jet p T is straightforward.
The generated events are processed by PYTHIA 6.425 [31] to add softer interactions and hadronization simulations. The default setting in PYTHIA is unchanged except for the settings of PARP(67) = 1.0 and PARP(71) = 1.0, as in the previous studies [22][23][24]. The event selection and the jet reconstruction are applied to the obtained hadron-level events according to the definition in the ATLAS measurement [30]. The events are accepted if the electron and neutrino from the W decay satisfy the following conditions on the transverse momentum ( p T ), pseudorapidity (η), and transverse mass (m T ): where m T is defined as m 2 T = 2 p e T p ν T − p e T p ν T . The Born-level momentum of the electron before applying the photon-radiation correction is used for this selection. The hadron jets are reconstructed by using FastJet 3.0.3 [32], with the application of the anti-k T algorithm with R = 0.4. All stable particles, except for the electron from the W decay and neutrinos, are used for the reconstruction. The reconstructed jets satisfying the condition p jet T > 30 GeV, |η jet | < 4.4, and R(e, jet) > 0.5 (14) are assumed to be detected, where R(e, jet) denotes the separation from the decay electron in R.
The cross section obtained from the simulation is presented as a function of the threshold jet multiplicity in Fig. 5. The statistical error of the simulation is illustrated using boxes. The result is compared with the corresponding ATLAS measurement at 7 TeV [30] plotted with error bars. It is to be noted that no error bar is drawn for the measurement in the first bin because the measurement value is not available in the data archive [33]. The plotted value, which corresponds to the total cross section, has been derived from the σ (≥ 1 jet) result and the result for the σ (≥ 1 jet)/σ (≥ 0 jet) ratio available in the archive.
We can see that the simulation is in good agreement with the measurement up to the threshold jet multiplicity of two and significantly underestimates the measurement for higher multiplicities. This result is reasonable since we have fully evaluated the jet production up to two jets in the simulation. Events with three or more jets are generated by PS radiation covering a limited Q 2 range. It should be noted that, since the MEs for the event generation are all evaluated at the tree level, there is no guarantee concerning the absolute value of the cross section. The data/simulation ratio in the total cross section that is shown in the first bin is the so-called K-factor. The agreement up to two jets is further improved if we correct the normalization of the simulation with this factor (= 1.11).   The matching in the combination of W + 0 and 1 jet events has already been discussed in previous studies. The soft-gluon correction applied to W + 1 jet events does not have a significant impact on the discussion. The matching of the newly included W + 2 jet process can be studied by investigating the property of the second jet, where the detected jets are ordered according to the p T value. Figure 6 shows the p T spectrum of the second jet. Together with the summed spectrum, the contributions from W + 0 jet (w0 j), 1 jet (w1 j), and 2 jet (w2 j) events are presented separately. We can see that the contributions from w0 j and w1 j, which are predominantly determined by PS, are dominant at low  p T , while w2 j is dominant at high p T . The long highp T tail of w2 j shows a characteristic powerlaw behavior expected from W + 2 jet MEs. The sum of the three contributions shows a smooth transition from the low-to highp T regions, to yield a spectrum which is in good agreement with the plotted ATLAS measurement. Although the smooth jetp T spectrum observed in the above is already satisfactory evidence for the ME-PS matching, a further test can be carried out by investigating the stability of the spectrum against variation in the energy scales. The second-jet p T spectra for the choice of μ F = m W /2 and μ F = 2m W are presented in Fig. 7, together with the default result with μ F = m W . The renormalization scale μ R for ME calculations is fixed to m W in this study. The other scales, μ ISR and μ FSR , are set equal to μ F to test the matching. We can see apparent μ F dependences in the separate w0 j + w1 j and w2 j contributions. These dependences compensate each other in the combined result to yield a smooth spectrum which is significantly less dependent on μ F .   The w2 j contribution has an apparent scale dependence since the inclusion of the divergent LL component varies according to μ F . This dependence may induce a scale dependence in the p T spectra of the leading jet and W . Figure 8 shows the p T spectrum of the leading jet, i.e., the highestp T jet. The contributions from w0 j, w1 j, and w2 j are shown separately together with the summed spectrum. We naively expected that the leading-jet p T would be predominantly determined by the w1 j contribution at high p T . However, as we can see in Fig. 8, the w2 j contribution is comparable to or larger than the w1 j contribution at p T μ F . Hence, the induced scale dependence may become significant in this p T spectrum.
The μ F dependence of the leading-jet p T spectrum is shown in Fig. 9. Although the simulation shows reasonable agreement with the measurement, the μ F dependence is not small at high p T , as we were concerned. A similar amount of μ F dependence is also observed in the W -boson p T

14/16
Downloaded from https://academic.oup.com/ptep/article-abstract/2015/5/053B04/1531060 by guest on 28 July 2018 spectrum since the W -boson p T is predominantly determined by the leading jet at high p T . In order to obtain a stable spectrum, the unavoidable μ F dependence of w2 j has to be compensated with the μ F dependence of the w1 j contribution, which can be merely induced by the μ F dependence of PDF. However, the w1 j contribution is very stable against μ F variation in the relevant p T range. This stability is accidental since the μ F dependence of PDF is determined by the combination of the QCD evolution and the functional form of PDF, i.e., the dependence on the momentum fraction of partons. The μ F dependence observed in Fig. 9 is thus unavoidable in our matching method.
As we can see in Figs. 7 and 9, the simulation tends to yield larger cross sections than the measurement at high p T . This tendency may be attributed to the inappropriate choice of the energy scales, particularly to the small μ R value (= m W ) compared to the jet p T . The capability of the simulation as a measurement tool will be discussed elsewhere by carrying out a detailed comparison with measurements. A more appropriate definition of the energy scales will be pursued in the course of the study.

Conclusion
We have extended the LLL subtraction method for ME-PS matching, which we have developed for combining 0-jet and 1-jet production processes in association with electroweak-boson production in hadron collisions, to 2-jet production processes. We have introduced an approximation that simultaneously reproduces the collinear divergence and the soft-gluon divergence in 2-jet production MEs. This approximation is used for the subtraction instead of the collinear approximation in the LLL subtraction in order to make the 2-jet MEs finite. The alteration in the subtraction can be compensated by applying an appropriate correction to 1-jet events to which the PS simulation is applied to generate additional jets. This matching method can be generalized for processes including three or more jets if the identity between the collinear limit of the soft approximation and the soft limit of the collinear approximation is proved.
The performance of the developed matching method has been tested using W + 0, 1, and 2 jet production processes, and the simulation result has been compared with the ATLAS measurement at LHC. The inclusion of the 2-jet process shows a remarkable impact on the second-jet p T spectrum. The 2-jet contribution exhibiting a long tail to high p T is smoothly connected to lowp T contributions from 0-jet and 1-jet processes, yielding a spectrum in good agreement with the measurement.
The contribution from the 2-jet process has an apparent factorization-scale (μ F ) dependence in our matching method. We have shown that this dependence is compensated by the μ F dependence of the 0-jet and 1-jet contributions in the second-jet p T spectrum. On the other hand, this μ F dependence induces a sizable μ F dependence in the leading-jet p T spectrum at high p T . This is caused by the fact that the μ F dependence of the W + 1 jet contribution at high p T , which is determined by PDF, is accidentally very small.