Angularity in DIS at next-to-next-to-leading log accuracy

Angularity is a class of event-shape observables that can be measured in deep-inelastic scattering. With its continuous parameter $a$ one can interpolate angularity between thrust and broadening and further access beyond the region. Providing such systematic way to access various observables makes angularity attractive in analysis with event shapes. We give the definition of angularity for DIS and factorize the cross-section by using soft-collinear effective theory. The factorization is valid in a wide range of $a$ below and above thrust region but invalid in broadening limit. It contains an angularity beam function, which is new result and we give the expression at $\mathcal{O}(\as)$. We also perform large log resummation of angularity and make predictions at various values of $a$ at next-to-next-to-leading log accuracy.


Introduction
Deep inelastic scattering (DIS) [1][2][3][4][5] is an important tool to probe internal structure of nucleus and to test our understanding of Quantum Chromodynamics (QCD) [6,7]. Newly proposed future electron-ion collider (EIC) [8,9] and electron-ion collider in China (EicC) [10] would provide another opportunities to test and to improve our understanding of the structure of nuclei as well as dynamics of QCD in collisions.
One of classic ways of studying QCD events is to measure event shape variables that are observables designed to characterize the geometric shape of hadron distribution in the event. Their values can tell us if the distribution is pencil-like, planar, and spherical. Thrust [11] characterising dijet events is one of classic event shapes predicted to very high accuracy [12][13][14][15][16][17][18] in e + e − collisions and analysis with thrust determined the strong coupling constant α s at 1% level precision [16,17]. The shape variables are also studied at hadron collisions to probe jet substructures [19][20][21][22][23][24][25][26][27]. They were also extensively studied in DIS [28]. Improved theoretical predictions [29,30] and new predictions [31][32][33][34] for event shapes that were not measured before should be valuable for future experiments such as EIC and EicC. Angularity [35] compared to other event shapes was defined relatively later and was not measured in DIS process. One of its features is that the angularity is rather a class of observables defined by a continuous parameter a, which puts weights on rapidity factor and changes relative contributions of particles with large rapidity compared to ones with smaller rapidity. The definition of angularity for a hadronic group X can be written as where the sum is over all hadrons i in the group X, p i ⊥ and η i are transverse momentum and rapidity defined with respect to a given axis. Q N is usually chosen to be a hard momentum scale of the system. The angularity parameter a can be any real number in region −∞ < a < 2. It is known to be infrared-unsafe for a ≥ 2 since particles are weighted rapidity too strongly and the angularity becomes sensitive to collinear splitting. [36]. Under change of a, the value of τ a is sensitive for the collinear with large rapidity while it is less sensitive for the soft with relatively smaller rapidity. So, one can control relative contribution between the collinear and soft with the value of a. The angularity reduces to well-known event shapes: thrust when a = 0 and broadening when a = 1. Factorization in soft-collinear effective theory (SCET) [37][38][39][40][41] near and below thrust region a < 1 was known a while ago [36,42] while broadening region is understood [43] more recently.
Angularity axis was originally chosen to be the thrust axis [35] in e + e − annihilation and the axis and its back-to-back make hadrons grouped into two-hemisphere regions, left H L and right H R , which sum up to give the total angularity τ ee a = τ a (H L ) + τ a (H R ). 1 In jet substructure studies, jet axis and its constituents are usually defined by jet algorithms. Variants of angularity with alternative axis choices [44] can be made out of different purposes.
In this paper, we study DIS angularity defined by two axes, beam (B) axis along proton-beam direction and jet (J) axis associated with a leading jet in the final state as in [29]. For the definition of jet axis, we consider 1-jettiness axis that minimizes the DIS 1-jettiness or, jet axis defined by typical jet algorithms such as cone, C/A, k t and anti-k t . Both axes are equivalent at leading power in small τ a expansion and the axis difference is suppressed by powers of τ a . Of course, one may consider alternative axes such as z-axis in the Breit or, CM frame and broadening axis [44], even multiple axes associated with multi-jet angularity. For a simplicity of discussion we focus on our research on two axes with proposed axis. We measure angularity of a global event obtains contributions from two regions 2 τ DIS a = τ a (H B ) + τ a (H J ) , (1.2) where H B and H J are beam and jet hemispheres each of which has respective jet and beam axes. The events are populated in small angularity region associated with soft and collinear radiations. We study this region with leading-power factorization by using SCET in a < 1 region and make the prediction for τ DIS a at next-to-next-leading log (NNLL) accuracy. The angularity was measured in e + e − annihilation by LEP [45]. It was not measured in DIS while its simplified versions thrust and broadening were measured and analyzed at HERA by the ZEUS and H1 collaborations [46][47][48][49][50][51]. Like other event shapes one can use the measurement to determine the strong coupling constant. Having multiple observables at different values of a would be beneficial in studying correlations between the value of coupling and non-perturbative effect systematically as discussed in [52]. For the given tension in the values of the strong coupling constant determined from the thrust [17] and C-parameter [53] deviated several standard deviations away from the world average significantly weighted by Lattice determinations (FLAG2019) [54,55]. New determination with the angularity from independent process such as DIS provides independent path to address this issue.
Our paper is organized as follows. We first give our definition of axes and express the angularity in terms of Lorentz invariant four-vector products in Sec. 2, show the factorization of angularity cross section in terms of hard, beam, jet, and soft functions derived by using SCET in Sec. 3, then resum logarithms of τ a in Laplace space in Sec. 4. Our angularity beam function at O(α s ) is given in Sec. 5 and numerical results for resummed cross section at NNLL accuracy is given in Sec. 6. Finally, we conclude in Sec. 7.

