Resummation prediction on the jet mass spectrum in one-jet inclusive production at the LHC

We study the factorization and resummation prediction on the jet mass spectrum in one-jet inclusive production at the LHC based on soft-collinear effective theory. The soft function with anti-$k_T$ algorithm is calculated at next-to-leading order and its validity is demonstrated by checking the agreement between the expanded leading singular terms with the exact fixed-order result. The large logarithms $\ln^{n} (m_J^2/p_T^2)$ and the global logarithms $\ln^{n} (s_4/p_T^2)$ in the process are resummed to all order at next-to-leading logarithmic and next-to-next-to-leading logarithmic level, respectively. The cross section is enhanced by about 23% from the next-to-leading logarithmic level to next-to-next-to-leading logarithmic level. Comparing our resummation predictions with those from Monte Carlo tool PYTHIA and ATLAS data at the 7 TeV LHC, we find that the peak positions of the jet mass spectra agree with those from PYTHIA at parton level, and the predictions of the jet mass spectra with non-perturbative effects are in coincidence with the ATLAS data. We also show the predictions at the future 13 TeV LHC.


I. INTRODUCTION
The substructure of jets produced at the Large Hadron Collider (LHC) has become one of the hot topics for both theorists and experimentalists. The particles such as massive electroweak bosons, top quark and other possible new resonances produced with transverse momenta much greater than their masses, i.e., p T ≫ m, can decay to hadronic products, which are almost collinear and may be recombined into a single jet by jet algorithms. Therefore it is necessary to find a way to distinguish the interesting signal jets from the purely QCD backgrounds.
The event generators such as SHERPA [14,15], PYTHIA [16,17] and HERWIG++ [18,19], can provide fully differential events, by which any observable can be predicted and compared with data. However, the various event generators employ different models for parton shower and non-perturbative effects, such as the hadronization and multiparton interactions. As a consequence, they might provide very different predictions. For instance, the jet mass spectra from the PYTHIA and the HERWIG++ do not agree with each other, as shown in Ref. [20]. Moreover, there is a type of color correlation between the initial and final colored particles that is not taken into account in these event generators.
In order to obtain more precise predictions and test the validity of the Monte Carlo tools, it is important to develop a theoretical framework to study the jet substructure.
The theoretical developments of prediction on jet mass spectrum at hadron colliders can be found in [28,32,33,35]. In Ref. [32], the jet mass was investigated with the pQCD resummation formalism by focusing on the processes independent jet function, where it was found that the nonperturbative effects are important at small jet mass. The author of Ref. [33] studied the distributions of m J /p J T in pp → dijet and Z + 1 jet processes at NLL, using the formula in Refs. [36,37], and including resummation effects of non-global logrithms (NGLs) in large-N c approximation. The jet mass spectrum with the Higgs + 1 jet process was discussed [35] in the N-jettiness global event shape [38]. The factorization formula and resummation prediction of the jet mass spectrum for direct photon production in the framework of SCET was provided in Ref. [28], where the soft function was factorized into two pieces with different scales. Thought the non-global logarithms were not resummed there, their contribution were estimated and it was found that the NGLs only affect the jet mass spectrum in the peak regions significantly.
Studies of the jet mass can not only help us understand QCD, but also be useful to search for new physics, especially in the complex QCD environment of the LHC. In particular, if we want to identify the mass peak of a highly boosted particle, the jet mass spectrum of QCD background must be calculated precisely. Actually, the jet invariant mass were explored in both ATALS and CMS collaborations at the 7 TeV LHC [20,39]. From these results, we can see that the jet mass m J peaks at about 50 GeV, which can be much smaller than the transverse momenta of jet p T . Therefore there exist large logarithmic terms α n s m 2 J ln m (m 2 J /p 2 T ) with m ≤ 2n − 1 in the perturbative calculations near the peak region, which need to be resummed to all order in order to give reliable predictions.
In this paper, we study one of the simplest jet substructures, i.e. the invariant mass of a jet, and investigate the factorization and resummation prediction on the jet mass spectrum in SCET for one-jet inclusive production at the LHC. Compared with direct photon process [28], the factorization formula for dijet process is more complicated due to the nontrivial color structure and associating soft radiation. The illustrative picture of this process is shown in Fig. 1. Since the soft radiation can either be inside or outside the cone of the measured jet, there are two kinematic variables which can lead to large logarithms at threshold limit: one is the invariant mass m J of the measured jet, and another is the invariant mass √ s 4 of the partonic system that recoils against the observed jet. In the threshold region m 2 J → 0 and s 4 → 0, both of the large logarithms ln n (m 2 J /p 2 T ) and ln n (s 4 /p 2 T ) need to be resummed to all order. In the threshold limits, the cross section can be factorized as where H, S, J, f P are the hard function, soft function, jet function and parton distribution function (PDF), respectively. Both of the hard and soft function are matrices in color space.
The hard function includes the short distance contributions arising from virtual corrections.
The jet function presents the collinear radiation in the jet. The indices "obs." and "rec." denote the observed jet and the recoiled one, respectively. The effects from soft gluon emission are incorporated in the soft function and its phase space is constrained by the jet algorithms. It is noteworthy that the large angle soft gluon arising from the initial state radiation (ISR) and recoiled final state radiation are taken into account in this formalism.
In contrast to the cone algorithm adopted in Ref. [28], we choose anti-k T algorithm [40] to calculate the jet and soft functions, which is boost-invariant and stable against the change of jet boundary [41]. Thus, our prediction can be valid for the jet with both small and large rapidity, and is more useful for phenomenological purposes.
This paper is organized as follows. In Sec. II, we analyze the kinematics of the one-jet inclusive production at hadron colliders and give the definition of the threshold region. In Sec.III, we derive the factorization formula. In Sec. IV and Sec. V, we show the results of hard function and jet function at NLO, respectively. We calculate the soft function at NLO and present its refactorization in Sec. VI. In Sec. VII, we give the final renormalization group (RG) improved cross section analytically. In Sec. VIII, we discuss the numerical results of the jet mass distribution for one-jet inclusive production at the LHC, including the leading singular distribution at threshold limit, scale uncertainties, R dependence, distinction between quark jets and gluon jets, and comparison between the RG improved predictions and ATLAS data. We conclude in Sec. IX.

