Analytic Calculation of 1-Jettiness in DIS at $\mathcal O(\alpha_s)$

We present an analytic $\mathcal O(\alpha_s)$ calculation of cross sections in deep inelastic scattering (DIS) dependent on an event shape, 1-jettiness, that probes final states with one jet plus initial state radiation. This is the first entirely analytic calculation for a DIS event shape cross section at this order. We present results for the differential and cumulative 1-jettiness cross sections, and express both in terms of structure functions dependent not only on the usual DIS variables $x$, $Q^2$ but also on the 1-jettiness $\tau$. Combined with previous results for log resummation, predictions are obtained over the entire range of the 1-jettiness distribution.


Introduction
In high energy colliders, jet production plays an important role in probing the strong interaction, hadron structure, dense media, and new particles beyond the Standard Model. Thus predicting jet production cross sections and jet structure is one of the important tasks of Quantum Chromodynamics (QCD). Jet algorithms [1][2][3][4][5][6] allow exclusive study of jets and definitions of cross sections with a definite number of jets. However, they also introduce various parameters like jet radii or sizes and jet vetoes, which require more effort to predict accurately in analytic calculations in QCD. Event shapes [7] provide a simple, inclusive way to identify final states that are jet-like, and can often be predicted to very high accuracy in QCD. Thrust in e + e − collisions [8] is a classic example of a two-jet event shape that has been extensively studied in both theory and experiment. Thrust cross sections in e + e − have been predicted to very high accuracy, N 3 LL+O(α 3 s ) in resummed and fixed-order perturbation theory [9][10][11][12][13][14], along with rigorous treatments of nonperturbative power corrections [14][15][16], that have led to unprecedented 1%-level precision in determinations of the strong coupling constant α s from e + e − event shape data [13,14,17].
Event shapes in DIS have also been studied but not as extensively as in e + e − , and the theoretical accuracy has yet to catch up to the same level. Two versions of DIS thrust have been defined and measured in H1 and ZEUS experiments at HERA [18][19][20][21][22][23] and they have been calculated up to next-leading-logarithmic accuracy (NLL) at resummed order and numerically to O(α 2 s ) at fixed order [24,25]. The measured DIS thrusts involve non-global logarithms (NGLs), which present a theoretical obstacle to higher order accuracy [25,26].
Versions of thrust such as e + e − thrust and the DIS thrust τ Q defined in [24] do not suffer from NGLs. A class of event shapes called N -jettiness τ N [27] is a generalization of these versions of thrust and are applicable in different collider environments, including e + e − , leptonhadron, and hadron-hadron collisions. τ N measures the degree of collimation of final-state hadrons along N light-like directions in addition to any initial-state radiation (ISR) along the incoming beam directions. In a number of recent papers [28][29][30], factorization theorems for various versions of 1-jettiness τ in DIS have been derived by using soft collinear effective theory (SCET) [31][32][33][34][35]. To date, this has enabled log resummation up to NNLL accuracy [28][29][30], which is one order higher in resummed accuracy than earlier results [24,25].
The SCET results [28][29][30] correctly capture and resum all logarithmic terms (singular), while non-logarithmic terms (nonsingular) can be obtained from fixed-order computations in full QCD. The full cross section is the sum of singular and nonsingular parts and can be written as σ full (τ ) = σ sing (τ ) + σ ns (τ ) . (1.1) The singular part is factorized in terms of hard, jet, beam, and soft functions each of which depends on the relevant energy scale for each mode [28][29][30]. This separation of scales and renormalization group (RG) evolution between them allows for resummation of the large logarithms in the fixed-order expansion of the cross section. When the RG evolution is turned off in the singular part, the full cross section reduces to the ordinary fixed-order result. The nonsingular part is obtained by subtracting the fixed-order singular part from the fixed-order cross section. For an accurate prediction over the entire range of an event shape distribution, both fixed-order and resummed calculations should be consistently improved. While NNLL resummation of the singular part in Eq. (1.1) has been performed for three different versions of DIS 1-jettiness in [28] and another version in [29,30], no analytic computations of the nonsingular part at O(α s ) or above have yet been performed. In [36] an O(α s ) result has been numerically obtained for a version of 1-jettiness that requires a jet algorithm to determine the jet momentum. Such a numerical approach is appropriate for such cases and allows for the flexibility of using different jet algorithms.
In this paper, we carry out the first analytic O(α s ) calculation for a DIS event shape. We choose the version of 1-jettiness called τ b 1 in [28], which groups final-state particles into back-to-back hemispheres in the Breit frame and is the same as the DIS thrust called τ Q in Ref. [24]. It can be written as where p i is the momentum of the ith particle in the final state, and Q 2 ≡ −q 2 is determined by the momentum transfer q in the event. The reference vectors are defined by q b B = xP and q b J = q +xP , where P is the proton momentum. In the Breit frame these vectors point exactly back-to-back. The second definition in Eq. (1.2) is valid in the Breit frame, and requires measuring the z components p i z of momenta of particles only in the jet hemisphere (current hemisphere) H J . The definition in Eq. (1.2) differs from the measured version τ H1 = 1−T ZEUS γ [20,23] in normalization (replacing 2/Q by 1/E hemi where E hemi is the hemisphere energy). We present our results in terms of fixed-order singular and nonsingular parts of the cross section as in Eq. (1.1). They can be put in a simple form which can easily be implemented in other analyses. The main new results of this paper are the nonsingular 1-jettiness structure functions given by Eq. (4.5).
We also show numerical results with perturbative uncertainties by varying scales at the HERA energy. Our results could be compared to existing HERA data [18][19][20][21][22][23] or to future EIC data [37]. In [28], by comparing our resummed singular cross section to the known fixed-order total cross section, we estimated that the nonsingular corrections would amount to several percent of the total cross section, and this expectation is borne out by our computations here.
The paper is organized as follows: In Sec. 2 we briefly review the relevant kinematic variables in DIS and our definition of 1-jettiness, and express the cross section in terms of structure functions. In Sec. 3, we outline the basic steps of the O(α s ) computation including the phase space for 1-jettiness and perturbative matching of the hadronic tensor onto parton distribution functions (PDFs). Sec. 4 contains our main results, analytic O(α s ) expressions for 1-jettiness structure functions. Details of the fixed-order calculation are given in App. A, App. B and App. C. In Sec. 5 numerical results are given for structure functions at O(α s ) fixed-order accuracy and cross sections at NLL + O(α s ) resummed accuracy. Basic details entering the resummation of the singular terms are reviewed in App. D and App. E for convenience. Finally, we will conclude in Sec. 6.