Angularity in DIS
In this section, we first review kinematic variables of DIS, then give the definition of DIS angularity with the beam and jet axes, finally express the angularity in term of four-vector dot products so that it remains unchanged in any frame under arbitrary boost.
In DIS, an electron with momentum k scatters off a proton of momentum P by exchanging a virtual photon with a large momentum transfer q. The space-like photon momentum can be written by a positive definite quantity Q 2 = −q 2 , where Q sets the momentum scale of the scattering Q Λ QCD . Björken scaling variable x = Q 2 /(2P ·q) ranges between 0 ≤ x ≤ 1 and inelasticity y = Q 2 /(xs) ranges between 0 ≤ y ≤ 1 where s = 2k ·P is the virtuality of electron and proton system. Jet productions are dominated by tree level process, where incoming quark of momentum xP from the proton is struck by the photon and propagates into the final quark with momentum xP + q. The jet momentum P jet should be close to the quark momentum but not the same due to soft radiations induces recoil and change transverse momentum of the jet. 3 q B = xP , where the axis momenta are massless q 2 B = q 2 J = 0. A particle momentum of p can be expressed in terms of these axes and transverse momentum orthogonal to these axes 4 We define a generalized rapidity by using q B and q J , which is not back-to-back and its conjugate rapidity η JB = −η BJ has an opposite sign. If the p is closer to beam axis, the value of η BJ is positive and if it is closer to jet axis, it is negative. This defines beam and jet hemispheres: One can define beam angularity τ B a (p) by replacing η with η BJ in Eq. (1.1) and jet angularity τ J a (p) with η JB . We take τ B a when p ∈ H B and τ J a when p ∈ H J : τ B a (p)θ(η BJ ) + τ J a (p)θ(η JB ). This is equivalent to the angularity with absolute rapidity |η BJ | for all particles in entire regions as in e + e − angularity. We can also use min/max operator and we take smaller one like min{τ B a , τ J a } for a < 1 but larger one like max{τ B a , τ J a } for a > 1. Inserting the transverse momentum in Eq. (2.3) and rapidity in Eq. (2.4) into Eq. (1.1) gives the angularity in term of Lorentz scalars. We choose (1.1) that reduces to the definition of 1-jettiness [29] in a → 0 limit. Applying this for global event, DIS angularity is expressed as The Breit frame makes the picture clear since photon and proton are aligned along z-axis: P = Q xn z 2 and q = Q( nz 2 −n z 2 ) where nz = (1, 0, 0, 1) andnz = (1, 0, 0, −1). Then, qB becomes Qn z 2 and qJ is close to xP + q = Q nz 2 and they are close to back-to-back axis. 4 For back-to-back unit vector n = (1,n) andn = (1, −n) with n·n = 2, p = n·pn 2 +n·p n 2 + p ⊥ .

Factorized cross section
The full QCD cross section for inclusive DIS is conventionally organized into leptonic and hadronic tensors. Here, the hadronic tensor is defined with an additional measurement of angularity and the angularity cross-section for inclusive DIS is written as where the leptonic tensor is given by where α em is the fine-structure constant of QED and Q f is the electric charge of a quark with flavor f . The hadronic tensor defined by QCD current J µ (x) =ψγ µ ψ(x) with a delta function that measures τ a as In second line we removed the sum over X and replace the measure τ a (X) by an operator τ a which acts on state X asτ The operatorτ a can be written in terms of momentum flow operators as discussed in [42] and we bravely omit the discussion by referring to [42]. Now we discuss power counting with a small parameter λ and factorization of angularity in the framework of SCET [37][38][39][40][41]. The parameter λ characterizes the scale of collinear and soft momenta constrained by an observable here τ a or, some cutoffs. To describe momentum scale of soft and collinear modes in light-cone coordinates along beam and jet axes, we define light-like unit vectors n B,J where n i = (1,n i ) with i = {B, J} andn i are unit vectors along beam/jet directions defined in Eq. (2.1), respectively. Conjugate vectorsn i with normalization n i ·n i = 2 can be defined asn Then, we obtain ω B =n B ·q B and ω J =n J ·q J . In the beam region a particle of momentum p can be expressed in terms of n B andn B : p = p +nB 2 + p − n B 2 + p ⊥ or, p = (p + , p − , p ⊥ ) = (n B ·p,n B ·p, p ⊥ ) and in the same way for the jet region. Collinear and soft modes p c,s and angularity τ B a (p c,s ) of each mode in the beam region can be expressed as This implies relevant soft mode contributing to τ B a has a scale of We take λ c as a main power counting parameter and simply denote λ by dropping the subscript c from now on. Power counting of soft mode changes with a. In the limit a → 1, virtualities of soft and collinear modes are of the same order p 2 c ∼ p 2 s and so is their transverse momenta which makes recoil effect of soft radiation important in SCET II . For a below and away from the limit (a < 1), virtuality and transverse momentum of soft mode are parametrically suppressed and this is so called ultra-soft mode, which makes SCET I factorization insensitive to the recoil effect. Here, we pay our attention to SCET I region.
With this observation one can express the collinear momentum as p c =p + k, where the label momentump =n B ·p n B 2 +p ⊥ or,p = (0 ,n B ·p ,p ⊥ ) and a residual momentum k ∼ λ 2 Q is the order of ultra-soft mode λ 2 Q. One can repeat the same power counting for the momenta in the jet region. Now momenta and fields can be expended as series sum of label momenta and fields which gives us basic building blocks of SCET. Neglecting the power correction O(λ 2 ), we match the current J µ (x) =ψγ µ ψ(x) in Eq. (3.3) onto the operators in SCET, where n 1 , n 2 are unit four vectors in two directions and those are to be aligned with n B and n J during their sum andp 1,2 are associated label momentap i = ω i n i 2 +p ⊥ i . The matching coefficient C carries spin indices α, β. The SCET operator O after BPS field redefinition [40] is given by Here we just write the quark part and drop the gluon part, which is not relevant in DIS [29]. The quark jet fields χ n,p (x) are composed of n-collinear quark fields ξ n (x) and Wilson lines W n (x) expressed as where P µ is a operator [39] that measures the label momentum of collinear fields P µ χ n,p = p µ χ n,p . When it is summed over all label momenta χ n (x) = p χ n,p (x). The n-collinear Wilson line is where A µ n (x) = p A µ n,p (x) is a n-collinear gluon field. The soft gluon Wilson line can be written in the fundamental representation, for incoming states along n i = n B , as (3.13) In case of n i = n J i, e., for outgoing states the path for Y n J , in Eq. (3.13), turns into 0 to +∞ [56,57].
After the BPS field redefinition [40], the measurement operatorτ a in Eq. (3.4) is also split up linearly into decoupled beam-collinear, jet-collinear and soft components each of them again constructed respective momentum flow operators of SCET (3.14) Henceforth we write the collinear operators only byτ J a ,τ B a dropping the c superscript. The soft part of the decomposition can be also split up into two operatorsτ S a =τ S J a +τ S B a depending on the jet and beam contributions. Similarly, the final state X in Eq. (3.4) is also decomposed into three sectors as |X = |X cB |X cJ |X S and each operator in Eq. (3.14) matches with states in each sector.
After inserting the matching result Eq. (3.9) into Eq. (3.3), the hadronic tensor has a factorized form. After having several simplifying steps as shown in App. A, we obtain To express the cross section in Eq. ). Now the differential cross-section can be written as The sum over q goes over light quark flavors q = {u, d, s, c, b} and the functions B q , J, S are quark beam, jet, and soft functions for angularity observables defined in Eqs. (A.7), (A.14), and (A.12). Jet and soft functions are computed analytically in [36] at one-loop and numerically in [58] for soft function and in [59] for jet function at two-loop. The angularity beam function is defined for the first time in this paper and we show one-loop result in Sec. 5.