II. ANALYSIS OF KINEMATICS AND FACTORIZATION
In this section, we introduce the relevant kinematical variables and the factorization formula needed in our analysis. We consider the process where J denotes the leading final jet, and m J is its invariant mass. The partonic channels include qq → qq, gg → qq, gg → gg and their various crossing ones. The Feynman diagrams at leading order (LO) are shown in Appendix A.
It is convenient to introduce two lightlike vectors n µ a = (1, 0, 0, 1) and n µ b = (1, 0, 0, −1) along the beam directions, and another lightlike vector n J = (1,n J ) along the measured jet direction. In the center-of-mass (CM) frame of the initial partons, for the one-jet inclusive production, the momentum of recoiling parton to the observed jet is along the direction n J = (1, −n J ). In the CM frame of the hadronic collision, the momenta of the incoming hadrons are given by Here E CM is the CM energy of the collider and we have neglected the mass of the hadrons.
The momenta of the incoming partons, with a light-cone momentum fraction of the hadronic momenta, are Any four vector can be decomposed along the light-like reference vector n i Hence the momentum p µ can be denoted by p µ = (p + , p − , p ⊥ ). The momentum modes relevant to our discussions are the collinear mode p µ n J ∼ √ŝ (λ 2 , 1, λ), anti-collinear mode p μ n J ∼ √ŝ (1, λ 2 , λ) and soft mode p µ s ∼ √ŝ (λ 2 , λ 2 , λ 2 ), where λ = m J / √ŝ is treated as a small expansion parameter. In the partonic threshold limits m J → 0 and s 4 → 0, the radiation is constrained to be either soft or collinear with the final-state partons.
In order to identify energetic cluster of radiation, the sequential recombination jet algorithms are used. The longitudinal boost invariant distance measures d ij and d iB are defined by where R is the jet radius parameter, y i and φ i are rapidity and azimuthal angle of the jet i, respectively. α = −1, 0 and 1 represent the inclusive anti-k T [40], Cambridge-Aachen [42,43] and k T [44,45] jet algorithms, respectively. The effects of jet algorithms on the resummation have been studied in Refs. [41,[46][47][48][49], among which Ref. [41] has shown that jet boundary can be changed significantly by boundary clustering for Cambridge-Aachen and k T algorithms, while the change of the phase space is power suppressed for anti-k T algorithm. In this paper, the anti-k T algorithm is adopted, and the jet boundary is just a circle of radius R in φ − y plane around the jet direction.
After clustering jets, the jet invariant mass m J receives contribution from the radiation inside the jet, whether from collinear and soft gluons. Thus we split the soft radiation k µ to two parts, denoted by k µ = k µ in + k µ out . Then, the partonic threshold variables take the form For later convenience, we write k in ≡ n J · k in and k out ≡n J · k out .
The hadronic threshold is defined as M 2 X → 0. In this limit ,the final state radiations and beam remnants are highly suppressed, which leads to final states consisting of two narrow jets, as well as the remaining soft radiations. For convenience, we introduce the In terms of m X , x 1 , x 2 and v, In the limit This expression is helpful when we derive the RG equation of the soft function by using the RG invariance in Sec. VI.

III. FACTORIZATION IN SCET
To derive a factorization formula for dijet process in SCET, we first have to match the full QCD onto the effective theory [50,51]. To illustrate the factorization in detail, we consider the process qq ′ → qq ′ . The initial partons are labeled by 1 and 2 and the final partons are labeled by 3 and 4, and the relevant operator in QCD is given by [52] O QCD where c I denotes a 4 order color tensor with color indices a i , and Γ (Γ ′ ) denote the chirality (P L or P R ). In SCET, the n-collinear quark field ψ n can be written as where W † n is the Wilson line, and χ n is the gauge invariant combination of W † n and collinear quark field ξ n in SCET. At the leading power in λ, only the n·A s component of soft gluons can interact with the n-collinear field χ n (x), which can be decoupled by a field redefinition [53]: with Then the effective Lagrangian can be expressed as with Here C Γ I is the hard matching coefficient. The scattering amplitude for the qq ′ → qq ′ can be written as where |C Γ is the vector of Wilson coefficient combination in color basis |c I , as following For qq ′ → qq ′ , the color basis is chosen as The differential cross section can be written as dσ dp T dydm 2 where Since the soft and collinear sectors are decoupled due to field redefinition, the matrix element in Eq. (25) can be factorized into a product of several matrices, where N init = 1/(4N 2 ) denotes the average over the colors and spin of the initial-state partons, and α1, β1, etc, are Dirac indices. The initial state collinear sectors match to the conventional PDFs: and the matrix elements of the collinear fields in the final state match to the quark jet function: The soft function can be defined as the matrix element associated with the soft Wilson line which can be decomposed in the color basis Now the matrix element appearing in Eq. (28) can be simplified as All the above components in the factorization form in Eq. (28) satisfy certain RG equations, which we will discuss in the following sections. Combining the different parts together, we get the factorized differential cross section in the threshold limits dσ dp T dydm 2 where C ij is the hard-scattering kernel And H IJ is the hard function, the details of which are shown in Sec. IV.
For other channels, such as gg → qq ′ or gg → gg, the formula of factorization is similar to the process qq ′ → qq ′ , except for the different jet functions and PDFs. The definitions of gluon PDF and jet function are given by and

IV. HARD FUNCTION
The coefficient C Γ I can be obtained by matching the full theory onto SCET. The one loop results for all partonic 2 → 2 process in QCD have been available in Ref. [52], which are derived in dimensional regularization and the MS renormalization scheme. In this section, we show the crossing relations for different channels and the RG evolution briefly. The explicit expressions of hard matching coefficients are shown in Appendix B.
The Wilson coefficients for the channel qq ′ → qq ′ are denoted by C Γ I (s, t, u) and the others can be obtained by crossing symmetries, as shown in Table I. For example, the Wilson coefficients for the channel qq → q ′q′ are C Γ I (u, s, t). Γ in Wilson coefficients denotes the chirality of the incoming and outgoing partons. In general, there are 16 possible chirality amplitudes. Actually, for the channel qq ′ → qq ′ , only 4 chirality amplitudes are non-zero. This is because that chiralities of massless particles 1 and 3 (2 and 4) must be the same.
We rewrite the Wilson coefficients as C λ 1 ,λ 2 I ≡ C λ 1 ,λ 2 ,λ 3 ,λ 4 I with λ 1,2 = L or R. In addition, since the chirality can be changed by charge conjugation, the only two independent chirality amplitudes for qq ′ → qq ′ are C LL I and C LR I . If the 4 quarks are identical, there are two additional non-vanishing chirality amplitudes C LRRL I and C RLLR I because of the contribution of u-channel. The interference between tchannel and u-channel also makes the results different from qq ′ → qq ′ case. The results for qq → qq can be expressed as where The results of other channel associated with qq → qq can be obtained by crossing symmetry as shown in Table I.
Next, we consider the Wilson coefficients for gg → qq channel and its crossing. There are six relevant channels The Wilson coefficients for the channel gg → qq are denoted by C λ 1 ,λ 2 ,λ 3 ,λ 4 I (s, t, u) and the others can be obtained by crossing symmetries as shown in Table I. According to parity invariance, we have In addition, C λ 1 ,λ 2 ,λ 3 ,λ 4 The explicit expressions of M Γ I and Q Γ I are listed in Appendix B for the convenience of the reader.

B. RG evolution of the Hard Function
The Wilson coefficients C Γ I satisfy the RG equation [54][55][56][57][58] where Γ H IJ can be expressed as with and where β(α s ) is the QCD beta function, γ cusp is the cusp anomalous dimension, and n q and n g is the number of external quarks and gluons involved in the process, respectively. M IJ (s, t, u) denotes the color mixing terms, and can be written as where s 12 = s 34 =ŝ, s 13 = s 24 =t, s 14 = s 23 =û, and L(x) is defined as The explicit expressions of M IJ for each channel can be found in Appendix B. M IJ can be diagonalized with eigenvalues λ K . For example, for qq ′ → qq ′ channel, we have where F (s, t, u) denotes the transform matrix, which can be calculated numerically. The Wilson coefficients in the diagonal basis are denoted byĈ Γ K ≡ F KI C Γ I , which satisfy the RG equation The hard function in the diagonal basis is denoted byĤ the RG equation of the hard function can be obtained, Solving the RG equation, we can get the resummed hard function where S(ν, µ) and A Γ (ν, µ) are defined as Up to NNLL level, we need three loop γ cusp and β function and two loop γ H , and their explicit expressions are collected in the appendix A of Ref. [59].

V. JET FUNCTION
The jet functions J(p 2 , µ), defined in Eqs. (30) and (38), describes a collinear quark or gluon with the invariant mass p 2 . It is process independent and has been calculated at NLO in Ref. [60] and NNLO in Refs. [61,62], respectively. The nonvanishing diagrams contributing to the NLO jet function in Feynman gauge are given in Fig. 2. The relevant diagrams corresponding to quark jet function are shown in the top row, and the ones corresponding to gluon are shown by the others. The RG evolution of the quark jet function is given by The gluon jet function is the same with C F replaced by C A and γ Jq replaced by γ Jg , respectively. This evolution equation can be solved by the Laplace transformation [63,64]: which satisfies the the RG equation Thus the jet function at an arbitrary scale µ is given by where η jq = 2C F A Γ (µ j , µ) and η jq = 2C A A Γ (µ j , µ). The µ-dependent part of the Laplace transformed jet function j(L, µ) is determined by the anomalous dimensions of the jet function as in Eq. (58), while the µ-independent part can be obtained by a fixed-order calcula- When jet algorithm is applied, the phase space for the collinear radiation need to be constrained. For anti-k T algorithm, the restriction of phase space is given by p and k denote the incoming and loop momenta, respectively. In the threshold limit p 2 → 0, the integral region of momentum k can be expressed as To avoid the double counting of the soft region covered by soft function, the zero-bin subtraction [65] should be considered. Taking the soft limit of the restriction in Eq. (61), the phase space of zero-bin region is given by [22,24] After integrating the phase space and taking zero-bin substraction, the jet functions with anti-k T algorithm are given by where J q (p 2 , µ) and J g (p 2 , µ) are traditional jet functions. J anti−k T q and J anti−k T g approach to the traditional ones when jet radius R → ∞, which means that there is no restriction to collinear radiation. In addition, it can be seen that the R-dependent terms in Eq. (64) are We have checked numerically that the corrections from R-dependent terms to jet mass spectra are below 1% for m J < 100 GeV, so we will use the traditional jet functions in the following numerical calculation.

VI. SOFT FUNCTION
The soft function defined in Eq. (31) describes soft interactions between all colored particles. When calculating the soft function, we need to consider jet algorithm, which imposes a restriction on the phase space of the soft radiation. In Ref. [28], the jet was defined as all the radiation in a cone of half-angle R around the jet direction, which is different from the one defined by the standard sequencial recombination jet algorithms at hadron colliders.
In this work, we choose anti-k T algorithm to calculate the boost-invariant soft function. In addition, as discussed in Ref. [28], there are multiple soft scales in soft function, which need to be refactorized. In this section, we first discuss the calculation of the NLO soft function with jet algorithm for all channels, and then show its refactorization. The details of the calculations can be found in Appendix C.

A. NLO calculation
As shown in Eq. (32), the soft function S(k in , k out , µ) can be decomposed in color space and calculated in the eikonal approximation. Eq. (24) has shown the color structures for 4-quark channels. For gg → qq and 4-gluon channels, the color structures are given by and , respectively. At LO, the soft functions is given by At NLO, the soft functions can be expressed as follows where i and j index the massless external partons, while I and J index the color structures.
According to Eq. (32), the color matrix (w ij ) IJ can be written as For I ij (k in , k out , p T , y, R, µ), we need to calculate the non-vanishing real emission diagrams in dimension regularization, as shown in Fig. 3, which is given by In the CM of partons, the measurement function M R (k in , k out , R, q) for anti-k T algorithm where y J and φ J presents the rapidity and azimuthal angle of the measured jet. And y and φ are the rapidity and azimuthal angle of the soft gluon with momentum q µ , respectively.
For convenience, we calculate the soft function in the CM frame of initial partons and take φ J to be zero. The results of function I ij are I 23 and I 24 can be obtained by the relations I 23 (k in , k out , y J , R, µ) =I 13 (k in , k out , −y J , R, µ) , In the calculation, we take the limit that R → 0, but we have kept all the terms up to O(R 2 ) so that our result can can be used for a wider range of R.

B. RG equation of the Soft Function
Now, we discuss the evolution of the soft function. Its RG equation can be derived by using the RG invariance of the cross section. In the hadronic threshold limit, we have [59] To transform the convolution form to a product form, using the Laplace transformation we can obtain where s is the Laplace transformed soft functioñ The RG invariance requires And the RG equation of PDF is with for beam N 1 and N 2 , respectively. Herê in threshold limit m 2 J 1 ,J 2 → 0 and w → 1. The gluon PDF equation is the same with C F → C A and γ fq → γ fg . With Eqs. (53), (58), (78) and (79), the RG equation of the soft functionˆ s K ′ K in the diagonal basis is given by where γ S = γ H + γ f i 1 + γ f i 2 − γ J 1 − γ J 2 , C i,j = C F and C A for quark and gluon, respectively.
The relation between the soft functionsˆ s and s iŝ We have checked that the NLO soft function in Eq. (68) satisfies the RG equation (82), which means our factorization is reasonable.

C. Refactorization of the soft function
As shown in Eq. (68), the soft function depends on two variables k in and k out , which are k in ∼ m 2 J /p T and k out ∼ s 4 /p T , in principle. It means that we should treat the two scales separately to control the convergence of perturbative expansion. However, at two-loop level, a complicated dependence on k in /k out will emerge [27,66], which represents the nonglobal structure of the soft radiation. Although we could not ideally factorize the soft function into separate two pieces which depend only on k in /µ and k out /µ, respectively, we can at least extract part of the soft function which depend only on a single scale [28]. We define an auxiliary soft function S in which only depends on k in In color basis, the NLO S in (k in ) can be expressed as where I in ij can be computed by the similar integration in Eq. (70), except for the measurement function replaced by Besides, it is necessary to introduce the residual soft function to describe the soft radiation excluded by S in (k in ) At one-loop level, we consider only one soft gluon emission, which is either inside or outside the jet. It means that S res describe the soft radiation outside the jet, which only depend on k out at O(α s ), and we rewrite it by notation S out (k out ) where I out ij (k out , y J , R, µ) can be calculated according to Eq. (70), through replacing the measurement function by M out (k out , R, q) Now, the soft function at O(α s ) in diagonal basis readŝ s in (L in , µ) andˆ s out (L out , µ) is the Laplace transformation ofŜ in andŜ out , respectively, and their RG equations are where γ in,out are anomalous dimensions depending on the jet radius R, which are given at one-loop level in Appendix C. Solving the RG equation, we get the resummed soft functions with As shown in Ref. [28], the above procedure, so-called refactorization, is an approximate factorization, because the residual soft function would depend on both k in and k out beyond O(α s ). At two-loop level, ln n (k in /k out ) would emerge due to the non-global structure, which has been widely studied at the e + e − colliders [27,30,[66][67][68][69][70][71], but rarely investigated at hadron colliders with a sequencial recombination jet algorithm [33]. A systematical discussion of them requires the complete two-loop results of the soft function with jet algorithms, and is beyond the scope of this paper.

VII. RG IMPROVED CROSS SECTION
From Eq. (35), using Eqs. (54), (59) and (92), we can obtain with And the resummed cross section (34) can be written as dσ NNLLp dp T dydm 2 where NNLL p denotes the approximate NNLL resummation, which means that the NGLs are ignored in this paper. Here, we have changed the integration variables from x 2 to s 4 , which have relation To give precise predictions, we resum the singular terms ln n (m 2 J /p 2 T ) and ln n (s 4 /p 2 T ) in threshold limits to all orders and include the nonsingular terms up to NLO. And the RG improved differential cross section is given by is the weight function, as defined in Refs. [72,73].

VIII. NUMERICAL RESULTS
In this section, we discuss the numerical results for the jet mass distribution in dijet process at the LHC. Throughout the numerical calculations, we use the MSTW2008 PDF sets [74] and associated strong coupling α s . In order to compare with Monte Carlo tools, we use PYTHIA8 [17] with its default "Tune 4C" input. FASTJET [75] is used to perform jet clustering, and the anti-k T algorithm is chosen unless specified otherwise.

A. Leading singular spectrum of jet mass
To verify the correctness of the factorizatoin formula, we expand the Eq. (35) to the leading singular terms (blue dashed line), and compare with the exact NLO results (red solid line), which are obtained from Ref. [76]. From Fig. 4, we can see that the leading singular terms of the jet mass distribution can reproduce the exact NLO jet mass spectrum in small m J region. As m J increases, the difference between the leading singular terms and the exact NLO results increase. We find that in both cases of R = 0.4 and 1, the expanded results agree with the fixed-order ones. This means that our soft function is applicable for not only small R.
where the relevant one-loop contributions get the extreme values.
However, the extreme points of the one-loop contributions of the observed jet function J 1 (µ j 1 ) and soft function S in (µ s in ) do not exist. It can be seen from their NLO corrections where A, B, C are scale independent coefficients. If we measure the jet mass m J , it should not be integrated so that there is no quadratic logarithm term of µ s in and µ j 1 in the one-loop corrections, which is different from the cases of S out (µ sout ) and J 2 (µ j 2 ). As shown from the red line in Fig. 5(c), the NLO corrections to J 1 always decrease with increasing µ j 1 . For where c R is an R-dependent parameter, p * T = 400GeV and µ * = 1.67m 1.47 J (m J in GeV) [28]. According to the extreme points of the variations of S in , c R is numerically determined as 14000 and 2400 for R = 1 and R = 0.4, respectively.
After all of the natural scales involved in this process have been chosen, we discuss the scale dependence of the resummation results of jet mass spectra in Fig. 6. At NNLL p level, three loop cusp anomalous dimension and two loop normal anomalous dimension are used.
For the R-dependent pieces, the one-loop soft anomalous dimensions are used. At NLL level, two loop cusp anomalous dimension and one loop normal anomalous dimension are used. Fig. 6 shows the scale uncertainties for variation of each scales by a factor of 2 about its default value. It can be seen that the scale uncertainties for µ h , µ j 1 , µ j 2 and µ sout reduce significantly from NLL to NNLL p . But for scale µ s in , the NNLL p bands are broader than the NLL ones at large m J region. The reason may be that in large m J region non-singular terms become important and the resummation results are unreliable. In addition, we can also see that the distribution is enhanced by about 23% from NLL to NNLL p at the peak region.
We confirm numerically that this enhancement mainly comes from the one-loop corrections of the hard function, which are included at NNLL p order, but not at NLL. This means that if we want to obtain accurate theoretical predictions, the high order corrections of the hard function must be included.

C. R dependence
In Fig. 7(a), the blue and red solid lines show the results of NNLL p resummation for R = 1 and 0.4 , respectively. We can see that the jet mass spectra shift to right with increasing R, and peak at about 20 GeV for R = 0.4 and 40 GeV for R = 1, respectively. This is due to the fact that when R increases, more large angle soft radiation can be combined into the jet, so that the invariant mass of jet m J = (p c + k s ) 2 become larger. The results from PYTHIA are shown as dashed histograms. Fig. 7(a) shows that the peak positions and shapes of our resummation results agree with the ones of PYTHIA at parton level.

D. The difference of jet mass spectra between quark and gluon
In order to study the difference between quark jet and gluon jet, we show the jet mass distributions for processes with quark and gluon final state separately. In Fig. 7(b), the blue and red solid lines correspond to qq → qq and qq → gg, respectively. The jet mass spectra for quark and gluon jet peak at about 30 GeV and 55 GeV , respectively, which is E. Phenomenological studies of jet mass spectrum at the LHC In this section, we give the RG improved predictions of jet mass spectra at the LHC, and compare them with the results of PYTHIA and the ATLAS data [20]. Fig. 8 shows the normalized jet mass distributions with R = 1 in four different p T bins. At NNLL p + NLO level, the jet mass spectra peak around 25-40 GeV, and shift to right with increasing jet p T .
The peak positions agree with the ones of PYTHIA at parton level. In addition, we can see from the results of PYTHIA that the additional hadronization and multiparton interaction shift the spectra to right by about 10 GeV and 20 GeV, respectively. This means that if we want to obtain predictions which are comparable to data 1 , the non-perturbative effects must be considered. Ref. [80] has computed the non-perturbative corrections to jet mass and their results have been used for Z+1 jet process in Ref. [33], where a shift m 2 J → m 2 J + 2ΩR p T for jet mass has been used to account for the non-perturbative effects. However, as discussed in Ref. [33], this shift in small jet mass region is not meaningful, so we truncate the spectrum in the left side of the peak. Fig. 8 also shows that the NNLL p + NLO results with a shift of Ω = 3.0 GeV (the black solid lines) are consistent with the ATLAS data [20] in all of four p T bins.
Here the shift accounts for the total effects of hadronization and multiparton interaction, so the value of Ω in this work is larger than the one in Ref. [33], where only hadronization is concerned. Notice that our treatment of the non-perturbative interaction effects here is just an approximation. A precise estimate of these effects require some modification of the resummation scheme and global fitting with certain precise data. The further discussion of the non-perturbative effects is beyond the scope of this work, and left in future study. In Fig. 9, we give our RG improved predictions at the 13 TeV LHC. Comparing with the 1 Here we have included the multiparton interaction to the non-perturbative effects for simplicity thought it is not necessarily true. results at the 7 TeV LHC, the jet mass spectra in the same kinematic region shift to right by about 5 GeV. The reason is that the dominated contributions is from qg → qg and gg → gg channel for 7 TeV and 13 TeV LHC, respectively, and the latter channel gives more gluon final states, the average jet mass of which is larger than the one of quark final states.

IX. CONCLUSION
We have studied the factorization and resummation of jet mass for the one-jet inclusive production at the LHC with SCET. The factorization formula is derived systematically.
The NLO soft function with anti-k T algorithm is calculated and its validity is demonstrated by checking the agreement between the expanded leading singular terms with the fixed order results. The soft function is refactorized into two pieces corresponding two different scales. The RG invariance of the cross section is checked at NLO for all channels, which demonstrates the correctness of the factorization. By ignoring the NGLs, we first carry out the resummation at approximate NNLL level. From the numerical results, we find that the jet mass spectrum is enhanced by about 23% from NLL to NNLL p at the peak region.
The enhancement mainly comes from one-loop correction of the hard function. The jet mass spectra shift to right with increasing jet radius R and transverse momentum p T . In addition, we show that there is a significant difference in jet mass spectra between quark and gluon jets. Finally, the normalized jet mass distributions with R = 1 are given in four different transverse momentum regions. We show that the NNLL p + NLO spectra peak at [25][26][27][28][29][30][31][32][33][34][35][36][37][38][39][40] GeV and shift to right with jet p T increasing. The peak positions agree with the ones of PYTHIA at parton level. Including the non-perturbative effects, our results are consistent with the ATLAS data. We also give the RG improved predictions at the 13 TeV LHC and find that the peak shift to right by about 5 GeV comparing with the results at the 7 TeV LHC. Our results are helpful to precisely study jet mass spectrum at hadron colliders and test the validity of the Monte Carlo tools. The Feynman diagrams for dijet process at LO are shown in Fig. 10. 10: Feynman diagrams contributing to the 2 → 2 process at leading order.

Appendix B: Explicit expressions of hard Wilson coefficients
All the expression of Wilson coefficients can be found in Ref. [52]. We list them below the convenience of the reader.
For qq ′ → qq ′ channels, the expressions of the Wilson coefficients in Eq. (40) are given by with If the 4 quarks are identical, the corresponding Wilson coefficients can be obtained by using Eq. (40). The other crossed channels, the Wilson coefficients can be obtained by using crossing relations shown in Table. I.
For gg → qq channel, the Wilson coefficients are given by  where For the other crossed channels, the Wilson coefficients can be obtained by using crossing relations shown in Table. I.
For 4-gluon channel, the Wilson coefficients can be obtained by Eq. (44). The LO matching coefficients M Γ I can be obtained in Table B. At NLO, we also need Q. They can be expressed in terms of A, B and F , the expressions of which are The color matrix of NLO soft function has been defined in Eq. (69). At tree level, the color matrix iss The NLO color matrix is For gg → qq channel, the color matrix at tree level is The NLO color matrix is For gg → gg channel the color matrix at tree level is with The NLO color matrix is with a = − 1 16

Calculation of I ij
Here, we show the detail of the calculation of the I ij function. First, in order to compute I out ij conveniently, we define an auxiliary function I aux ij (k out ) with the measurement function M aux (k out , R, q), which is the same as M in in Eq. 86 except for the delta function. Then I out ij can be obtained by where I full and then For I in 12 , we can get This integration can be computed analytically by approximation at small R sin φ ≈ φ = r sin ϕ , From Fig. 4, we can see that the approximation is validity at even larger R, i.e. R = 1. The other I in ij and I aux γ