1-Jettiness in DIS
In this section we review DIS kinematic variables that will be used throughout the paper and the definition of the 1-jettiness τ cross section in DIS, whose computation will be the main prediction of our paper.

Kinematic Variables
In DIS, an incoming electron with 4-momentum k scatters off a proton with momentum P by exchanging a virtual photon 1 with a large momentum transfer q = k − k , where k is the momentum of the outgoing electron. Because the photon has spacelike momentum it has a negative virtuality, and one can define the positive definite quantity (2.1) Q sets the momentum scale of the scattering. We will be interested in hard scattering, where Q Λ QCD . A dimensionless quantity x called the Björken scaling variable is defined by which ranges between 0 ≤ x ≤ 1. Another dimensionless quantity y is defined by y ≡ 2P ·q 2P ·k , which ranges between 0 ≤ y ≤ 1. This variable y represents the energy loss of the electron in the proton rest frame. The three variables x, y, and Q 2 are related to one another via Q 2 = xys, where s = (P + k) 2 is the total invariant mass of the incoming particles. The total momentum of the final state X is p X = q + P and the invariant mass is given by For large x very near 1, the final state consists of a single tightly collimated jet of hadrons. This region has been analyzed in SCET in, e.g., [38][39][40][41][42]. We will instead be interested in different region where two or more energetic jets can occur. This occurs in the "classic" region where x has a generic size x ∼ 1 − x ∼ 1 such that p 2 X ∼ Q 2 . Although the cross section we compute is frame independent, there is a convenient frame in which to perform the intermediate steps of the calculation. This is the Breit frame, where the virtual photon with momentum q µ is purely spacelike, and collides with the proton with momentum P µ along the z direction. In this frame the virtual photon and the proton have momenta where n z = (1, 0, 0, 1) andn z = (1, 0, 0, −1).

1-Jettiness
To probe the number of jets in the final state produced at a given value of x and Q, an additional measurement needs to be made. A simple event shape that accomplishes this is the N -jettiness [27], a generalization of the thrust [8]. It is defined by the sum of projections of final-state particle momenta onto whichever axis is closest among N jet and N b beam axes, where N b = 0 for e + e − collisions, 1 for ep DIS, and 2 for pp collisions. The N -jettiness τ N is designed so that it becomes close to zero for an event with N well-collimated jets in the final state away from any hadronic beam axes. For example, 1-jettiness in DIS is defined by one jet and one beam axis: where q B , q J are lightlike four-vectors along the beam and jet directions. It is natural to choose q B along the proton direction. One can consider several options for choosing q J . In [28], we defined three versions of 1-jettiness τ a 1 , τ b 1 , and τ c 1 distinguished by different choices for q J : (a) q a J aligned along the jet axis determined by a jet algorithm, (b) q b J along the z axis in the Breit frame, and (3) q c J along the z axis in the center-of-momentum (CM) frame. In this paper we consider τ b 1 for which q b B and q b J are given by As shorthand, we drop both superscript and subscript in τ b 1 throughout the remainder of the paper.
In the Breit frame, the vectors q b B,J point exactly back-to-back with equal magnitude: and divide particles in the final state into two equal hemispheres. One is the "beam" or "remnant" hemisphere H B in the −z direction and the other is the "jet" or "current" hemisphere H J in the +z direction. The 1-jettiness τ in Eq. (2.6) has an experimental advantage in that it can be determined by measuring only one of the hemispheres, namely H J . This avoids having to measure the whole final state including the beam remnants, a technical difficulty in experiments such as H1 and ZEUS at HERA. By using q and P in the Breit frame in Eq. (2.3), the 1-jettiness can be written in the form The definition Eq. (2.8) directly corresponds to the thrust τ Q in DIS defined in [24]. We can obtain the physical upper limit on τ using the kinematic constraints that the jet momentum p J z ≥ 0 has to be positive, and that the beam momentum's z component is negative, so that . These conditions imply the upper limits on τ : (2.9)

1-Jettiness Cross Section
The 1-jettiness cross section can be expressed in terms of leptonic and hadronic tensors: where the lepton tensor for a photon exchange is given by where k and k are incoming and outgoing electron momenta and α ≡ α em . The hadronic tensor is the current-current correlator in the proton state, whereτ is a 1-jettiness operator that measures 1-jettiness when it acts on the final states, which we defined in [28], based on the construction of event shape measurement operators from the energy-momentum tensor in [43][44][45][46]. In this paper, we consider only the vector current J µ = f Q fqf γ µ q f . Previously we worked with both vector and axial-vector currents, see [28] for the appropriate generalizations. 2 Because the hadronic tensor depends only on the two momenta P and q, it can be decomposed into products of tensors constructed with g µν , P µ , q µ and structure functions depending on x, Q, and τ . In our conventions, where the two tensor structures that appear are: which arise from parity conservation and the Ward identity q µ W µν = W µν q ν = 0. If we considered parity-violating scattering, e.g. with neutrinos, a third tensor T µν 3 = −i µναβ q α P β would also appear.
In terms of the structure functions appearing in Eq. (2.13), the cross section Eq. (2.10) can be expressed where F L ≡ F 2 − 2xF 1 . We use calligraphic font for the structure functions in the differential τ cross section. We will use Roman font F 1,L for the structure functions in the integrated cross section, see Eqs. (4.1) and (4.2).
The structure functions F i can be obtained by contracting the hadronic tensor with the metric tensor or the proton momentum P µ : We choose to always work with expressions in D = 4 − 2 dimensions for the vector indices µ and ν, so the factor of 1/(1 − ) in F 1 comes from taking the contraction g µν T µν 1 . The contraction P µ P ν W µν turns out to be finite as → 0. The standard structure functions depend just on x and Q 2 , while those in Eq. (2.16) are additionally differential in τ . The structure functions can be written in terms of singular and nonsingular parts as we did for the cross section in Eq. (1.1). We will present singular and nonsingular parts of the structure functions in Sec. 4, from which one easily obtains the corresponding parts of the cross section via Eq. (2.15).

