From Underlying Event Sensitive To Insensitive: Factorization and Resummation

In this paper we study the transverse energy spectrum for the Drell-Yan process. The transverse energy is measured within the central region defined by a (pseudo-) rapidity cutoff. Soft-collinear effective theory (SCET) is used to factorize the cross section and resum large logarithms of the rapidity cutoff and ratios of widely separated scales that appear in the fixed order result. We develop a framework which can smoothly interpolate between various regions of the spectrum and eventually match onto the fixed order result. This way a reliable calculation is obtained for the contribution of the initial state radiation to the measurement. By comparing our result for Drell-Yan against Pythia we obtain a simple model that describes the contribution from multiparton interactions (MPI). A model with little or no dependence on the primary process gives results in agreement with the simulation. Based on this observation we propose MPI insensitive measurements. These observables are insensitive to the MPI contributions as implemented in Pythia and we compare against the purely perturbative result obtained with the standard collinear factorization.


Introduction
Modern experimental and theoretical studies of processes in hadron colliders are often limited by our understanding of the underlying event which describes all that is seen by the detectors that does not come directly from the primary hard process. In hadronic collisions understanding the various contributions to the underlying event is crucial not only for testing quantum chromodynamics (QCD) but also in searches for new physics and precision measurements. The bulk of underlying event activity comes from multiparton interactions (interactions between the proton remnants from the hard process) and initial and final state radiation (ISR and FSR). Although the contribution to the underlying event from initial and final state radiation can be calculated in perturbation theory, contributions from multiparton interactions are more challenging to estimate. Currently the most effective way for integrating MPI with the hard process and partonic initial and final state showers are through models implemented in Monte Carlo simulations.
In experimental and Monte Carlo studies a class of observables known as MPI sensitive observables are used to probe the underling event activity in hadronic colliders. Transverse energy, E T , is such an observable and is defined as where p (i) T is the scalar transverse momentum of particle i and η (i) is its pseudo-rapidity. For the CMS and ATLAS experiments at the large hadron collider (LHC) the cutoff parameter, η cut , is typically chosen to be ∼ 2 − 2.5 (see, for example, Refs. [1][2][3][4][5]). Other examples of such observables are the beam thrust [3,[6][7][8] and the transverse thrust [9] but in this paper we will focus on transverse energy.
MPI sensitive observables take large contributions from spectator-spectator interactions and it was shown in Refs. [10][11][12] that these contributions are related to the violation of the traditional factorization due to Glauber gluon exchanges. In this paper we do not attempt to prove a factorization formula but we rather adopt an alternative approach where we include multiparton interactions through a model function convolved with the perturbative calculation from the collinear and soft factorization. We study the dependence of the model on the hard scale of the process using Pythia simulations and we find that (for the LHC) below the TeV scale the MPI distribution is independent of the hard scale. The same result was found in Ref. [13] using Herwig++ by studying different primary processes (Higgs, Z, and W ± production). The effect of MPI in Higgs transverse energy distributions was also studied in Ref. [14].
In Ref. [15] it was shown that the factorization of the cross section depends on the region of phase-space under study, even for relatively large rapidity cutoff. Particularly, two regions of phase-space are identified, where r = exp(−η cut ) is the cutoff "radius" and Q the partonic center-of-mass energy. In this work we review the analysis of Ref. [15] and we illustrate how within the framework of soft-collinear effective theory [16][17][18][19] (SCET) we can study the effects of rapidity cutoff on resummed transverse energy distributions measured in hadronic collisions. We use the factorization of Ref. [15] and demonstrate that in the limit E T Qr and with the appropriate choice of dynamical scales, this factorization reduces to the one introduced in Refs. [13,20] for global measurements of transverse energy. In this limit the cross section is independent of the rapidity cutoff up to power corrections of O(Qr/E T ). To simplify the discussion we focus on the Drell-Yan process pp → γ * (→ + − ) + X, where the measurement of transverse energy is imposed on X.
In region I since the cross section is independent of the rapidity cutoff, logarithmic enhancements from non-global effects are not important. However, in region II such effects are expected to become important for E T Qr, we find that for the values of r we are interested resummation of global logarithms alone is sufficient to describe transverse energy distribution where it has significant support.
Since our formalism allows us to calculate the transverse energy spectrum for a wide range of the rapidity cutoff parameter, it can be used for understanding the rapidity dependence of the MPI. For example we found using Pythia that the mean transverse energy from MPI increases linearly with η cut for 1.5 < η cut < 3.5.
Relying on the observation that the model function is insensitive to the hard scale in process we propose an observable defined by MPI-sensitive transverse energy but designed to be MPI-insensitive such that we can make predictions for this observable using the standard soft and collinear factorization formula. We refer to this observable as the subtracted transverse energy and is defined as the difference of the mean transverse energy for two different hard scales, Comparing measurements of this observable against our analytic calculations we can determine if the assumptions made in order to build the model are reasonable. We demonstrate that this observable is independent of MPI contributions for phenomenologically relevant regions, when MPI are calculated using Pythia. Measurement of the mean traverse energy as a function of the hard scale was already performed for various processes in Refs. [4,5]. This subtraction method can be generalized to other additive quantities and as an example, we show that for beam thrust [7,9,20,21] also is insensitive to MPI effects as generated by Pythia.
The factorization of the cross section for regions I and II (see Eq.(1.2)) within SCET is discussed in Sections 2.1 and 2.2, respectively and in Section 2.3 we discuss the merging of the corresponding factorizations with the use of profile scales. We describe the matching onto the fixed order result in QCD in Section 2.4. Furthermore in Section 2.4 we give the assumptions made on the contribution of MPI which leads to the convolution of the perturbative result and a model function for the form of the true cross section. Including the MPI contribution using Pythia we construct a model function for the MPI that gives an accurate description of the the simulation data. The model we construct is independent of the partonic invariant mass. Based on the assumptions that lead to the convolutional form of the cross section we introduce an observable insensitive to MPI in Section 3.1. We confirm that these observables are MPI independent by comparing our purely perturbative results to Pythia simulations. We conclude in Section 4.

Factorization
In this section we illustrate how within the framework of SCET we can reliably describe the transverse energy distribution for the process qq → γ * + X for phenomenologically interesting values of the transverse energy and the rapidity cutoff. It was shown in Ref. [15] that when a rapidity cutoff is imposed the transverse energy distribution is insensitive to the cutoff parameter only in the region E T Qr. In this region the transverse energy spectrum can be described with the factorization theorem for the global case and for this reason we begin in section 2.1 presenting a factorization theorem for the global definition of E T . In section 2.2 we review the factorization of the cross section for E T Qr and in section 2.3 we show how both regions can be described in a single factorization theorem with appropriate choice of dynamical scales which we refer to as profile scales. Finally we discuss the matching onto the fixed order QCD result which describes the region E T ∼ Q in section 2.4

Transverse energy as a global observable: region I
Here the transverse energy is defined as a global observable by, where i extends over all the particles in the event other than the di-lepton pair from the decay of the virtual photon. In the region where E T is parametrically smaller than the invariant mass of the di-lepton pair, the relevant modes to the measurement are the soft and collinear modes and their corresponding scaling is, where p ± and p ⊥ are the light-cone and perpendicular components of momenta with respect to the beam axis. The effective field theory that describes the dynamics and interactions of these modes is SCET II . The hard scaling modes, p µ h ∼ (Q, Q, Q), have been integrated out during the construction of the effective theory. The cross section can then be factorized into hard, soft, and collinear functions [13,20]: where Q and y are the invariant mass and rapidity of the virtual photon and the parton momentum fraction is given by x 1,2 = Qe ±y / √ s. The hard process, qq → γ * (→ + − ) + X is described through the hard function H, which is the the product of matching coefficients from matching QCD onto SCET. The initial state radiation (ISR) from soft and collinear emissions is incorporated within the soft, S s , and beam functions, B q/P , respectively. In addition, the beam functions contain information regarding the extraction of a parton, a, from the proton. The operator definition of the beam function is [6,22], where P n (k) is the proton with momentum k µ = (0 + , k − , 0 ⊥ ), and x B = p − /k − is the fraction of the proton momentum carried by the quark field. Although the beam function is a non-perturbative object for E T Λ QCD it can be matched onto the (also non-perturbative but well known) collinear parton distribution functions (PDFs). This is achieved through a convolution of perturbative calculable matching coefficients and the PDFs evaluated at a common scale, µ, [6]: where I G j/i are the matching coefficients and we use the superscript G to denote that these are the matching coefficients for the case of global measurement in contrast to the case where rapidity cutoff is implemented. We analyze the latter case in the following section. The nextto-leading order (NLO) matching coefficients are given in Appendix A. The soft function can be calculated order by order in perturbation theory using the operator definition and the NLO result is given in Ref. [20]. The Born cross section, σ 0 is defined by All elements of Eq.(2.3) depend on the factorization scale, µ, and thus need to be evaluated at a common scale before combining them to construct a scale independent cross section. For this reason we use renormalization group (RG) methods that allow us to evolve each function from its canonical scale up to an arbitrary scale. This will result in a transverse energy distribution with resummed logarithms of ratios of E T and Q, up to a particular logarithmic accuracy. In this work we will study the next-to-leading logarithmic prime (NLL') accuracy.
The perturbative expansions of the beam matching coefficients and the soft function suffer from rapidity divergences that are not regulated with pure dimensional regularization. For this reason we use the rapidity regulator introduced in Refs. [23,24]. Although the rapidity regulator dependence cancels at the level of cross section, the rapidity scale, ν, introduced during the regularization procedure allows us to resum the complete set of logarithms of E T /Q. It is only after solving the rapidity-renormalization-group (RRG) equations that we may resum all logarithms of E T /Q up to a particular accuracy. Thus the final result for the resummed distribution is where U H and V ss are defined in Appendix C as the solutions to the following RG and RRG equations, More details regarding the RG and RRG properties of the transverse energy or broadening dependent functions can be found in Refs. [20,23]. The canonical scales µ H , µ ss , and µ B are used as the initial conditions for the solutions of the differential equations in Eqs. (2.8) and are chosen such that they minimize the logarithms in the perturbative expansion of the corresponding functions: Similarly for the rapidity scales we have, As mentioned earlier in Ref. [15] it was shown that the the cross section in region I is well described by the global factorization. That means, We use the above equation to describe the spectrum in region I and later in section 2.3 to show that we can describe both regions I and II with a single factorization theorem.

Transverse energy with rapidity cutoff: region II
The transverse energy with rapidity cutoff is defined by, where η (i) is pseudo-rapidity of the i-th particle. As in the global case, we sum over all the particles in the event excluding the di-lepton pair and η cut is the cutoff parameter 1 . Region II, E T Qr, is discussed in detail in Ref. [15]. Here we summarize only the main results necessary for the analysis relevant to this work. As was illustrated in Ref. [15], in this region we can identify an additional soft scale which is collinear enough to resolve the boundary of the rapidity cutoff. This mode was first introduced in Ref. [25] 2 in the context of jet-radius resummation. Adopting the naming scheme of Ref. [25] we refer to this mode as soft-collinear. Thus all the relevant modes are: (u-)soft, collinear, and soft-collinear. The corresponding scaling is, (2.14) These collinear and soft-collinear modes are associated with the direction of one of the beams, similar modes exist for the direction of the other beam. The effective theory that describes these modes is SCET ++ and in this region the cross section factorizes in the following way, where S s (E T ) is the same global soft function that appears in Eq.(2.3) and S n (E T , r) is the soft-collinear function describing the contribution from soft-collinear modes near the cutoff boundary. In the first line we combined soft and soft-collinear functions into the total soft function, S, We note that the soft-collinear functions depend on the rapidity scale ν. This is due to the fact that, compared to the global case, the rapidity divergences (and thus the rapidity scale dependence) appear in the soft-collinear function rather than in the beam function.
The beam functions are rapidity-finite and take, contributions from radiation within two distinct regions of phase-space, below the rapidity cutoff (η cut < η) and beyond the cutoff (η cut > η). Radiation below the cutoff contributes only to the so-called unmeasured beam function which is proportional to δ(E T ) and contributes only to the zeroth bin of transverse energy. Radiation beyond the cutoff will contribute to the beam function through a power corrections of O(E T /Qr). These power corrections could be ignored in the small transverse energy limit but are important in the regime where E T ∼ Qr. Thus the beam function can be written as, B II a/P (x, E T , r; µ) = B q/P (x, r; µ)δ(E T ) + ∆B q/P (x, E T , r; µ) . (2.17) The beam function can be matched onto the collinear PDFs when E T Λ QCD , where the matching coefficient I II a/i can be written as, The first term, I j/i , is the term that determines the unmeasured beam function. The perturbative expansion of this term contains UV divergences that need to be regulated and renormalized. This procedure determines the RG anomalous dimension and evolution of the beam function, B II j/P . The second term in Eq. (2.19), ∆B j/i , gives the contribution to the power corrections that appear in the beam function. This term requires zero-bin subtraction and is finite. The implicit dependence on the factorization scale µ in ∆B j/i is due to the strong coupling constant. The operator definition of the beam function for region II and the one-loop result for the corresponding matching coefficients are given in Section 3 of Ref. [15].
The resummed distribution involves evolving each term in the factorization theorem from its canonical scale to a common scale both in virtuality and rapidity. Similarly to the global case this is achieved through the solution of the corresponding RG equations. For the final result we have where the virtuality scales are We note that, in contrast to the global measurement, the two beams are evaluated at two distinct scales. For central events the two scales are of the same order of magnitude but have different values depending on the rapidity of the virtual photon. This is a consequence of the rapidity cutoff since imposing such a constraint breaks boost invariance. This can be avoided by choosing a dynamic value of the cutoff parameter in a boost invariant way, i.e. η cut (y) = η cut ±y. 3 Our one loop results are then modified with the replacement η cut → η cut (y) and this gives us a boost invariant scale µ II B/B = Qe ηcut . With this choice we ensure that the jet scale is always parametrically smaller than the hard scale for all values of the virtual photon's rapidity. Although a boost invariant definition of the rapidity cutoff is phenomenologically preferred, experimentally fixed cutoff is used and therefore here we proceed with the same choice. The rapidity scales are, In the next section we discuss how modifying these scales and using the factorized cross section in Eq.(2.20) lead to a result that can describe both region I and II with a smooth interpolation in the intermediate regime.

Profile scales and merging
The goal of this section is to show that in the limit r → 0 and E T Qr the factorization for Region II (i.e. Eg. (2.20)) matches onto that for the global measurment (i.e., Eq.(2.7)) with the appropriate choice of dynamical scales which we refer to as profile scales. That is, The exact form of the the profile scales is not important but they need to satisfy the following asymptotic behavior, To see why this set of scales is appropriate for the matching of the two regimes consider evolution kernels that appear in Eqs.(2.7) and (2.20) which are, respectively. In region II where E T Qr, the transverse energy distribution should be described by dσ (II) and thus the introduction of profile scales has no influence in the form of the factorization theorem. This is true since in that region the profiles reduce to the scales that they replace (see first line of Eq.(2.24)). Therefore we have, In the other region, E T Qr, the beam profiles equal the global soft scale, µ ss = E T , thus the beam evolution kernels, U B , reduce to the identity, Since the soft-collinear rapidity profile scale, ν pf sc , is asymptotically reaching the beam rapidity scale for the global measurement, we have, (2.29) For the rest of this section we demonstrate that up to power corrections Eq.(2.29) can be extended to the NLL, NLL', and NNLL cross section. Since this has been shown for the evolution kernels we only need to show the same holds for the fixed order terms at O(α 0 s ) for NLL and O(α 1 s ) for the NLL' and NNLL cross-sections. At O(α 0 s ) this is trivial since both cases reduce to the Born cross section. At O(α 1 s ) we note that the hard function, H(Q; µ), and the global-soft function, S s (E T ; µ, ν), appear in both factorization theorems, therefore is sufficient to show S n (E T , r; µ, ν) ⊗ B II a/P (x, E T , r; µ) Since this task is more technical we leave the details for Appendix B and here we give a phase-space based argument using the corresponding operator definitions. For example the operator definition of the soft-collinear function is [25], which is proportional to δ(E T −E Xsc T (η cut )) with E Xsc T (η cut ) evaluated using Eq.(2.13). Taking the limit r → 0 (or equivalently η cut → ∞), where on the r.h.s. E Xsc T is defined globally. The above equation can be understood in the following way: contributions to the regions of phase-space where particles are emitted within the cone are proportional to the size of available phase-space volume to the power of the number of particles in the cone region (i.e., (V cone ) #-of particles in cone ). Thus, in the small cone limit these corners of phase-space will be suppressed compared to the regions where all particles in |X sc are emitted within the measured region and may be ignored. This corresponds to a global definition of transverse energy. Working in the MS scheme any higher order correction gives scaleless integrals and thus This corresponds to no contribution to the measurement from soft-collinear modes, and suggests that in that limit the soft-collinear modes are redundant. This should be expected since if we define z sc ≡ p − sc /p − c ∼ E Xsc T /(Qr) and demand z sc 1, then as we take the limit r → 0 we unavoidably have E Xsc

Matching onto fixed order
In order to describes the transverse energy spectrum in the region E T ∼ Q we need to match the resummed distribution to the fixed order (FO) result from the full theory. This is necessary in order to include power corrections of E T /Q not described by the effective theory. Furthermore, in this region logarithms of E T /Q are not large and thus the FO result correctly describes the transverse energy spectrum. A smooth interpolation for the intermediate regime can be achieved by adding to the resummed distribution the difference of the full theory FO and effective theory FO result, where we obtain dσ/dE T by integrating over dy in the region y ∈ (−y max , y max ) and over dQ in the region Q ∈ (Q min , Q max ), where Q max , Q min , and y max , define the kinematic cuts for the photon's rapidity and invariant mass.. In order to remain within the central region we need to impose y max < η cut . In Eq.(2.35) dσ FO /dE T is the O(α s ) full QCD result where no rapidity cutoff is imposed. Alternatively one can use, where now dσ FO (η cut ) is the full QCD result where the rapidity cutoff is imposed. The difference between Eq.(2.35) and Eq,(2.36) are power corrections which we already neglected during the construction of the factorization theorem. Note that dσ G, FO /dE T does not depend on the rapidity scale ν since rapidity divergences cancel at fixed order in the convolution of the soft function and the beam functions. The µ scale dependence of dσ G, FO /dE T and dσ FO /dE T comes from the running of the strong coupling and the scale dependence of PDFs. The choice of µ need to be the same for both so that detailed cancellation of the two is achieved in the region E T Q where the resummed distribution describes the spectrum. Detailed cancellation also needs to be achieved between dσ G, FO /dE T and dσ II /dE T in the region E T Q where the fixed order result describes the spectrum. For this reason we need to turn-off evolution at E T Q. This can be easily done choosing µ = Q and using the profile scales in Eq.(2.24) and replacing µ ss , µ B → µ pf ss (E T ) and Comparing our matched NLL' result against a simulation using MadGraph [27]+Pythia [28,29], we find good agreement. The hard process pp → γ * is performed in MadGraph and then showered by Pythia. We use Pythia build-in matrix element (ME) corrections for describing the distribution in the far tail. As discussed in Refs. [28,29] this corresponds up to one additional hard emission from the initial state partons. This is sufficient for our case since we are matching only to NLO corrections in that region. 4 In Figure 1, we show the comparison in the peak region (left) and tail region (right) for the choice η cut = 2.5.
The error band is estimated by varying all scales by a factor of two and one-half around their canonical values. The total error is calculated by adding in quadrature all variations. Caution is necessary here since the scales choice is implemented through the profile functions 4 As discussed in the introduction, we will in Section 3.1 consider moments of this distribution, particularly the first moment. For phenomenological applications that require studies of higher moments, where contributions from the far tail are further enhanced, one needs to consider higher hard parton multiplicity. In Monte Carlo simulations this is achieved through a "Match and Merge" procedure [28,29]. Since here the simulation analysis is included only for purposes of comparison with the analytic NNL'+NLO result the default build in ME correction of Pythia are sufficient.  in order to transition from one region to the other. That requires that a single profile function will change with any of the scale variations in order to ensure the proper transition without double counting the variations. For example, the global-soft and beam profiles will change accordingly when we consider hard scale variation in order to freeze evolution in the far tail but should remain unchanged for E T Q. We collected all the details on the choice of profile functions and scale variation in Appendix D.
In order to compare our analytic result with the partonic distributions in Pythia we turned off the multi-parton interactions and hadronization. The non-perturbative/hadronization effects on the resummed distributions can be studied using the operator definition of the soft and collinear functions [30][31][32][33]. Usually the hadronization effects are included through a convolution of the soft function or the cross section with a model function (which needs to be determined from experiment). The convolution is over the measured observable and thus the model function depends on the observable. The form of the model function is usually determined using the operator product expansion to get the first few moments. This was done for various event and jet-shape observables such as thrust, event-shape angularities, jet mass, groomed-jet mass, D 2 , e.t.c. [34][35][36][37][38]. Ref. [39], studies the non-perturbative effects in transverse momentum dependent (TMD) distributions and jet broadening in e + e − which are most closely related to the measurement presented in this paper.
In contrast, contributions from multi-parton interactions (MPI) are not very well understood theoretically, and a systematic approach for describing these effects has yet to be developed. The subject of MPI and how our formalism can be used to study its effect is discussed in the next section.

Multiparton interactions
The origin of MPI is from secondary interactions of the beam remnants through Glauber exchanges. These interactions are known to break factorization in measurements of global observables but cancel in inclusive cross-sections. A variety of MPI sensitive observables are used in experimental studies for understanding the properties of underlying event (UE) but a comparison to the theory is currently impossible. In this paper we propose a prescription to describe MPI contributions to transverse energy with a rapidity cutoff. In experimental measurements of UE common choices for the rapidity cutoff parameter are η cut = 2 and η cut = 2.5 (for example, see Refs. [1][2][3][4][5]). Our prescription is based on the following two conjectures: • Contributions to underlying event from MPI can be modeled by a convolution of a model function with perturbative results, • The MPI model function is insensitive to hard scale Q.
These assumptions lead to the following expression for the transverse energy spectrum including MPI, where f MPI (E T , η cut ) is the model function that needs to be fitted to the experiment. Similar approach was used in Refs. [38,40,41] in order to incorporate for contribution from UE to jet substructure observables. We allow the model function to depend on η cut to properly incorporate the change in phase-space for different experiments. The dependence on η cut can give us useful information regarding the pseudo-rapidity distribution of MPI in hadronic collisions. Note that the second conjecture can be relaxed allowing the model function to vary slowly with the hard scale. Then instead of Eq.(3.1), the transverse energy spectrum is given by This approach might be more appropriate in studies over an extended range of Q. In this work we consider Q ∈ (100, 1000) GeV and in this region it is sufficient to use the model of Eq.(3.1). For the parameterizations of the MPI model function we used the half normal distribution, where N = 2/(α(η cut ) π) fixes the normalization of the model function to unity and α(η cut ) controls the first moment of the model function, If the conjectures above can be shown to be true up to power corrections, within the effective theory, then α(η cut ) can be written in terms of universal non-perturbative functions such as multiparton distribution functions. Here α(η cut ) can be fixed directly from experimental measurements using, Since no experimental data are available for this measurement (see Refs. [1][2][3][4][5] for relevant experimental studies) we use Monte Carlo simulation data. We find that in the region 1.5 < η cut < 3.5, α(η cut ) can be well described by a linear fit, where A is a parameter that describes the mean transverse energy deposited in the central region from MPI, and depends on the hadronic invariant mass, √ s. Since in this work we are considering only √ s = 13 TeV we treat A as a constant. Fitting to the simulation data, we find A = 22.7 GeV.
In Figure 2, we illustrate the effect of MPI interactions to measurements of transverse energy within a pseudo-rapidity region, as described by Pythia. Once MPI effects are included, the transverse energy distribution differs significantly from the perturbation calculation. On the other hand, by includingthe contribution in Eq.(3.1) with the model in Eq.(3.3) we were able to accurately describe the simulation data.
We emphasize here that the aim of this section is to illustrate that a relatively simple model can describe the contribution of MPI for a large spectrum of the partonic invariant mass. More flexible models can achieve even better agreement, for example one can deviate from the linear fit in Eq.(3.6) allowing A to depend on η cut . Also we could deviate from he functional from of f MPI of Eq.(3.3) (see also the work in Ref. [13] where the dependence of f MPI on s for fixed η cut = 4.5 is discussed).

MPI-insensitive observables
In this section we show how we can use measurements of transverse energy to construct observables independent of MPI. Our proposal depends on the conjunctures above Eq.(3.1) thus the observables we propose can be used either in order to validate these conjunctures or for phenomenological studies, e.g., one can test the conjectures in Drell-Yan and use them in phenomenological studies of Higgs production.
We define the subtracted moments as follows where, E n T (Q) is the n th moment at the hard scale Q and is defined by the following where σ(E T , Q) refers to the differential cross section in E T and Q. Assuming the MPI contribution can be modeled by a function f MPI (E T ) convoluted with the perturbative cross section in Eq. (2.15) we have, The numerator in Eq. (3.8) can be written as where in the second line we changed the integration variable to ω = E T − E T and n C k are the binomial coefficients. With this, Eq. (3.8) can be written as where the perturbative average · · · pert and MPI average · · · MPI are defined by Eq.
T,pert is defined in similar way to Eq. (3.7) in terms of the perturbative cross section. By taking the difference at different values of Q, the first term at k = 0, which is E n T MPI , cancels. Thus for the first few values of n we have, Therefore, at n = 1 the MPI contribution precisely cancels and the difference can be predicted by purely perturbative results. We refer to ∆E T ≡ ∆E (1) T as the subtracted transverse energy. The difference of higher moments includes MPI contributions that can be used to determine parameters of MPI models. Note that the results in Eq. (3.13) are obtained from the two conjectures above Eq.(3.1) but are independent of the model function, f MPI .
We demonstrate the cancelation of the MPI contribution in ∆E T using Pythia simulation in Figure 3. We are comparing the observable ∆E T (Q, 110 GeV) evaluated with the default MPI model of Pythia and the purely perturbative result. The uncertainty here is evaluated by using the maximum and minimum values of the error bands from the perturbative results. It is clear that, within the uncertainty, the observable ∆E T is independent of the MPI contributions, as implemented in Pythia, for that range of the hard scale. In Figure 3 shows this for the cases η cut = 2.5 and η cut = 3.5.

Generalization to other observables
The approach of subtracting the mean at different hard scales to obtain MPI-insensitive measurements can be implemented in other observables as well. This is true for additive observables for which the contributions from MPI can be achieved through a convolution. A characteristic example of this is the beam thrust [3,[6][7][8], B, defined as, 5 (3.14) The corresponding subtracted observable is In Figure 4 we demonstrate that ∆B(Q, Q 0 ) is also insensitive to the MPI contributions for a large range of the virtual photon's invariant mass. An obvious advantage of using beam thrust is that the contribution of each particle is weighted by exp(−η) and thus contribution from particles in the forward region, close to the rapidity boundary, is exponentially suppressed. Therefore the measurement is insensitive to the rapidity cutoff 6 . On the other hand this means that we cannot use beam thrust as our observable if we aim to study the pseudo-rapidity dependence of MPI through the model function f MPI .  Figure 4. The observable ∆B(Q, Q 0 ) (subtracted beam thrust) for Q 0 = 100 GeV and η cut = 4.5. The uncertainty bands correspond to statistical uncertainty due to finite sample data.

Conclusion
In this paper we demonstrated an effective field theory approach for calculating the transverse energy spectrum when a rapidity cutoff is imposed. In contrast to Ref. [15], in this work we use dynamical (profile) scales in order to interpolate between our result from the small transverse energy region (region II) to large values of transverse energy (region I). Finally we used a subtraction scheme to match onto the fixed order result of the full theory (QCD) at E T ∼ Q. Although the cross section in the far tail region where E T ∼ Q is highly suppressed compared to regions I and II, it is important to know the spectrum in all ranges of E T since the far tail region gives significant contributions when we calculate moments of the transverse energy.
As an example we choose to study the process pp → γ(→ + − ) + X away from the Z-pole region. Comparing our results with Pythia (ISR only) simulations we find excellent agreement. We then proceed to introduce a prescription to include the effect of multiparton interactions (MPI). The prescription we propose, which is based on the conjectures above Eq.(3.1), is to simply convolve the perturbative spectrum with a model function. The model function should not depend on the hard process but could have small variations with the change of the hard scale in the process. Comparing our result with Pythia (ISR+MPI) we find that for the range 100-1000 GeV of the photon's invariant mass the MPI model function has little or no hard scale dependence.
Assuming independence of the model function on the hard scale of the problem we introduce an observable, which is independent of MPI effects. This observable we are considering is the subtracted first moment of the transverse energy at two different scales. We compare our purely perturbative calculations of this observable against Pythia (ISR+MPI) and we show that the two agree very well.
Although in this paper we considering only transverse energy measurements, in the last section we discuss generalizations of this MPI insensitive observable to other event shapes, such as beam thrust and transverse thrust. An advantage of using transverse energy as a probe to the MPI effects is that we have strong sensitivity to the rapidity cutoff parameter, η cut . Measurements at different values of η cut can give an insight into the pseudo-rapidity dependence of MPI.
The independence of the MPI to the hard process using Monte Carlo simulations was also investigated in Ref. [13] in Higgs, Z, and W ± production processes implemented in Herwig++ [42]. A future application of our work is to evaluate the Higgs or Z/W ± transverse energy spectrum, using the prescription and model function for MPI we propose in this paper, and compare with results from Monte Carlo simulations.
in part by the Director, Office of Science, Office of Nuclear Physics, of the U.S. Department of Energy under grant numbers DE-FG02-05ER41368. DK would like to thank the Nuclear Theory group at LANL for hospitality during portions of this work.

A The beam function matching for region I
We evaluate the bare global beam function O(α s ) terms using Eq. (25) and (26) from Ref. [43], Expanding first in η and then in we have, where P i/j (x) are the QCD splitting kernels [44,45] andγ q = 3/2. The divergent P q/q (x)/ term cancels during the matching with the collinear parton distribution functions (PDFs) and therefore should not be included in the renormalization kernel. Thus in the MS scheme this yields, and the corresponding renormalization function where we omitted the rapidity and virtuality scale the arguments for simplicity of notation. To evaluate the off-diagonal element I G q/g we start with the corresponding partonic beam function given in Eq.(26) of Ref. [43], Expanding first in η and then in we have, Except the term P q/g (x)/ which cancels during the matching there is no other divergent term. Thus the contribution from the off-diagonal element at this order does not contribute to the renormalization function or the corresponding anomalous dimension. Thus the matching coefficient is, The rapidity and virtuality anomalous dimensions, γ ν and γ µ respectively are evaluated using the following, and thus we have

B Merging fixed order
Here we perform an explicit calculation to show that the product S n ⊗ B II a/P reduces to B G a/P in the large transverse energy limit, E T Qr. To this end, we work with the partonic-level functions, B G q/i in Eq. (2.3) and the combination of soft-collinear and beam function S n ⊗ B II therefore we need to show This task becomes much easier if we work with cumulant bare functions defined where the cumulant of the one-loop beam function in Eq. (2.3) is given by integrating Eq. (A.1), (B.5) The one-loop cumulant soft-collinear function is given integrating Eq.(2.12) of Ref. [15] S (1) where Note that the zero-bin contribution on the last line is precisely cancelled against the soft-collinear function in Eq. (B.6). .

(B.9)
In the last step, we use the identify: Note that the second term on last line is simply power corrections in the limit p cut T p − r: (B.12) There is no zero-bin contribution in the gluon channel thus using Eqs.(B.4) and (B.18) of Ref. [15] I where (B.14) By comparing Eq. (B.12) to Eq. (B.13), we find and we have .

(B.16)
As in Eq. (B.9), the second term on RHS is the power correction.

C Evolution and resummation
In this appendix we give the details for the solutions of the renormalization group and rapidity renormalization group equations. This section is divided into two subsections. In Section C.1 we discuss the virtuality renormalization group equations and the solutions of those equations and in Section C.2 the rapidity renormalization group evolution is described. All elements of factorization (hard, soft, soft-collinear, and beam) satisfy renormalization group equations, in contrast, only transverse energy dependent quantities have rapidity RGE.

C.1 Renormalization group evolution
The RGEs we encounter in this work belong to the same category of what was referred to in Ref. [15] as unmeasured evolution equations. In this paper we do not discuss the evolution of measured quantities and therefore such a distinction is redundant. Also for the processes we are considering the hard and soft function have trivial color structure and therefore we do not address the complications that appear when one considers multi-jet processes in hadronic collisions. The RGEs we consider have the following form, where γ F µ is the virtuality anomalous dimension. We refer to the first term in the square brackets as the cusp part since Γ F µ [α s ] is proportional to the cusp anomalous dimension, and the second term, ∆γ F µ [α S ], as the non-cusp part. Both the cusp and the non-cusp terms have an expansion in the strong coupling. For the cusp term we have, and similarly the non-cusp part is given by, The solution to the RGE in Eq.(C.1) is Since in this work we are interested only in the NLL and NLL' result we may keep only the first two terms in the perturbative expansion of the cusp part (i.e., Γ F 0 , Γ 0 cusp , and Γ 1 cusp ) and only the first term form the non-cusp part (γ 0 F ). Performing this expansion we get, where r = α(µ)/α(µ 0 ) and β n are the coefficients of the QCD β-function, Table 1 the expressions for all ingredients necessary to perform the evolution of any function that appears in the factorization theorems we considered in this paper are given in.

C.2 Rapidity renormalization group evolution
In this section we summarize the solution for the rapidity renormalization group equations for the global soft, S s , soft-collinear, S n , and global beam, B G a/P functions. Even though the unmeasured beam function of region II has transverse energy dependence, it does not have rapidity divergences and thus does not acquire rapidity scale dependence. The RRG equation for transverse energy measurements of the function F (E T ) ∈ {S s , S n , B G a/P } takes the following form, The solution of this equations is where ν 0 is the characteristic scale (for each function) from which we start the evolution. This scale is chosen such that rapidity logarithms are minimized. For the global-soft, softcollinear, and beam functions these scales are given in Table 1. The first term in the rapidity anomalous dimension in Eq.(C.11) is proportional to the cusp anomalous dimension and the proportionality constant we denote with ξ F , (Γ F ν = ξ F Γ cusp ). We define the the plusdistribution in Eq. (C.13) through its inverse Laplace transform,

D Profile functions and scale variation
In this section we give the details on the choice of profile functions and how scale variation is performed. There are four profile functions. They assist either for the transition from region I to region II or for turning off evolution in the far tail (in order to match onto the fixed order result). The asymptotic behavior in all regions for all profile functions is collected in Eqs.(2.24) and (2.37). The relations in Eq.(2.24) need to be always satisfied in order to reproduce the correct result when transitioning from region II to region I. That is, for example, when we vary µ ss (in order to explore uncertainty due to scale variation), µ pf B should be varied accordingly.
There are five scale variations that we need to consider: the hard scale, µ H , which extends to all three regions (I, II, and far tail), the global soft scale, µ ss , which extends only to region I and II, the beam scale, µ II B , which only applies in region II, the global soft rapidity scale, ν ss , for both regions I and II, and finally the collinear rapidity scale which corresponds to ν sc in region II, and to ν B in region I. In practice we will control these variations through five corresponding parameters, H , ss , B , ζ ss , and ζ cs (= +1, −1, 0). The four profile functions and the hard scale are then defined as follows. We choose n = 10. In Figure 5 we illustrate how each of the profile functions changes when we consider the variation for one of the five control parameters: H , ss , B , ζ ss , and ζ cs . In each plot we separate the three different regions with vertical lines. One can also explore the sensitivity of the differential cross section to the choice of the function g(µ 1 , µ 2 ). In our formalism this can be performed by setting all variation control parameters to zero and varying the parameter n. We find that for 4 < n < 12 the result falls within the error bands of the scale variation. We find that n = 10 gives the best numerical stability for the central values. We do not explore different parameterizations of the function g(µ 1 , µ 2 ).