Resummation in Laplace space
The cross section in Eq. (3.16) contains logarithms of τ a , which shows singular behavior in small τ a limit and spoils the convergence of perturbation theory and these logarithms can be resummed in the Fourier space [17,60] or, in the Laplace space [61,62]. Here, we work in the Laplace space, where the transformation is defined as where ν is the variable in Laplace space conjugate to τ a . After resummation, we come back to τ a space by performing the inverse transformation : (2πi) −1 ν 0 +i∞ ν 0 −i∞ dν exp(ντ a )G(ν) Because the resummation procedure is pretty standard, in the beginning of the section we make a quick summary of resummation procedure and for those who want to follow each step closely, we give more details in following subsections. The cross section in Eq. (3.16) written as convolutions of beam, jet and soft functions in the momentum space is simply written as products of the functions in Laplace space denoted as B, J, S where σ q is simplified cross section and we should take the flavor sum q and multiply by the born cross section dσ 0 /(dxdQ 2 ) to obtain Eq.
where evolution kernels K G and η G are integration of the anomalous dimensions in Eq. (4.20), j G are constants in Eq. (4.14), the characteristic logarithm L G contains ν in its argument as shown in Eq. (4.12) and the function G(ν, µ G ) contains single and double log terms L G , L 2 G at order α s . Inserting Eq. (4.3) into Eq. (4.2) gives resummed cross section in Laplace space. To make inverse transformation easy, we replace the log terms L G in the function G(ν, µ G ) in Eq. (4.3) by a derivative operator ∂ η G /j G while we retain the log term on the exponent so that the derivative turns into the log when it hits the exponent.
where g = {b,j,s} are the functions in terms of the derivative operator. This replacement makes these functions independent of variable ν hence unchanged under the inverse transformation. The evolution term exp[j G η G L G ] ∼ ν −η G only changes as in Eq. (4.27). After rewriting σ(ν) as products ofb,j,s and collecting all L G terms from evolution factors in each function we havẽ where {µ i } means µ H , µ J , µ B , µ S and κ and Ω are sums of the evolution kernels Taking the inverse Laplace transformation we can write the cross-section in the momentum space as Up to including the flavor sum and the born cross section, Eq. (4.8) is the resummed version of Eq. (3.16). A cumulative cross-section is another conventional way to express the resummed results In rest of the section, we discuss more details about Laplace transformation and resummation. Readers familiar with technical step may skip following subsections.

Renormalization group evolution
In this subsection, we discuss RGE in the Laplace space and its solution that evolves functions of factorization from one to another scales and resum large logarithms. The RGE of S, J, B in the Laplace space and also the hard function H can be written as where G = H, S, J, B and γ G is their anomalous dimensions. Although the hard function H is not transformed it is included in Eq. (4.10) since its RGE structure is the same as other functions. The jet and beam functions are defined by the same collinear operator and their anomalous dimensions are the same γ J (µ) = γ B (µ). We just give the anomalous dimension of the jet function and do not separately give that of the beam function. The γ G takes following structure where Γ cusp (α s ) and γ G (α s ) are the cusp and non-cusp anomalous dimensions. The characteristic logarithm L G is defined as The consistency relation followed by scale independence of cross section dσ(µ)/dµ = 0 is given by γ H (µ) + γ S (µ) + 2γ J (µ) = 0, which is valid for any values of Q, µ, ν in Eq. (4.11) and it turns into three consistency relations The constants j G and κ G are given by The universal cusp anomalous dimension Γ cusp (α s ) and non-cusp anomalous dimension γ G (α s ) are expressed in powers of α s as where Γ n are given in App. D and one-loop result for γ G n are given in [36] which again satisfies the consistency in Eq. (4.13) at the order α s . The two-loop hard anomalous dimension is well known [61,63] and available up to three-loops [64] γ The two-loop soft anomalous dimension is computed in [58] and a simplified expression for the angularity is given in [59]. , , (4.18) where ∆γ C A ,n f vanishes in the thrust limit (a → 0) and γ S 1 becomes that of thrust. The one for jet function is given by the consistency . The solution of RGE is given by Eq. (4.3) where the terms K G , η G on the exponent are defined by integration of the anomalous dimension in Eq. (4.11) where we split L G (µ ) in Eq. (4.11) into two logs above and then define three integrals after replacing dµ µ by dα s /β(α s ) In App. D we give results of integration done order by order in α s . By collecting K Γ and K γ G together as below