Setup of the Computation
In this section we outline the basic steps in the O(α s ) computation of the 1-jettiness cross section in DIS. First, we describe the standard perturbative matching procedure for the hadronic tensor onto PDFs, which allows us to compute the matching coefficients using partonic external states. Then we set up the phase space integrals in the Breit frame in which the intermediate steps of the computation are simpler. The final results are frame independent. The reader who wishes to skip these details may turn directly to the final results in Sec. 4 and Sec. 5.

Perturbative Matching
Here, we describe the matching procedure to determine the short-distance coefficients that match the hadronic tensor W µν (x, Q 2 , τ ) onto PDFs. By using the operator product expansion (OPE) the hadronic tensor can be written in the factorized form where f i/h is the PDF for a parton i ∈ {q,q, g} in a hadron h, and w i µν is the short-distance coefficient that we will determine by perturbative matching. 3 The sum over q,q goes over all light flavors f ∈ {u, d, s, c, b} at the collision energies we consider. On the left-hand side of Eq. (3.1), the superscript h specifies the hadron in the initial state. The coefficients w i µν 3 The first power correction ∼ Λ QCD /(Qτ ) in Eq. (3.1), as well as higher-order terms ∼ [Λ QCD /(Qτ )] k , are all described by the leading-order soft function in the small-τ factorization theorem [28]. however can be computed in perturbation theory using any appropriate initial state including partonic ones. This is what we shall describe in this subsection. The factorization theorem for an initial parton j is given by where the arguments are implicit for simplicity and the convolution integral over ξ in Eq. (3.2) is replaced by the symbol ⊗. By comparing Eq. (3.2) to Eq. (3.1), all implicit notations can easily be recovered. We determine the coefficients w i by computing the W j µν , which are defined by Eq. (2.12) but with a quark or antiquark j = q,q or gluon j = g in the initial state, and subtracting out the partonic PDFs, which are IR divergent and require a regulator. We perform this computation using dimensional regularization and defining PDFs in the MS scheme.
Working to O(α s ), and denoting the order α n s piece of each function by a superscript (n) , Eq. (3.2) becomes Using MS and using to regulate IR divergences, the partonic PDFs to O(α s ) are given by (see, e.g., [47] for related discussion) where the color factors and splitting functions are given by with P qq (z) and P qg (z) given in Eq. (4.7) below. There are no contributions containing the splitting functions P gq , P gq , P gg since the tensor W j µν in Eq. (2.12) we compute contains only the quark current.
For the 1-jettiness structure functions Eq. (2.16), we only need the projections −g µν W µν and P µ P ν W µν of the hadronic tensor in Eq. (3.1); hence, we obtain the projected coefficients −g µν w i µν and P µ P ν w i µν . We compute the contracted tensors −g µν W i µν explicitly in App. B. Including the factor of (1− ) coming from the tensor contractions in D dimensions, we obtain: where w q,g G are finite as → 0. For the quark tensor we consider one flavor f at a time, since it will get convolved with a different PDF for each quark flavor, while for the gluon tensor we include the sum over all flavors. The contractions P µ P ν W i µν begin at O(α s ), and take the form and are finite as → 0. Plugging these forms into the matching conditions Eq. (3.3), we find the 1/ IR divergences cancel between the PDFs in Eq. (3.4) and the computed tensors in Eq. (3.7), leaving the finite matching coefficients We compute the finite coefficients w q,g G,P explicitly in App. B, and they are given in Eqs. (B.12b), (B.20b), (B.22b), and (B.25b).

Phase Space
In this section, we evaluate some of the phase-space integrals for 1-and 2-body final states. In the partonic computation of the tensor W µν given in Eq. (2.12) or Eq. (3.1), we sum over all the possible n-body final partonic states, where the 1/s i factor is from averaging over the spins or polarizations of the initial parton j: s q = sq = 2 and s g = 2(1− ), and in the last equality we define the n-body contribution W to W . We sum over the spins or polarizations of all the final-state partons. In Eq. (3.12), is the amplitude for the initial parton j with momentum P to scatter off the current J µ and produce the final-state partons with momenta p 1 . . . p n , and the n-body phase space integration measure is given by (3.14) The n = 1 term in the sum in Eq. (3.12) is given by where the arguments of M µ are implicit. The n = 2 term is given by where the lightcone components are (p + 2 , p − 2 ) = (n z · p 2 ,n z · p 2 ). We have chosen to do the integrals using the momentum-conserving delta function in Eq. (3.16) and the mass-shell delta functions in Eq. (3.14) in such an order that the p − 2 integral is left over in Eq. (3.16) to be done last. Where p 2 and p 1 appear in Eq. (3.16), they take values given by the formulas: The 1-jettiness τ is now expressed as a function of x and v. Two particles in the final state can be assigned in four different ways to the two hemispheres, and the formula for τ (x, v) differs in each of these regions. These four regions (a) to (d) are illustrated in Fig. 1. The function τ (x, v) can be broken down into four pieces, where the two-dimensional step function Θ (i) covers each region (i) and τ (i) is the value of 1-jettiness in the corresponding region:

Analytic Results for DIS 1-Jettiness Structure Functions
In this section we present our analytic results for the structure functions in Eq. (2.16) that determine the 1-jettiness cross section in Eq. (2.15). We will present our results in terms of the structure functions appearing in the cumulative (integrated) cross section, where the F i are given by Eq. (2.16). The results for the integrated structure functions F i are more compact to write down than for F i . We give the results for the differential structure functions F i in App. C. As the cross section in Eq. (1.1) is written in terms of singular and nonsingular parts, we express the structure functions as: The fixed-order structure functions are obtained from the calculation of projected hadronic tensors in Eq. (2.16) that are calculated in App. B and App. C. The singular part of the cross section was calculated in [28]. Our main new results here are the nonsingular parts of the structure functions that are obtained by subtracting off the known singular parts from the full expressions. We will present our final expressions for the singular and nonsingular parts of F 1 and F L in Eq. (4.3) in the following form: The singular parts of these can be extracted from the singular cross section in [28], and are given in Eq. (4.8). Our main new results here are for the nonsingular parts. The functions A ns i and B ns i are given by the nonsingular parts of P µ P ν W i µν and −g µν W i µν , respectively. They are obtained by integrating the differential structure functions in Eq. (C.6). We find which is one of our main results. Here, we have defined the theta function which turns on inside the physically-allowed region 0 < τ < τ max given by Eq. (2.9) and turns off outside. The plus distribution L n (z) is defined in App. A. The standard splitting functions P qq and P qg are given by The formulas for B q and B g in Eqs. (4.5c) and (4.5d) appear to contain terms which are still divergent as τ → 0, but these divergences cancel in the sum of all terms. Formulas for B q,g given as sums of explicitly nonsingular terms can be found in Eq. (C.9). One may recognize that the Θ 0 terms in Eq. (4.5) introduce a discontinuity in the cumulative cross section at τ = 1. This feature is associated with asymmetric initial momentum in the z direction, which can give rise to an event with one of the hemispheres containing all final-state particles and the other being empty. As illustrated in Fig. 1, this occurs in regions (a) and (d), where τ takes on its maximum allowed values in Eq. (2.9), 1 for x < 1/2 and (1 − x)/x for x > 1/2. For x < 1/2, this appears at τ = 1 as a delta function in the differential structure functions Eq. (C.6) and a discontinuity in the integrated structure functions Eq. (4.5). However, this feature is not seen for x > 1/2 at τ = (1 − x)/x, because we see that the integrals proportional to Θ 0 in Eq. (4.5) go to zero for τ = (1 − x)/x, the range of integration shrinking to zero.
The singular part of the cross section has been computed in [28], from which the singular part of the structure functions can be extracted. F sing 1 is simply half of the cumulant cross section given in Eq. (174) in [28], and F sing L = 0. The singular parts A sing i and B sing i of the functions in Eq. (4.4) are given by The sum of Eqs. (4.5) and (4.8) gives the complete fixed-order O(α s ) result for the DIS 1jettiness structure functions. When we take values of τ beyond the physical maximum, where Θ 0 terms are turned off, the result reproduces the standard inclusive structure functions in x and Q 2 , which are given by (e.g. [48]) where we have defined the two functions In this section we have presented the complete O(α s ) results for the fixed-order structure functions in the DIS 1-jettiness cross section. The expressions Eq. (4.5) for the nonsingular contributions to the structure functions in Eq. (4.4) are the primary new results of this paper.

Numerical Results
In this section, we present numerical results for the structure functions F 1,L that appear in the differential 1-jettiness cross section in Eq. (2.15) and the corresponding F 1,L in Eq. (4.2) that appear in the integrated cross section Eq. (4.1). We computed these structure functions to O(α s ) in Sec. 4 and App. C. We also present predictions for the τ cross sections themselves. For structure functions, we show the fixed-order O(α s ) results for the singular part (in τ ), the nonsingular part and their sum. For the cross section, we show resummed results at NLL + O(α s ) accuracy as well as the pure fixed-order results. At this order of accuracy we have the fixed-order parts of the hard, jet, beam, and soft functions in the singular part Eq. (E.1) at the same order in O(α s ) as in the nonsingular part. 4 4 Resummation of the singular terms in the τ cross section is in fact available up to NNLL accuracy [28].
For simplicity, we choose to illustrate results only at NLL resummed accuracy in this paper (see [14,49] for definition of primed accuracy). As described in Ref. [49], formulae for resummed differential and integrated cross sections at unprimed orders of accuracy may suffer from a mismatch in the actual logarithmic accuracy achieved, depending on how the formulae are written. One can ensure that the differential distribution at N k LL matches the accuracy of the corresponding integrated cross section by differentiating the integrated cross section including the τ dependence in the scales µi(τ ). However, in the large τ ("far tail") region, Ref. [14] observed that this procedure leads to unrealistically large uncertainties, and recommends that the τ dependence in µi(τ ) not be differentiated in going from the integrated to the differential cross section. It is possible to write the differential cross section in a way that interpolates between the two approaches for small For our numerical results plotted here, we set the collision energy to be √ s = 300 GeV, which corresponds to the H1 and ZEUS experiments, and choose Q = 80 GeV. We adopt MSTW2008 PDF sets at NLO [52] with five light quark and antiquark flavors and run α s (µ) with the 2-loop beta function in Eq. (E.5) starting at the values α s (m Z ) = 0.1202 used in NLO PDFs. Fig. 2 shows the components of the fixed-order results for the structure function F 1 (x, Q 2 , τ ) in the integrated cross section, given by Eqs. (4.4a), (4.5), and (4.8), and of the F 1 (x, Q 2 , τ ) in the differential distribution, given by Eqs. (C.3), (C.5), and (C.6), at two values x = 0.2 and 0.7. We set all scales to be µ = Q = 80 GeV. In the integrated structure function F 1 , the sum of singular (dashed line) and nonsingular (dotted line) contributions give the full result (solid line). The full result approaches the total result F 1 (x, Q 2 ) (horizontal dashed line) in Eq. (4.9a) as τ approaches 1. For x = 0.2, the singular part alone undershoots the total, and and large τ , but this task does not lie within the scope of this paper. As observed in [49], equivalent accuracy between differential and integrated cross sections is in fact maintained if one works at primed orders, whether one differentiates µi(τ ) or not. Thus we will work here at NLL accuracy and evaluate the differential cross section by not differentiating µi(τ ) in the integrated cross section, see Eq. (E. 19). This avoids the potential negative issues pointed out in both [14] and [49]. Some recent progress (e.g. [50]) has been made in obtaining ingredients needed for NNLL or N 3 LL accuracy [51] for the related version of 1-jettiness τ a 1 defined in [28,30]. 0. 2.
Τ Τ max the nonsingular part makes up the difference. For x = 0.7, the singular part overshoots the total, and the corresponding nonsingular part is mostly negative. Although it is imperceptible in Fig. 2, there is actually a small discontinuity in the x = 0.2 plot at τ = 1, and the total (solid red) F 1 does not reach the full result (dashed black) until above τ = 1. We will zoom in on this feature in Fig. 4. For the differential F 1 in Fig. 2, we plot the absolute value on a log scale. The results illustrate that there is a large cancellation between the singular and nonsingular pieces in the large τ region so that the total goes to zero in this tail. This same cancellation was discussed for e + e − thrust in Ref. [14], and appears in various other cross sections that have singular and nonsingular components. The tail falls faster for larger x because τ dependence enters into PDFs in a form like f q (x(1 + τ )), as seen in Eq. (C.6), which falls faster as x increases. The overall normalization also becomes smaller for larger x due to the PDFs falling off. Fig. 3 shows the fixed-order results for the longitudinal structure function F L (x, Q 2 , τ ) for the integrated cross section, given by Eqs. (4.4b) and (4.5), and F L (x, Q 2 , τ ) for the differential distribution, given by Eqs. (C.3) and (C.6), at x = 0.1 , 0.2, and 0.7. These are purely nonsingular in τ . The plots are normalized to the total F L (x, Q 2 ) in Eq. (4.9b). Note that F L is finite at τ = 0 at O(α s ). The distribution monotonically decreases with τ . For the left plot at x = 0.1, 0.2, there is a perceptible gap from the total (straight dashed line) at τ = 1 before the curves reach the value 1. This jump is explored in Fig. 4. Fig. 4 illustrates the discontinuities in the cumulative F 1 and F L near τ max . The jump is smaller than 1% in F 1 and is about a few percent in F L . These discontinuities are reduced for increasing x and disappear at x = 1/2 and beyond. As described in Sec. 4, these discontinuities are associated with events where the jet hemisphere is empty and the beam hemisphere contains all final-state particles as seen in the Breit frame, so whole regions of phase space end up contributing to the same fixed value of τ (see Fig. 1). Such events do not occur in the observables defined in the partonic CM frame such as e + e − thrust. This discontinuity is infrared safe, and though its magnitude is very small, it is in principle measurable.  The cross section in Eq. (1.1) with all the scale dependencies made explicit in singular and nonsingular parts can be written and is given in Eq. (E.1). The singular part depends on the scales µ H , µ J , µ B , µ S associated with hard, jet, beam, and soft radiation, respectively, and the nonsingular part depends on µ ns as in conventional fixed-order results. For the full calculation, all scales should be specified.
In the region τ 1, there are large logarithms in the singular part and the logarithms can be resummed by RG evolution of the functions between µ and their individual canonical scales: (For more details on resummation of the singular part, see [28]. Basic results are reviewed in App. E.) However, µ S cannot be arbitrary small and it should freeze above the nonperturbative regime that lives below 1 GeV. On the opposite end, where τ ∼ O(1) and logs of τ are not large, the resummation should be turned off by setting all µ i ≈ Q. In [28] we used profile functions µ i (τ ) satisfying above constraints and estimated perturbative uncertainties by varying parameters in the profile functions [14,53,54]. However, the profile defined in [28] has scales away from the canonical scales when x increases. Here we use improved profiles given in App. D, which set canonical scales in the resummation region that are independent of x. Fig. 5 shows the soft scale µ S (τ ) as a function of τ at x = 0.2 and 0.7 as well as the canonical choice τ Q (dashed line).
For the central values of µ ns and its variations, we make the same choice as [14], The scales are chosen to estimate theory uncertainties from un-resummed subleading logarithms in the nonsingular part. Fig. 6 shows resummed integrated and differential cross sections at NLL + O(α s ) as well as the purely fixed-order result at O(α s ). We plot normalized cross sections defined aŝ where σ 0 = 2πα 2 [1 + (1 − y) 2 ]/Q 4 . We multiply the differential distribution by τ for ease of displaying the whole τ region.
The uncertainty bands for the resummed results are obtained by summing all scale vari-ations described in Eqs. (5.2) and (D.5) in quadrature. The uncertainty bands for the fixedorder results are obtained from varying the central scale µ = Q up and down by factors of 2. As seen in Fig. 2 the tail of the distribution becomes shorter with increasing x. Relative uncertainties about the central value are larger for larger x because of slower convergence of the perturbative corrections associated with the PDF for increasing x (as can be seen from the fact that the residual scale dependence of the PDF increases with x). Fig. 6 only includes purely perturbative results. Nonperturbative effects in 1-jettiness are power suppressed by Λ QCD /(τ Q) for τ Λ QCD /Q, and the leading power correction can be expressed in terms of a single nonperturbative parameter Ω 1 . The parameter is universal for different versions of 1-jettiness in DIS defined in [28], and even appears in the power corrections for certain jet observables in pp → H/Z + jet with a small jet radius [55]. Alternatively, a shape function that takes nonperturbative behavior into account in the nonperturbative region as well as the power correction region [56], can be used as in [28]. In this paper, we omit implementing these nonperturbative effects.