Fixed-order coefficients
In this subsection we write fixed-order structure and coefficients of functions in factorized cross section at O(α s ). Hard function in [29] and jet and soft functions in [36] have the same logarithmic structures, which can easily be reconstructed by using RGE in Eq. (4.10). We can perturbatively solve the REG order by order in α s and truncate terms higher-order than O(α s ) In other words, we take the scale µ G in Eq.
The fixed-order beam function further factorized into short-distance coefficient and the PDF is given in Sec. 5. The coefficients j G , κ G , γ G 0 of logarithmic terms are given in Sec. 4.1 and the one-loop constant c G 1 are known from [65,66] for H and [36] for J and S.
where the jet function constant contains the integral f (a) Although two-loop constants [59,61,63,67] of hard, soft, and jet functions are known, we obtain the beam function at the one-loop level and at this moment we can achieve NNLL accuracy, for which all ingredients listed in Table 1 in App. D are available. Now we can obtain resummed functions by replacing G(ν, µ 0 ) in Eq. (4.3) by its fixed-order expression Eq. (4.22) and the resummed cross section in the Laplace space as shown in Eq. (4.2).

Inverse transformation
Let us discuss the inverse transformation to go back momentum space. We first discuss the transformation of individual functions G are essentially similar to that of the cross section. Then, We remark on the change in the case of cross section.
As discussed in the beginning of Sec. 4 we rewrite the fixed-order expression Eq. (4.22) in terms of derivative operators as where g =j,s represents jet and soft functions and analogous beam functionb is given in Eq. (5.6). Then, RG evolved result equivalent to Eq. (4.3) is where ν dependency only appears through L G in Eq. (4.12). For a moment, we use shorten conventions K G and η G for K G (µ G , µ) and η G (µ G , µ) and recover at the end of this subsection. With the following identity of inverse transformation The functions in momentum space is given by where the logarithm in τ a is defined by The cross section is products of beam, jet, and soft functions expressed in a form of Eq. (4.26) then, η G in Eq. (4.27) is replaced by the sum Ω = η S + 2η J . Finally we obtain Eq. (4.8).
We can further move all the exponents in Eq. (4.28) in front of g(∂ η G ) by shifting the derivative by where scales entering the function g and L G are µ G , while both µ and µ G enter K G and η G . The operator applied to the Gamma function turns into poly-logarithms , where ψ(x) = Γ (x)/Γ(x) and ψ (1) (x) = dψ(x)/dx. Note that in cross section, Γ(η G ) on left side of Eq. (4.31) replaced by Γ(Ω) hence, all η G on right side of Eq. (4.31) should be replaced accordingly. We finally obtain the cross section as . (4.32)

Angularity beam function
In this section we review structure of beam function, summarize our result for angularity beam function at O(α s ) and show numerical results at different values of a. We concentrate to the region a < 1 as we are restricted in SCET I region. We perform the details of computation in App. B and those who likes to know details may read the appendix.
The angularity beam function is the collinear matrix elements with initial proton state and it measures angularity τ a as well as momentum fraction x of off-shell quark having hard scattering. If the scale of beam function Qτ where the symbol ⊗ represents convolution integral defined as a ⊗ b = 1 x dξ/ξ a(x/ξ) b(ξ). We do not compute the gluon beam function B g/P (τ a , x) which is not relevant in the process. In the Laplace space it has similar structure Here, we give the expression of matching coefficient in Laplace space and the expression in momentum space is given in App. B. To the order α s , we have where leading-order coefficients 1 qj = δ qj δ(1 − z) and I (1) qj represents the one-loop correction. Structure of the matching coefficient can be obtained by solving RGE in Eq. (4.10) at fixed order in α s as done in Eq. (4.22). But unlike soft and jet functions, one has to take the factorized form in Eq. (5.2) into account when solving the RGE in Eq. (4.10). Then, we obtain where C qj = C F , T F for j = q, g. One of the logarithmic terms L B is associated with PDF with the splitting functions P qj In Eq. (5.4), the logarithmic terms are known since given Eqs. (4.14) and (4.16) whereas the constant c qj 1 (z) is new result and not given in any literature to our best knowledge. The beam function in terms of the derivative operator similar to Eq. (4.25) is given by where the distribution function L n (1 − z) is defined as The constants c qq 1 , c qg 1 with a = 0 reduces to those in thrust beam function [68]. Note that our factorization is valid in a < 1 region. In the region, they are well-defined and smooth functions of a. Their momentum-space expressions are given at the end of App. B. Now we present the numerical results of our angularity beam function. We compare the cumulative beam function that captures delta function δ(τ a ).
In Fig. 1 we show resummed results at LL, NLL, and NNLL accuracy as well as fixed-order result at NLO 5 . For resummed result we evolve the function from canonical scale µ B = Qτ to a hard scale µ = Q and for the fixed-order we take µ = Q. 5 Here, by NLO we mean fixed-order results at O(αs) from the factorized cross section in Eq. (3.16), which captures all singular terms of full QCD cross section. The full QCD result is not included in this work and will be considered in future work.  To estimate perturbative uncertainties the scale µ B is varied by a factor of two up and down for resummed result while the scale µ is held fixed. In the fixed-order the scale µ is varied by factors of two. We use MMHT2014 PDF data set [69] for the plot. The results are shown for u quark at Q = 30 GeV. The rows are for different values of angularity parameter a = −0.5, 0, and 0.5 and the two columns are corresponding to the different x = 0.05, 0.2. The numerical results show, the shape of τ a distribution is sensitive to a while less sensitive to x. The second row indicates the thrust results (a = 0 limit) for angularity beam function. Perturbative uncertainties tends to smaller at x = 0.05 due to typical pitching behavior of PDF around this value. At NLO, the beam function is related to fragmenting jet function (FJF) [70][71][72] via crossing symmetry at the amplitude level and both are factorized in similar way in terms of perturbative matching coefficients and non-perturbative (NP) matrix elements: PDF and fragmentation functions. Therefore, their matching coefficients I ij (τ a , z) and J ij (τ a , z) share many terms in common. It would be interesting to study a quantitative comparison coefficients between beam function and FJF.
The FJF includes the radiation around the measured fragmenting hadron at the final state. At first order in α s in fragmenting jet an off-shell quark splits into two partons as i * → kj whereas, in case of beam function an on-shell particle provides radiation i → k * j.
One can obtain one from the other by using following crossing relation between splitting functions [72] P where C qj = C F , T F for j = q, g and the splitting function P qj (z) is given in Eq. (5.5).
All the other terms in the beam function and FJF are cancelled out. Note that, the difference between the coefficients in the momentum space is the same as ∆ qj (z, a) up to multiplication of δ(τ a ). The differences is sensitive to a and z only and remains constant with respect to τ a . To understand the impact in the beam function we define a relative change D q normalized by the leading-order beam function, which is simply the quark PDF