Conclusions
Events with one or more jets plus initial state radiation dominate the population of final states in DIS for typical values of x. These events can be further probed by the inclusive event shape 1-jettiness τ . Events with small values of τ contain only one non-ISR jet, while multiple jets populate the large τ region. In this paper, we obtained analytically the O(α s ) cross section for all values of τ , and combined it with NLL resummation of the singular terms at small τ to obtain results accurate over the entire range of τ . This is the first analytic calculation of a DIS event shape at this order.
We wrote the results in terms of structure functions F 1 (x, Q 2 , τ ) and F L (x, Q 2 , τ ) which generalize the usual DIS structure functions F 1,L (x, Q 2 ). We gave structure functions for both the cumulative or integrated τ distribution as well as the differential τ distribution. Our predictions for the cumulative distribution agree with the total F 1,L (x, Q 2 ) for τ > τ max .
The cumulative cross section displays an interesting feature, a small discontinuity at τ = 1, which is a consequence of asymmetric initial momentum that can lead to one of hemispheres (in the Breit frame) being empty in the final state. This does not happen in e + e − thrust defined in the partonic CM frame.
We presented numerical results with perturbative uncertainties by varying scales at the HERA energy. In general the uncertainties grow with x due to the convergence of perturbative corrections in the cross section that are connected with the PDFs through their scale dependence. The tail of the τ distribution falls off faster as x grows. The size of the nonsingular terms is consistent with our expectations from [28] where we compared the resummed singular cross section with the total QCD cross section at x, Q 2 .
Our results represent a significant improvement in precision in the prediction of DIS event shape cross sections. The groundwork is in place to go to higher resummed [51] and fixedorder accuracy which we will pursue in the near future, and bring the science of event shapes in DIS to the same level of precision as has been achieved in e + e − . These predictions can be tested with existing HERA data and future EIC data, which should yield determinations of the strong coupling and hadron structure to unprecedented accuracy.

Acknowledgments
The work of DK and IS is supported by the Office of Nuclear Physics of the U.S. Department of Energy under Contract DE-SC0011090, and the work of CL by DOE Contract DE-AC52-06NA25396 and by the LDRD office at Los Alamos. We thank the organizers of the 2013 ESI Program on "Jets and Quantum Fields for LHC and Future Colliders" in Vienna, Austria, where part of this work was performed. DK would like to thank the Nuclear Theory group at LANL and CL would like to thank the MIT Center for Theoretical Physics and the KITP at UCSB for hospitality during portions of this work. This research was supported in part by the National Science Foundation under Grant No. PHY11-25915.

A Plus Distributions
In this section, we define plus distributions that we use and collect some useful identities involving them. The standard set of plus distributions L n (z) are defined by (see, e.g., [53]) Integrating against a well-behaved test function g(z) gives the familiar rule, We also define a distribution function with two arguments, which can be used when the presence of the divergence in a variable z is controlled by the value of a second variable z 0 , where g(z) is a test function. In the standard distribution L n (z) the subtraction of the singularity occurs at the singular point z = 0, while in L(z, z 0 ) the subtraction occurs at z = z 0 even if there is no singularity when z 0 = 0. L(z, z 0 ) reduces to L(z) at z 0 = 0. Evaluating phase space or loop integrals at O(α s ) or higher in dimensional regularization, we encounter singular terms like where the O( 0 ) terms are where Li(z) is the dilogarithm, defined by The function I s (τ, z) is singular both in τ and z, depending on both L(τ ) and L(x), while I ns (τ, z) is not singular in τ (though still singular in z). Note that the term on the last line of Eq. (A.8), which has a δ − τ + z 1−z , will not contribute to any of our perturbative structure functions because the expression in brackets that it multiplies and its derivative respect to τ are both zero at τ = z 1−z .

B Hadronic Tensor at Parton Level
In this section we calculate the hadronic tensor W µν defined in Eq. (2.12) where the proton initial state is replaced with a partonic (quark or gluon) state. Such a computation allows us to extract the short-distance matching coefficients w q,g µν in Eq. (3.1) onto PDFs, as described in Sec. 3.1. We denote the tensor for a quark initial state as W q µν and for a gluon initial state as W g µν . Up to O(α s ), W q µν involves a tree-level contribution and the one-gluon diagrams in Fig. 7, and can be decomposed into Meanwhile, W g µν is given just by tree-level real diagrams at O(α s ), shown in Fig. 8. The partonic tensor W i µν can be computed from Eq. (3.12), which is a phase space integral over the squared amplitude. In this section we compute the squared amplitudes; in the next section we will evaluate the complete phase space integrals. Figs. 7 and 8 represent the O(α s ) amplitudes for initial quark and initial gluon states. In App. B.1 and App. B.2 we evaluate the squared amplitudes built from these diagrams. The real diagrams have two-body final states with momenta p 1 and p 2 and as described in Fig. 1, they can enter two back-to-back hemispheres in four different ways, and the formula for the 1-jettiness τ in terms of p 1,2 in each configuration differs.