Angularity distributions at NNLL
We now present our prediction for the angularity distribution in Eq. (3.16), which resums large logarithmic terms to NNLL accuracy by using formula given in Eqs. (4.8) and (4.32) and α s coefficients of hard, jet, soft, and beam functions in Sec. 4 and Sec. 5. We are interested for the angularity parameter space well below a = 1 and do not go to the further a → 1 where power correction to the SCET I factorization cannot be neglected for such value of a. Ref. [43] discusses accessibility of this region of a, using recoil-sensitive factorization which is beyond scope of the work. Note that the cross section is differential in x and Q 2 as well as in τ a and we will show the dimensionless cross section normalized by the born cross section in Eq. (A.22) where dσ q /dτ a is given by Eqs. (4.8) and (4.32). As the differential cross-section falls rapidly in small τ a region, we also multiply by a weight factor τ a for better visibility. We choose √ s = 140 GeV as EIC plans to achieve √ s = 45 and 140 GeV [8,9]. The coupling constant is taken as α s (m z ) = 0.118 to have consistency to the central value of the MMHT2014 PDF data set [69]. Throughout this paper we use the PDF data set MMHT2014 at NNLO and include five quarks and antiquark flavors and three-loop beta functions are used for consistency. Alternative choices such as NNPDF31 [74], CT18 [75], HERAPDF20 [76], ABMP16 [77] show the similar result and their differences are much smaller than our perturbative uncertainties at NNLL accuracy. As shown in the analytical form Eq. (4.8), the angularity differential-cross section is sensitive to the parameters longitudinal momentum fraction x, momentum transferred in the process Q 2 and the angularity parameter a. The sensitivity to the distribution on these parameters x, Q 2 and a are studied numerically.
In the numerical calculations we estimate the perturbative uncertainty generated by varying the scales µ H,S,J,B given by the profile functions presented in e + e − angularity [59]. The profile functions are given in App. C. To estimate the uncertainties we specifically vary three parameters in the profile function Eq. (C.2) by a factor of 2 as e H = 2 ±1 , e S = ±1/2 and e J,B = ±1/2. It is noticed that, the profile function in [59] is designed for Z-pole or higher scale and not so suitable for low Q and large positive a region-show discontinuity in µ s , µ J,B for Q < 25 GeV and a > 0.3 region. We discuss a possible modification in the profile function to access the the low Q and positive a in App. C. This modification provide access to the low Q ∼ 10 GeV and large angularity parameter a = 0.5 region. For fixed-order uncertainty, µ is set to be Q and varied up and down by a factor of two.     figure contains four plots: the NLO result with perturbative uncertainty is illustrated in gray and the resummed results at the LL, NLL and NNLL accuracy are shown in green, blue, and red curves with uncertainty bands respectively. The difference between the NLO and NNLL results shows the effect of the resummation. The whole cross-section can be characterized in three physical regions: the peak region (τ a ∼ 2Λ QCD /Q 1), the tail region (2Λ QCD /Q τ a 1) and the far-tail region (τ a ∼ 1). In the peak region, nonperturbative effect is dominated and our perturbative approach is invalid. One can adopt a NP model [29] to take it into account. Since we do not include the NP effect and corresponding region τ a ∼ 0.01 in Fig. 3 should be modified significantly. In the absence of NP effect, one sees that NLO result blows up due to singular terms ln(τ a )/τ a , 1/τ a , while the LL, NLL, NNLL results resum those terms and well-behave in the region. In the tail region, we find the resummation effect is still significant by comparing unresummed NLO and resummed NNLL. The deviation of NLO from NNLL implies that conventional scale variations of NLO underestimates perturbative uncertainty of a pure fixed-order result in this region. We also find a reasonable overlapping between resummed results that imply perturbative convergence from LL to NNLL. The NP effect in the region can be estimated as a power correction of a NP matrix element as shown in [78] and recently it is proposed that the power corrections can be determined in the framework of Eikonal Dressed Gluon Exponentiation [79]. In the far tail region, the singular term are not dominated and resummation effect is not significant hence the nonsingular terms of α s , which are neglected in our calculation will give significant contributions. The region is not displayed in Fig. 3. The sub-figures in Fig. 3 indicate that the peak moves to right with the increasing value of a and also the peak value increases. The results for a = 0 represents the differential cross-section for DIS 1-jettiness with jet axis in [29]. For better understanding of uncertainty sensitivity on the parameters a, we present a relative percentage distributions in Fig. 4. The perturbative uncertainties increase with a from 10% to 30% at NNLL accuracy in Fig. 4. The a dependency is associated with the fact that the anomalous dimension in Eq. (4.11) largely depends on a through the constants j G and κ G multiplied by the cusp anomalous dimension. The large uncertainty at positive a due to the factor 1/ (1 − a) in κ G and the singularity at a = 1 implies that our result is invalid near that value as discussed in Sec. 3. Once we manually change a dependency in those constants while retaining the consistency in Eq. (4.13) still valid, we observed that the a dependency in uncertainties accordingly changes. Therefore, the anomalous dimension sensitivity to a leads to strong a dependency in the perturbative uncertainty.
We also present the angularity cross-section for fixed a = −0.5 and different x = 0.05, 0.5 in Fig. 5, where the colors code are same as the colors code discussed above for Fig. 3. The figures show that perturbative uncertainty increases with the increasing value of x. We also compare the uncertainty at NNLL for different a = −0.5, 0, 0.5 and low x = 0.05 as well as for the high x = 0.5 in Fig. 6. The uncertainty (at NNLL) is minimum and less sensitive to a for low x = 0.05 unlike higher x = 0.5.