B.1 Squared Amplitudes for γ * + q
For the process γ * (q µ )+q(P µ ) → q(p µ 1 ) at tree level, the amplitude isū(p 1 )γ µ u(P ). To obtain the structure functions Eq. (2.16) one needs projected squared amplitudes as where we have also summed over all quark spins. The projection in Eq. (B.3) is zero because of the Dirac equation P / u(P ) = 0. Here, Q f is the electromagnetic charge of quark with flavor f . We do not sum over flavors for the quark tensors until we convolve with PDFs. The virtual contribution can be extracted from the literature, see, e.g., Eq. (14.19) in [57]. At O(α s ) we obtain the cross terms between the tree-level and the virtual diagram shown in Fig. 7: again summed over all quark spins. We have kept the factor (1 − ) out front because it is to be cancelled by the same factor in Eq. (2.16). In the second step of Eq. (B.4), we converted to the MS scheme, making the replacement: Note that the finite part in Eq. (B.4) is the α s term of the hard function, which already appeared in our discussion in Ref. [28]. For discussion of the hard function in SCET see [38,58]. Eq. (B.5) is zero again by the Dirac equation P / u(P ) = 0. The real contribution to W q µν in Eq. (B.1) at O(α s ) comes from two diagrams for γ * (q µ )+ q(P µ ) → q(p µ 1 ) + g(p µ 2 ) shown in Fig. 7. The projected amplitudes for the diagrams, summed over quark spins and gluon polarizations, are given by where v = p − 2 /Q as in Eq. (3.18) or Fig. 1. Eq. (B.7) can be found from Eq. (14.23) in [57].

B.2 Squared Amplitudes for γ * + g
The tree-level process with an initial gluon γ * (q µ )+g(P µ ) → q(p µ 1 )q(p µ 2 ) starts at O(α s ), illustrated in Fig. 8. The projected amplitudes for the process, summed over gluon polarizations and quark spins, are given by Note that these are symmetric under v → 1 − v, which just switches the final-state quark and antiquark in Fig. 8. Note also that the projection in Eq. (B.10) is independent of v, making the phase-space integral in Eq. (3.18) particularly simple. Here we go ahead and include the sum over quark flavors f ∈ {u, d, s, c, b} since the gluon PDF with which we will convolve these results is independent of quark flavors produced in the final states in Fig. 8. Since both possibilities of the photon interacting with the quark or the antiquark are already included in the sum of the two diagrams in Fig. 8, we need sum over the five light flavors only once, and not repeat the sum for antiquark flavorsf .

B.3 Projected Hadronic Tensor
In this subsection we obtain hadronic tensors by integrating the squared amplitudes obtained in App. B.1 and App. B.2 using the phase space integrations in Eqs. (3.15) and (3.18). The latter goes over the four regions in Fig. 1 with a different formula for τ (x, v) depending on which hemispheres the two final-state particles enter.

B.3.1 Quark tensor
For P µ P ν W q µν , only the real emission contribution Eq. (B.8) is nonzero, and it contains no IR divergence, so we can safely set = 0. Using Eq. (3.18) to integrate Eq. (B.8) over the four regions in Fig. 1, we obtain the contributions where the generalized theta function Θ 0 is defined in Eq. (4.6). The sum of the four contributions in Eq. (B.11) gives the result: The contribution from the real diagrams in Eq. (B.7) is more involved. We must integrate Eq. (B.7) over the two-body phase space using Eq. (3.18). We consider in turn the four contributions −g µν W real (a,b,c,d) µν corresponding to the four regions in Fig. 1.
In region (a), where x < v < 1 − x and x < 1/2, the integrand in Eq. (B.7) is finite and we can set = 0 in Eqs. (3.18) and (B.7), giving In where we converted to the MS scheme using Eq. (B.6), and the singular and nonsingular parts of the finite terms are given by Now we collect all pieces contributing to −g µν W q µν and sum them together.
Here we separately write IR divergent and finite terms in Eq. (B.20a) in order to clearly show the structure of the result, which we anticipated above in Eq. (3.7). From this result we extract the matching coefficient w q G in Eq. (3.11). The functions S q i are coefficients of singular terms in τ , R q is regular in τ , and ∆ i are coefficients of delta functions. They are given by

B.3.2 Gluon tensor
The calculation of the hadronic tensor for the gluon state follows the same steps as for the quark state. For the projection P µ P ν W g µν , we insert Eq. (B.10) into the two-body phase space integral Eq. (3.18), and obtain The integration in Eq. For the projection −g µν W g µν , we insert Eq. (B.9) into Eq. (3.18), and obtain in the four different regions in Fig. 1, where we again separately write the IR divergent and finite terms to reflect the structure anticipated in Eq. (3.7). This result gives the matching coefficient w g G in Eq. (3.11). The functions R g and ∆ g are defined by