Conclusions
In DIS, events with one final jet with initial state radiation dominate in final states at typical values of x. These events can be studied by using event shapes. Angularity would be useful in DIS because it allows us to access a class of event shapes interpolating between thrust and broadening and extrapolating beyond the region by varying an angularity parameter a, and to perform event-shape analysis systematically in the a space.
In this paper DIS angularity τ a is defined in terms of axes which take into account one-final jet and initial-state radiation along the proton-beam. The value of τ a is small for such events and one can describe the region using soft and collinear degrees of freedom. By using SCET we drive a factorization formula in a < 1 region and it is composed of hard, beam, jet and soft functions. The jet and soft functions are the same as those appearing in e + e − angularity and hard function appears in typical SCET factorization such as endpoint region or 1-jettiness in DIS.
Our factorization contains an angularity beam function, which is new result and we give the expression at O(α s ). The beam function has logarithmic terms which are known from jet function anomalous dimension and newly presented constant terms, which are functions of a as well as momentum fraction x and are smoothly changing with a. Our beam function at a = 0 reproduces well-known thrust beam function. We compare the beam function with FJF which shares many terms in common. The difference between their perturbative matching coefficients is expressed in a single constant term proportional to the splitting function. We find that in positive a region the difference increases relatively fast compared to negative a region due to a contribution from FJF while the contribution from beam function is steady in both regions. We also show τ a distributions and perturbative uncertainties of beam function and find that they are sensitive to the value of a while its generic behavior is similar to that of thrust beam function.
We also resum large logarithms of angularity by renormalization group evolution of the functions in our factorization and give our prediction for DIS angularity at next-tonext-to-leading log accuracy. We adopted scale profile function used in e + e − angularity and estimated perturbative uncertainties by varying three scales in hard, jet, beam, and soft functions separately. We find the distributions are sensitive to the values of a and per-turbative uncertainties tends to be smaller in negative a region but not further decreasing a around −1 and below.
Studies of this paper include the factorization, beam function and resummed cross section for DIS angularity. For application to future experiments such as EIC and/or EICC, there are several things to be considered in the future. Our results contains singular terms and are correct in small τ a region. Nonsingular terms obtained from fixed-order cross section will improve accuracy in large τ a regions and also be useful to determine the scale profile function, which at this moment we borrow from that of e + e − angularity. Non-perturbative effect to angularity should be included in future studies. For higher resummation, all ingredients except for two-loop constant of beam function are known to achieve so called NNLL , which includes two-loop hard, beam, jet, and soft functions on top of NNLL accuracy and three-loop non-cusp anomalous dimension of soft or, jet function is necessary to achieve N 3 LL accuracy.

Acknowledgments
The work of D.K. and T.M. is supported by the National Natural Science Foundation of China (NSFC) through Grant No. 11875112 and the China Postdoctoral Science Foundation through Grant No.KLH1512104. All authors contributed equally to this work.

APPENDIX A Intermediate steps of factorization
In this section, we summarize intermediate steps where we perform summing over label indices, applying label-momentum conservation, and simplifying soft and collinear matrix elements to get the final form of factorization formula in Eq. (3.16). Now using matched QCD current onto SCET in Eq. (3.9) and decomposition ofτ a in Eq. (3.14), the hadronic tensor of Eq. (3.3) can be written explicitly along with the matching currents and matching coefficients as The complex conjugate to the matching coefficient C µ (p 1 ,p 2 ) = γ 0 C † µ (p 1 ,p 2 )γ 0 and indeed it only depends on Lorentz scalar expressed as C † µ (p 1 −p 2 ) 2 . The proton is considered as n B -collinear state indicated and the state can decomposed into proton state in n B sector and vacuum states associated with independent sectors or, modes. In Eq. (A.1), the collinear sectors (along n 1 and n 2 ) and the soft sectors are decoupled and it can be factored out in terms of collinear and soft matrix elements of associated vacuum states. By taking into account that the fields within each collinear matrix element are zero if they are not aligned with the collinear sector, we sum over all n 1,2 and n 1,2 and come up with n J or n B . Now absorbing the integrals overp 1,2 into the unlabeled fields χ n 1,2 and decomposing the proton matrix element into the soft and collinear matrix elements one can write the hadronic tensor as Since the Wilson lines are space-like separated and the time ordering is same as the path ordering [65,80], in the soft matrix element. Two terms in the last line come from the fact that we have two ways to chose a pair of collinear fields in the proton matrix element and they turn into quark and anti-quark beam and jet functions. Now let us perform the x integration. During the integration we can ignore x dependence in the collinear field χ n,p (x) and soft Wilson line Y n (x) because these fields describe fluctuations of low energy/momentum modes and in exponent terms like e iks·x or, e in·kn·x/2 of these fields, the momenta are of residual scale suppressed compared to label or hard momenta in the exponent in first line of Eq. (A.2).
By expressing dot product in x · q, x ·p 1 , x ·p 2 in terms of n B ,n B or, n J ,n J coordinates in Eqs. (3.5) and (3.6) Note that ω 1,2 are large component of the label momentum:p i = ω i n i 2 +p ⊥ i . Rewriting the differential Then, we havep where we dropped subscript inp 2⊥ =p ⊥ and in the first line, we set the transverse momentum to zero q ⊥ +p ⊥ = 0 according to definition of the jet axis n J which we adjust such that it is aligned with jet momentum and no transverse momentum is left inp 1 . The hadronic tensor reads as where in the second equality, the matching coefficients C depending on q 2 is moved out of the integral and transverse momentum delta function in the proton matrix element is integrated. Now one soft and two collinear matrix elements are convolved by τ a integrals and corresponding measurement functionsδ S,J,B . Matrix elements in the above expression can be related to the known hard, soft, jet, and beam functions and by using these functions we further simplify the expression. Next we express the matrix elements of Eq. (A.6) in terms of the hard, jet, soft and beam function defined in the SCET formalism. The collinear vacuum matrix elements in Eq. (A.6) can be expressed in terms of the angularity jet function [36]. The jet function integrating over + and over x in Eq. (2.17a) [36] is given by where we specify large momentum component w in the argument. The e + e − angularity defined in [36] differs from DIS angularityτ J a in normalization and for a particle of momentum p, they are related where we find w =n J · q by comparing Eq. (A.7) and Eq. (A.6).
where we use approximations w = Q 2 w J /(2q B ·q J ) and 2q B · q J ≈ Q 2 , which are correct up to power corrections of order λ 2 . Then, the jet function that measures τ a can be rewritten in terms of Eq. (A.7) after rescaling of the measurement delta function by A.
Now the collinear matrix element for vacuum state in Eq. (A.6) can be expressed as in the last step we used the one-loop expression in [36] and identified that a term of the form µ/(n J ·q) becomes . From now on, we drop flavor index q since we have just quark-jet function, which is the same for any flavor q. Similarly for the soft function we also omit the index.
For the proton matrix element we can define a angularity beam function similar to jet function definition: We Where in the last step, we used the one-loop expression in Eq. (B.3) and the scaling factor A is absorbed in the way that replaces −n B · q by Q as in Eq. (A.11). The angularity soft function was obtained in [36] S(e s , Q) It is expressed with Wilson lines along n and alongn directions with a normalization n·n = 2 while n J and n B of Wilson lines in Eq. (A.6) is different n J ·n B = 2. We rescale them and rewrite soft matrix elements in Eq. (A.6) in terms of The Wilson line Y n is invariant under re-scaling of n by a constant factor α: Y n = Y αn then, we have Y n B = Y n and Y n J = Yn . When a particle enters the hemisphere n·p <n·p the observable e s is given by While our angularity for q B ·p < q J ·p equivalently, n ·p <n ·p is given where Q R = Q 2 / √ 2q B ·q J and A = Q/Q R is the scale factor for the soft function. This is still valid for multi-particle final state.
The soft matrix element from Eq. (A.6) can be expressed in terms of angularity soft function as where we rewrote the matching coefficient . We finally obtain the expression in Eq. (3.15). The hadronic tensor in Eq. (3.15) is contracted with the lepton tensor in Eq. (3.2) where the hard function is defined as .21) and the born level cross section is given by in the last step we use approximation q J = xP + q + O(λQ).