C Separating Singular and Nonsingular Parts of Hadronic Tensor
Here, we isolate the singular and nonsingular parts of the projections of the hadronic tensor for quark and gluon initial states computed in App. B. The tensor is obtained by convolving short distance coefficients determined by perturbative matching in Sec. 3.1 with PDFs as in Eq. (3.1). The nonsingular part is obtained by subtracting singular part of the W µν tensor that has been already calculated by using SCET in [28].
One can also separate singular and nonsingular parts by isolating the structures δ(τ ) and L n (τ ) that encode the most singular terms in the τ → 0 limit in Eqs. (B.20a) and (B.25a). The nonsingular part is then obtained by subtracting these terms from Eqs. (B.20a) and (B.25a). There is no singular term in Eqs. (B.12a) and (B.22a). We can separately carry out perturbative matching for singular part and nonsingular part and determine the short distance coefficients of each part.
We write hadronic tensors in terms of three pieces associated with PDFs for q,q, g In Eq. (C.1) the factor 1/x 2 is factored out to clarify that it comes from the product of proton momenta P µ P ν .
As we promised we present the results in terms of singular and nonsingular parts The singular parts A sing i and B sing i can be extracted from the calculation of the singular cross section in [28], giving A sing q,g = 0 , (C.5a) where P qq and P qg are given in Eq. (4.7). The antiquark contributions A sinḡ q and B sinḡ q are obtained by simply replacing q →q in Eqs. (C.5a) and (C.5b). We now include the sum over flavors in both the quark and gluon contributions.
The parameters e B,J,S in Eq. (D.1) are used to perform variations of the scales µ B,J,S to estimate uncertainties from omitted higher-order corrections to beam, jet, and soft functions. By default e B,J,S = 0, and are varied away from zero according to Eq. (D.5) below. The function g(τ ) = θ(t 3 − τ ) (1 − τ /t 3 ) 2 is designed to go to zero beyond τ = t 3 , where the resummation is turned off with µ H = µ B = µ J = µ S = µ ns , and it no longer makes sense to have an individual variation of the scales µ B,J,S . This parameterization maintains the relations µ J = √ µ H µ S and µ B = √ µ H µ S for the default values (e B,J,S = 0). Theoretically the function µ run (τ ) must be chosen to satisfy several key properties to ensure the proper treatment of different regions of τ : 1. In the region ln τ α −1 s where logs of τ need to be resummed, it follows "canonical" scaling µ S ∼ Qτ and µ B,J ∼ Q √ τ .
2. For very small τ ∼ Λ QCD /Q it reaches a plateau at a constant value µ 0 where µ 0 1 GeV (above Λ QCD ). This is the nonperturbative regime where a shape function becomes necessary.
3. For larger τ ∼ 1 (where τ < 1) it becomes equal to a constant value µ independent of τ . This is the region where the resummation is turned off and the prediction reverts to fixed-order. 4. It must smoothly interpolate between each pair of regions.
Various parameters are varied to account for the residual ambiguity in satisfying these criteria. One choice that satisfies these criteria is the profile function, This is what we use in the singular part of the cross section in Eq. (5.1), with the corresponding µ S (τ ) illustrated in Fig. 5. Other choices for the profile function are also possible, see, e.g., [28]. The function µ run (τ ) in Eq. (D.2) is linear in τ with a slope r from t 1 to t 2 so that the value of µ run sets µ B,J,S to be canonical via Eq. (D.1). The function approaches µ 0 below t 1 , and µ above t 2 via a smoothly rising function ζ. The requirement of continuity for µ run (τ, µ) and its first derivative at t 1 , t 2 , and t 3 , determine the parameters α, β and constrain the function ζ(τ, t 2 , t 3 ) at t 2 and t 3 , for which we choose two connected quadratic polynomials: The default central values of the parameters that we choose are: The central values of structure function and cross section results plotted in Sec. 5 correspond to the use of these parameters. Above t 2 the resummation effect is being gradually turned off, and near t 3 the fixed order contribution dominates. We choose t 3 to be roughly the size of τ max . For t 2 we require that it well separated from t 3 by more than 0.3 for smooth turn-off of the resummation, and that it be close to the region where the nonsingular and fixed-order singular parts are of the same size. The value of t 2 determined in this way depends on x, and is well approximated for x 0.7 with a linear fit as in Eq. (D.4).
To estimate theoretical uncertainties in the cross section Eq. (5.1) due to missing higher order terms in fixed-order and resummed perturbation theory, the scales µ H , µ B,J , and µ S are varied by changing µ and e B,J,S Eq. (D.1). We also vary the points t 1 , t 2 , and t 3 and µ 0 . Each parameter is separately varied one by one while keeping the others at their default values. The variations we perform around the central values are as follows:

E Resummed Singular Cross Section
Here, we collect expressions for the resummed singular part of the cross section in Eq. (5.1) that were obtained in [28] using SCET. We provide the expressions that are necessary to obtain the resummed results in Sec. 5 at NLL accuracy. For further details on the factorization and resummation procedure see Ref. [28].
(E.6) At NLL, we only need γ i to one loop and Γ q to two loops [62], as well as the two-loop beta function β. In the MS scheme the coefficients in Eq. (E.6) used in Eq. (E.4) are given by The anomalous dimension for the soft function is obtained from the consistency relation γ S = −γ q H /2 − γ q B . In the cross section Eq. (E.1), individual factors on the right-hand side depend on the overall factorization scale µ, but in the combination of all terms, this depends cancel out completely at any fixed order in either fixed-order or resummed perturbation theory. In contrast, the dependence of Eq. (E.1) on µ H , µ B , µ J , and µ S only cancels out order-by-order in resummed perturbation theory. So at any given order there is always residual dependence on these four variables that is cancelled by higher-order terms. This residual dependence is utilized as a measure of the remaining theoretical uncertainty.
When all these scales are set to be same µ H = µ J = µ B = µ S , Eq. (E.2) reduces to zero, the resummation factors on the first line of the right-hand side of Eq. (E.1) become unity, and Eq. (E.1) reduces to the fixed-order singular part which is given in Eq. (4.8). The fixed-order parts in Eq. (E.1) are given by [δ jq C F P qq (z) + δ jg T F P qg (z)] ln z , (E.8b) where H(Q 2 , µ H ) is hard function and J n , I qq n , I qg n , S n are the coefficients of jet, beam, and soft functions and we need the function and coefficients at O(α s ). Note that the coefficients functions contain logarithms of their last argument and the hard function also depends on the logarithm ln(Q 2 /µ 2 where the coefficients G n (α s , λ) are expressed in terms of G n (α s ) in Eq. (E.10) as G −1 (α s , λ) = G −1 (α s ) + (E.18) The resummed differential distribution can be written in similar pattern to Eq. (E.1), which we do not write out explicitly here. Alternatively, the differential distribution can be obtained by numerically differentiating the cumulant in Eq. (E.1) dσ sing dτ = lim →0σ c sing (τ + ; µ i (τ )) −σ c sing (τ − ; µ i (τ )) 2 , (E. 19) which corresponds to differentiating the explicit τ dependence inσ c but not the dependence inside µ i (τ ). See footnote 4 on why we choose this procedure.