B Beam function computation
Here we compute the constant terms c qq 1 and c qg 1 in the one-loop correction I (1) qj given in Eq. (5.4). The bare beam function for beam thrust is presented in several works [68,72,[81][82][83]. At one-loop, we just have a single parton emission from initial parton and by using relation between thrust and angularity for a single particle we can obtain the angularity beam function from known thrust beam function in [72].
In Breit frame, the initial parton with momentum P − = Q/z splits into struck quark with the large momentum component zP − = Q and the emitted gluon with (1 − z)P − = 1−z z Q. The angularity is written in terms of the gluon momentum k as Note that, we showed in App. A that this is equivalent to the definition given in Eq. (2.6) up to corrections suppressed by powers of λ. For a single particle, the angularity τ B a is related to the beam thrust t as (B.2) Now we rewriting the one-loop expression given in [72] in terms of angularity as where the function h q is given by In the Laplace space, quark beam function B q/P (ν, z, µ) is factorized into the Laplacespace coefficient I qj (ν, z/ξ, µ) and proton PDF f j/P (z), where ν is conjugate variable to τ a as in Eq. (5.2). The perturbative beam function is computed for quark or gluon state instead of proton. At one-loop, the function B q/k for a parton k = q, g can be written as By comparing Eqs. (B.5) and (5.2) one finds proton PDF f j is replaced by the parton PDF f j/k with a flavor k while the matching coefficient I qj remains the same. The parton PDFs to order α s are where 1 jk is defined below Eq. (5.3), the color factors are C qq,qq = C F and C qg,qg = T F , and the splitting functions P ik are given in Eq. (5.5). Therefore, one-loop contribution from parton can be written as The above expression is compared to renormalized beam function and the coefficients I qq,qg are determined. The term 1/τ 1+ a in the bare function Eq. (B.3) turns into Γ(− ) ν by the Laplace transformation defined in Eq. (4.1). The bare function in the Laplace space is given by After expanding Eq. (B.8) in powers of , we identify IR divergence that is matched to one-loop quark PDF in Eq. (B.6). The beam function have the same IR divergence with PDF [68] and leftover divergences are UV and they are absorbed into a renormalization factor Z qq . The finite term is the matching coefficient I qq . In the case of jet function, all IR divergences are cancelled as shown in [36]. Up to finite terms, we have The finite term at one loop is the matching coefficient with the form where the constant term c qq 1 is given by Eq. (5.7). The constant c qq 1 is newly computed for the first time while the coefficients of L 2 B and L B are the anomalous dimensions known from jet function which can be found easily comparing to Eq. (4.22).
The beam functions for gluon initial state in momentum space [72] and Laplace space are given by The bare function matches to bare gluon PDF and the matching coefficient as where the matching coefficient is given by and the constant c qg 1 is given by Eq. (5.8). For completeness we give the one-loop result in momentum space. By using the identity 1/τ (1+ ) a ↔ Γ(− )/ν between momentum and Laplace space and expanding in powers of we have (B.14) We obtain , (B.16) where L µ = ln(µ 2 /Q 2 ).

C Profile functions in high and low Q regions
For completeness, we rewrite the adopted profile functions presented for e + e − angularity in [59], which is aimed for high Q region around the Z pole and above, and needs a proper modification to apply for lower Q region. We discuss the modification to access the low Q and large positive a region in detail.
The factorized cross-section in Eq. (4.8) has four scales-hard scale µ H , soft scale µ S , jet scale µ J and beam scale µ B , each of which is related to corresponding hard, soft, jet, and beam functions. By choosing canonical scales µ H = Q, µ J,B = Qτ 1/(2−a) a , µ S = Qτ a , we can minimize the logarithms in each function at fixed-order in α s then resum large logarithms by the RG evolution from the canonical scale to a common scale µ as in Eq. (4.3). In small and perturbative region Λ QCD /Q τ a 1, the canonical choice works well. However, the scale profile as a function of τ a should be properly modified from the canonical scales in the fixed-order region τ a ∼ O(1) and the non-perturbative region τ a ∼ Λ QCD /Q. For this purpose we divide the distribution into three regions In the tail region, the canonical scales work well. As we approach the peak region we needs to freeze the scales before soft scale reaches non-perturbative region, where our prediction is not valid and as we approach the far-tail region we make all scales merge the hard scale smoothly. The profile function satisfying these constraints has the following forms [53,59] where the running scale µ run is given by The ζ function is designed to ensure the continuity from initial region {t i , y i , r 1 } to final region {t f , y f , r f }. The arguments y j and r j stands for the intercept and slope in that region. The transition between the non-perturbative, resummation and fixed-order regions of the distributions are controlled by the parameters t j as  In this work, τ sph a ranges from τ sph −1 ≈ 0.356 to τ sph 0.5 ≈ 0.616. The soft, jet and beam scales µ S , µ J,B merge with the hard scale µ H to have matching with the fixed-order perturbation theory in the far-tail region. The a-dependence in the t 0,1,2 are chosen based on the empirical observation, location of the peak scales as 3 a and t 2 is chosen as a point at which the singular and nonsingular contribution become equal in magnitude. The values of t j would be different in DIS and needs to be determined by using an fixed-order result for DIS angularity. However it is absent in our paper and we use the same values used e + e − angularity [59].
Theoretical uncertainty of our prediction is probed in band method [59], where the parameters n 0,1 are taken as constant n 0 = 1 GeV, n 1 = 10 GeV. We chose n 0 = 1 GeV ensuring a comparable deviation from the canonical scale at point t 0 . In Eq. (C.3), r = 1 and µ 0 = 1 GeV and explicit form of the ζ function is given [53,59] as where the coefficients are This profile function shows reasonable results at high Q > 30 GeV. According to the white paper, EIC is going to cover a large range in Q and we may need a access low Q  Table 1. Resummation accuracy and corresponding order of individual ingredients: cusp and non-cusp anomalous dimensions, beta function, and fixed-order hard, jet, beam, soft functions region experimentally. We noticed that this profile function shows discontinuity in µ s , µ J,B for the low Q and positive a region. The reason for this discontinuity is as follows: from Eq. (C.4) one can easily see that at low energy and positive a, t 1 becomes larger than t 2 , t 3 and the ζ-function in µ run (τ a ) failed to provide continuity at the transitions. Note that this discontinuity varies with a as t 1 ∝ 3 a indicates, if a increases the point of discontinuity moves towards large τ a (one can easily see by plotting the profile functions for Q ≤ 30 GeV and a ∼ 0.5).
To have access to the lower Q region we modify the above profile function of Eq. (C.3) incorporating the following conditions: if t 1 ≥ t 2 , the t 2 is set to be equal to t 1 and so on for the t 3 . The result for this modification is shown in Fig. 7, where first second and third columns are corresponding to Q = 60, 30, 10 GeV and two rows are for a = −0.5, 0.5. Each sub-plots represents the modified jet scale µ J in orange continuous line and soft scale µ S in blue dashed line. The colored (orange and blue) bands represent uncertainty in the corresponding scale variation. These scale variation is operated by varying e H = 2 ±1 , e S = ±1/2 and e J,B = ±1/2 in Eq. (C.2). To have a better visibility, we do not show uncertainty in the hard scale as it simply varies with single power in Q. Fig. 7 shows that our modification in profile function provides access to minimum scale Q min ∼ 10 GeV and a = 0.5.

D Anomalous dimensions and related integral
Here, we summarize our convention for logarithmic accuracy and give the expressions of anomalous dimensions and beta function coefficient used in the resummation.
Solution of RGE in Eq. (4.10) contains the exponent of integrated anomalous dimensions as shown in Eq. (4.19), which essentially resums logarithms. Here, our power counting for a large log L is α s L ∼ O(1) and leading log is a term like α n s L n+1 which is of order L ∼ 1/α s . Then, next-to-leading log (NLL) is like α n s L n and N k LL is like α n+k s L n . With this log counting, if we insert 1-loop cusp result ∝ α s Γ 0 into Eq. (4.19), we can find that the leading log term α s ln 2 (µ/µ G ) at order α s is obtained without α s evolution and all leading log terms are captured correctly with α s evolution with 1-loop beta function coefficients ∝ α s β 0 . In RGE, the non-cusp term is suppressed by one power of log and begin to contribute from NLL accuracy. The fixed-order functions with no large log have ordinary α s counting then, α n s term is comparable with N n+1 LL accuracy. Table 1 summarizes ingredients needed at each log accuracy.
To the NNLL order we need [84,85], The beta function expanded in powers of α s is given by The coefficients [86,87] are where X ≡ 1 + α s (µ 0 )β 0 ln(µ/µ 0 )/(2π). In our numerical calculations we take the full NNLL results in Eq. (D.4) for K Γ,γ , η Γ and in Eq. (D.5). To be consistent with the value of α s (µ) we take the NNLO PDFs in our numerical results.