Large Transverse Momentum in Semi-Inclusive Deeply Inelastic Scattering Beyond Lowest Order

TL;DR

This paper revisits second-order calculations of large transverse momentum in semi-inclusive deep inelastic scattering, confirming previous results and highlighting discrepancies with existing parton distribution and fragmentation functions.

Contribution

It provides a detailed second-order calculation of the transverse momentum dependence, supporting the existence of tension with current parton functions.

Findings

Confirmation of previous $O(\alpha_s^2)$ calculations

Identification of tension with current parton distribution functions

Reinforcement of the need for improved theoretical models

Abstract

Motivated by recently observed tension between $O\left(\alpha_s^2\right)$ calculations of very large transverse momentum dependence in both semi-inclusive deep inelastic scattering and Drell-Yan scattering, we repeat the details of the calculation through $O\left(\alpha_s^2\right)$ transversely differential cross section. The results confirm earlier calculations, and provide further support to the observation that tension exists with current parton distribution and fragmentation functions.

Figures

Tables

Table 1. Table 1: Correspondence between real and virtual graphs at order α s ​ ( Q ) 2 subscript 𝛼 𝑠 superscript 𝑄 2 \alpha_{s}(Q)^{2} . The superscript of the hard part Γ ∈ { g , p ​ p } Γ 𝑔 𝑝 𝑝 \Gamma\in\{g,pp\} is supressed. For virtual graphs, the additional gluon is inside one of the blobs.

virtual		real		combined
${\mathscr{H}}_{gq;\overline{q}}$		${\mathscr{H}}_{gq;\overline{q}g}$		${\mathscr{H}}_{gq}$
${\mathscr{H}}_{qq;g}$		${\mathscr{H}}_{qq;gg}$ $+$ ${\mathscr{H}}_{qq;q\overline{q}}$ $+$ ${\mathscr{H}}_{qq;q^{\prime }{\overline{q}}^{\prime }}$		${\mathscr{H}}_{qq}$
${\mathscr{H}}_{qg;q}$		${\mathscr{H}}_{qg;qg}$		${\mathscr{H}}_{qg}$
N/A		${\mathscr{H}}_{gg;q\overline{q}}$		${\mathscr{H}}_{gg}$
N/A		${\mathscr{H}}_{q\overline{q};qq}$		${\mathscr{H}}_{q\overline{q}}$
N/A		${\mathscr{H}}_{q{q}^{\prime };q{\overline{q}}^{\prime }}$		${\mathscr{H}}_{q{q}^{\prime }}$

Table 2. Table 2: The first column of real and virtual are the hard parts as in Table  1 . The second column gives the hard parts needed other than those in the first column. The superscript of the hard part Γ ∈ { g , p ​ p } Γ 𝑔 𝑝 𝑝 \Gamma\in\{g,pp\} is supressed. e q 2 superscript subscript 𝑒 𝑞 2 e_{q}^{2} , e q ′ 2 superscript subscript 𝑒 superscript 𝑞 ′ 2 e_{q^{\prime}}^{2} , and e q ​ e q ′ subscript 𝑒 𝑞 subscript 𝑒 superscript 𝑞 ′ e_{q}e_{q^{\prime}} are various combinations of quark electric charges. A charge combination as a superscript denotes terms proportional to that particular charge structure in a hard part.

virtual		real
${\mathscr{H}}_{gq;\overline{q}}$	${\mathscr{H}}_{g\overline{q};q}={\mathscr{H}}_{gq;\overline{q}}$	${\mathscr{H}}_{gq;\overline{q}g}$	${\mathscr{H}}_{g\overline{q};qg}={\mathscr{H}}_{gq;\overline{q}g}$
${\mathscr{H}}_{qq;g}$	${\mathscr{H}}_{\overline{q}\overline{q};g}={\mathscr{H}}_{qq;g}$	${\mathscr{H}}_{qq;gg}$	${\mathscr{H}}_{\overline{q}\overline{q};gg}={\mathscr{H}}_{qq;gg}$
		${\mathscr{H}}_{qq;q\overline{q}}$	${\mathscr{H}}_{\overline{q}\overline{q};q\overline{q}}={\mathscr{H}}_{qq;q\overline{q}}$
		${\mathscr{H}}_{qq;q^{\prime }{\overline{q}}^{\prime }}$	${\mathscr{H}}_{\overline{q}\overline{q};q^{\prime }{\overline{q}}^{\prime }}^{e_q^2}={\mathscr{H}}_{qq;q^{\prime }{\overline{q}}^{\prime }}^{e_q^2};{\mathscr{H}}_{\overline{q}\overline{q};q^{\prime }{\overline{q}}^{\prime }}^{e_{q^{\prime }}^2}={\mathscr{H}}_{qq;q^{\prime }{\overline{q}}^{\prime }}^{e_{q^{\prime }}^2}$
${\mathscr{H}}_{qg;q}$	${\mathscr{H}}_{\overline{q}g;\overline{q}}={\mathscr{H}}_{qg;q}$	${\mathscr{H}}_{qg;qg}$	${\mathscr{H}}_{\overline{q}g;\overline{q}g}={\mathscr{H}}_{qg;qg}$
N/A		${\mathscr{H}}_{gg;q\overline{q}}$	N/A
N/A		${\mathscr{H}}_{q\overline{q};qq}$	${\mathscr{H}}_{\overline{q}q;\overline{q}\overline{q}}={\mathscr{H}}_{q\overline{q};qq}$
N/A		${\mathscr{H}}_{q{q}^{\prime };q{\overline{q}}^{\prime }}$	${\mathscr{H}}_{\overline{q}{\overline{q}}^{\prime };\overline{q}{q}^{\prime }}={\mathscr{H}}_{q{q}^{\prime };q{\overline{q}}^{\prime }}$
			${\mathscr{H}}_{\overline{q}{q}^{\prime };\overline{q}{\overline{q}}^{\prime }}^{e_q^2}={\mathscr{H}}_{q{q}^{\prime };q{\overline{q}}^{\prime }}^{e_q^2};{\mathscr{H}}_{\overline{q}{q}^{\prime };\overline{q}{\overline{q}}^{\prime }}^{e_{q^{\prime }}^2}={\mathscr{H}}_{q{q}^{\prime };q{\overline{q}}^{\prime }}^{e_{q^{\prime }}^2};{\mathscr{H}}_{\overline{q}{q}^{\prime };\overline{q}{\overline{q}}^{\prime }}^{e_q{e}_{q^{\prime }}}=-{\mathscr{H}}_{q{q}^{\prime };q{\overline{q}}^{\prime }}^{e_q{e}_{q^{\prime }}}$
			${\mathscr{H}}_{q{\overline{q}}^{\prime };q{q}^{\prime }}^{e_q^2}={\mathscr{H}}_{q{q}^{\prime };q{\overline{q}}^{\prime }}^{e_q^2};{\mathscr{H}}_{q{\overline{q}}^{\prime };q{q}^{\prime }}^{e_{q^{\prime }}^2}={\mathscr{H}}_{q{q}^{\prime };q{\overline{q}}^{\prime }}^{e_{q^{\prime }}^2};{\mathscr{H}}_{q{\overline{q}}^{\prime };q{q}^{\prime }}^{e_q{e}_{q^{\prime }}}=-{\mathscr{H}}_{q{q}^{\prime };q{\overline{q}}^{\prime }}^{e_q{e}_{q^{\prime }}}$

Equations

$$l(l)+\text{Proton}(P)\to l^{\prime }(l^{\prime })+\text{Hadron}(P_H)+X,$$

$$4P_{\mathrm{H}}^0{E}^{\prime }\frac{\mathrm{d}{\sigma }_H}{\mathrm{d}{}^3\mathbf{l}{}^{\prime }\mathrm{d}{}^3{\mathbf{P}}_{\mathrm{H}}}=\frac{2{\alpha }_{\mathrm{em}}^2}{S_{lP}{Q}^4}{L}_{\mu \nu }{W}^{\mu \nu },$$

$$\frac{\mathrm{d}{\sigma }_H}{\mathrm{dx}\mathrm{dy}\mathrm{dz}\mathrm{d}{P}_{\text{H,T}}^2}=\frac{{\pi }^2{\alpha }_{\mathrm{em}}^{2}y}{2Q^{4}z}{L}_{\mu \nu }{W}^{\mu \nu },$$

$$L_{\mu \nu }\equiv 2(l_{\mu }{l}_{\nu }^{\prime }+l_{\mu }^{\prime }{l}_{\nu }-g_{\mu \nu }l\cdot l^{\prime }),$$

$$W^{\mu \nu }(P,q,P_H)\equiv \frac{1}{(2\pi {)}^4}\sum _X\int \mathrm{d}{}^{4}z{e}^{iq\cdot z}⟨P,S|j^{\mu }(z)|P_H,X⟩⟨P_H,X|j^{\nu }(0)|P,S⟩,$$

$$W^{\mu \nu }=\left(-g^{\mu \nu }+\frac{q^{\mu }{q}^{\nu }}{q^2}\right)F_1+\frac{(P^{\mu }-q^{\mu }P\cdot q/q^2)(P^{\nu }-q^{\nu }P\cdot q/q^2)}{P\cdot q}{F}_2.$$

$$\frac{\mathop{\mathrm{d}\mskip-0.59999mu\sigma_{H}}\nolimits}{\mathop{\mathrm{d}\mskip-0.59999mux}\nolimits\mathop{\mathrm{d}\mskip-0.59999muy}\nolimits\mathop{\mathrm{d}\mskip-0.59999muz}\nolimits\mathop{\mathrm{d}\mskip-0.59999muP_{\text{H,T}}^{2}}\nolimits}=\frac{\pi^{2}\alpha_{\rm em}^{2}}{zxyQ^{2}}\Biggl{[}xy^{2}F_{1}+(1-y)F_{2}\Biggr{]}\,.$$

$${\mathrm{P}}_g^{\mu \nu }=g^{\mu \nu },{\mathrm{P}}_{PP}^{\mu \nu }=P^{\mu }{P}^{\nu }.$$

$$F_1={\mathrm{P}}_1^{\mu \nu }{W}_{\mu \nu },F_2={\mathrm{P}}_2^{\mu \nu }{W}_{\mu \nu }$$

$${\mathrm{P}}_1^{\mu \nu }=\frac{1}{1-ϵ}\left(-\frac{1}{2}{\mathrm{P}}_g^{\mu \nu }+\frac{2x^2}{Q^2}{\mathrm{P}}_{PP}^{\mu \nu }\right),{\mathrm{P}}_2^{\mu \nu }=\left(\frac{3-2ϵ}{1-ϵ}\right)\frac{4x^3}{Q^2}{\mathrm{P}}_{PP}^{\mu \nu }-\frac{x}{1-ϵ}{\mathrm{P}}_g^{\mu \nu }.$$

$${\mathbf{q}}_{\text{T}}=-\frac{{\mathbf{P}}_{\text{H,T}}}{z}.$$

$$4P_{\mathrm{H}}^0{E}^{\prime }\frac{\mathrm{d}{\sigma }_H}{\mathrm{d}{}^{3-2ϵ}\mathbf{l}{}^{\prime }\mathrm{d}{}^{3-2ϵ}{\mathbf{P}}_{\mathrm{H}}}={\int }_x^1\frac{\mathrm{d\xi }}{\xi }{\int }_z^1\frac{\mathrm{d\zeta }}{{\zeta }^2}\left(4k_1^0{E}^{\prime }\frac{\mathrm{d}{\hat{\sigma }}_{ij}}{\mathrm{d}{}^{3-2ϵ}\mathbf{l}{}^{\prime }\mathrm{d}{}^{3-2ϵ}{\mathbf{k}}_1}\right)f_{i/P}(\xi )d_{H/j}(\zeta ).$$

$$\hat{x}\equiv \frac{Q^2}{(2p\cdot q)}=\frac{x}{\xi },\hat{z}\equiv \frac{k_1\cdot p}{(p\cdot q)}=\frac{z}{\zeta },k_{1\text{T}}\equiv \frac{P_{H,\text{T}}}{\zeta }.$$

$$4k_1^0{E}^{\prime }\frac{\mathrm{d}{\hat{\sigma }}_{ij}}{\mathrm{d}{}^{3-2ϵ}\mathbf{l}{}^{\prime }\mathrm{d}{}^{3-2ϵ}{\mathbf{k}}_1}=\frac{2{\alpha }_{\mathrm{em}}^2}{\hat{s}{Q}^4}{L}_{\mu \nu }{\hat{W}}^{\mu \nu },$$

$${\hat{W}}_{\mu \nu ,ij}\equiv \frac{1}{2}\frac{1}{(2\pi {)}^4}\sum _{s,X}\int \mathrm{d}{}^{4}w{e}^{iq\cdot w}⟨p_i,s|j_{\mu }(w)|k_{1j},X⟩⟨k_{1j},X|j_{\nu }(0)|p_i,s⟩.$$

$$W^{\mu \nu }(P,q,P_H)={\int }_x^1\frac{\mathrm{d\xi }}{\xi }{\int }_z^1\frac{\mathrm{d\zeta }}{{\zeta }^2}{\hat{W}}_{ij}^{\mu \nu }(q,x/\xi ,z/\zeta )f_{i/P}(\xi )d_{H/j}(\zeta ).$$

$${\hat{W}}_{ij}^{\mu \nu }=\left(-g^{\mu \nu }+\frac{q^{\mu }{q}^{\nu }}{q^2}\right){\hat{F}}_{1,ij}+\frac{(p^{\mu }-q^{\mu }p\cdot q/q^2)(p^{\nu }-q^{\nu }p\cdot q/q^2)}{p\cdot q}{\hat{F}}_{2,ij}.$$

$$F_1(x,z,Q^2,{\mathbf{q}}_{\text{T}})$$

$$F_2(x,z,Q^2,{\mathbf{q}}_{\text{T}})$$

$$\{{\mathrm{P}}_g^{\mu \nu }{\hat{W}}_{\mu \nu }^{(N)};{\mathrm{P}}_{PP}^{\mu \nu }{\hat{W}}_{\mu \nu }^{(N)}\}$$

$$\equiv \{{\mathrm{P}}_g^{\mu \nu }{\hat{W}}_{\mu \nu }^{(N)};{\mathrm{P}}_{PP}^{\mu \nu }{\hat{W}}_{\mu \nu }^{(N)}{\}}_{\mathrm{unsub}}-\mathrm{Subtractions}.$$

$${\gamma }^{\ast }(q)+\text{parton}(p)\to \text{parton}(k_1)+\text{(N−1) spectator partons},$$

$${\mathrm{P}}_{\Gamma }^{\mu \nu }{W}_{\mu \nu }(P,q,P_H)=\sum _N{\int }_x^1\frac{\mathrm{d\xi }}{\xi }{\int }_z^1\frac{\mathrm{d\zeta }}{{\zeta }^2}\{{\mathrm{P}}_g^{\mu \nu }{\hat{W}}_{\mu \nu }^{(N)};{\mathrm{P}}_{PP}^{\mu \nu }{\hat{W}}_{\mu \nu }^{(N)}\}{f}_{i/P}(\xi )d_{H/j}(\zeta ),$$

$$s=$$

$$t_i=$$

$$u_i=$$

$$s_{ij}=$$

$$s=$$

$$t=$$

$$u=$$

$$\mathrm{dt}\mathrm{du}=$$

$$\int \mathrm{d\xi }\mathrm{d\zeta }\cdots ={\int }_A^1\mathrm{d\xi }{\int }_0^B\mathrm{d}{s}_{23}\frac{{\hat{x}}^2{P}_{H,\text{T}}^2+\hat{x}{z}^2{Q}^2(1-\hat{x})}{z(Q^2(1-\hat{x})-s_{23}\hat{x}{)}^2}\cdots ,$$

$$\zeta =\frac{\hat{x}{P}_{H,\text{T}}^2+z^2{Q}^2(1-\hat{x})}{z(Q^2(1-\hat{x})-s_{23}\hat{x})},$$

$$A=x+\frac{x{P}_{H,\text{T}}^2}{z(1-z)Q^2},$$

$$B=Q^2\left(\frac{1}{\hat{x}}-1\right)(1-z)-\frac{P_{H,\text{T}}^2}{z}.$$

Full text

††preprint: JLAB-THY-19-2897

I Introduction

In a previous article Gonzalez-Hernandez et al. (2018a), we discussed the semi-inclusive deep inelastic scattering (SIDIS) process:

[TABLE]

and we highlighted the challenge of finding agreement between $O\left(\alpha_{s}^{2}\right)$ calculations and existing SIDIS data in the very large transverse momentum transverse momentum limit where standard collinear factorization is expected to be valid. One motivation is that obtaining a description of the small transverse momentum behavior associated with nucleon structure requires a good understanding of the matching to large transverse momentum where a transverse momentum dependent (TMD) factorization description fails. Given the current focus on using deeply inelastic hadro-production to access nucleon structure sensitivity sensitivity, it is imperative to examine theoretical framework for the full $q_{\text{T}}$-range in greater detail. In the language of Ref. Gonzalez-Hernandez et al. (2018a), we are interested in this paper in what was there called “region 3” behavior, corresponding to where $q_{\text{T}}$ is so large that the small $q_{\text{T}}/Q$ approximations associated with TMDs are not reliable, but where ordinary collinear factorization should be applicable and reliable. Calculations to $O\left(\alpha_{s}^{2}\right)$ have existed for some time Daleo et al. (2005); Kniehl et al. (2005). The main observation of Ref. Gonzalez-Hernandez et al. (2018a) was that, while the $O\left(\alpha_{s}^{2}\right)$ correction gives an order of magnitude increase over the leading order, that still is not sufficient to achieve reasonable agreement with data for $Q$ in the region of one to several GeVs, transverse momentum of order $Q$, and for moderate Bjorken-$x$. Combined with similar observations concerning the Drell-Yan pointed out in Bacchetta et al. (2019), this points to general tension between transverse momentum dependent cross sections at $q_{\text{T}}\sim Q$ and collinear factorization.

There are a number of potential explanations or solutions, including a direct re-tuning of collinear parton distribution and/or fragmentation functions to transversely differential SIDIS cross sections. (A proposal to constrain gluon PDFs in transversely differential Drell-Yan cross sections was made already in Berger et al. (1998).) But before proceeding to consider these directions, it is important to validate the $O(\alpha_{S}^{2})$ in Daleo et al. (2005); Kniehl et al. (2005) that led to our conclusions in Gonzalez-Hernandez et al. (2018a). Therefore, we have in this paper repeated the large transverse momentum order $O\left(\alpha_{s}^{2}\right)$ calculation following a slightly different formal framework. We are able to reproduce the results in Daleo et al. (2005) very closely, thus bolstering the observations we made earlier in Gonzalez-Hernandez et al. (2018a).

The specific purposes of this article are as follows: i.) to lay out the logical steps of the calculation with enough detail, we hope, to lead to clues as to how to improve phenomenological agreement, ii.) to make available a convenient numerical implementation for the large $q_{\text{T}}$ region and, iii.) to present some quantitative results relevant to current experimental programs such as those at COMPASS and Jefferson Lab 12 GeV. Overall, our results support the general observations made in Gonzalez-Hernandez et al. (2018a).

In our calculations we use qgraf Nogueira (1993) to generate Feynman graphs and FORM Kuipers et al. (2013) to carry out spinor and color traces for the amplitudes. The renormalization counter terms are computed in MATHEMATICA with the packages FeynArts Kublbeck et al. (1990); Hahn (2001) and FeynCalc Mertig et al. (1991); Shtabovenko et al. (2016). Further analytic manipulations were performed in MATHEMATICA. The code to compute the cross sections are publically available at Gonzalez-Hernandez et al. (2018b).

In Sec. II, we explain our setup and notation. In Sec. III we summarize the organization of our calculations. This includes a classification of the partonic Feynman graphs needed at order $\alpha_{s}(Q)^{2}$ and a discussion of the phase space integrals for multiparton final states. We discuss the results of the calculation in Sec. IV, and give concluding remarks in Sec. V.

II Notation and Conventions

We will express quantities in terms of the conventional kinematical variable $z\equiv P_{H}\cdot P/(P\cdot q)$. $\bm{P}_{H,\text{T}}$ is the Breit frame transverse momentum of the produced hadron, and $P$ and $q$ are the four-momenta of the incoming target hadron and the virtual photon respectively. We will focus on the unpolarized and azimuthally-independent cross section since this is the most straight forward observable to calculate. (Although the results have potential implications also for polarization dependent observables.) Many general treatments of SIDIS as a process are available Meng et al. (1992); Levelt and Mulders (1994); Meng et al. (1996); Mulders and Tangerman (1996); Nadolsky et al. (1999, 2001); Barone et al. (2002); Ji et al. (2004); Bacchetta et al. (2004); Koike et al. (2006); Bacchetta et al. (2007, 2008). Our notation is consistent with trends among these, with modifications as needed for our current purposes.

The unpolarized differential cross section is

[TABLE]

or

[TABLE]

where $\sigma_{H}$ is the unpolarized hadronic cross section. $S_{lP}=(l+P)^{2}$. The hadron transverse momentum $\bm{P}_{H,\text{T}}$ is defined in a frame where the photon and incoming hadron are back-to-back (a “photon” frame). The Bjorken $x$ and the $y$ variable are the usual definitions $x\equiv Q^{2}/2P\cdot q$ and $y\equiv P\cdot q/P\cdot l$. $L_{\mu\nu}$ is the usual leptonic tensor

[TABLE]

and $W^{\mu\nu}$ is the hadronic tensor for SIDIS:

[TABLE]

where we have omitted spin and azimuthal angle dependent terms since we do not consider these in this paper. Also, we assume that $x$ and $1/Q$ are small enough that both the proton and lepton mass can be dropped in kinematical and phase space factors. The normalization convention in Eq. (4) is so that the prefactor on the right hand side of Eq. (2) is similar to the unpolarized case.

The unpolarized structure functions $F_{1}$ and $F_{2}$ are defined by the usual gauge invariant decomposition

[TABLE]

Then, the cross section is

[TABLE]

In calculations, it is convenient to work with Lorentz invariant structure function extraction tensors ${\rm P}_{\Gamma}^{\mu\nu}$ with $\Gamma\in\{g,PP\}$ where

[TABLE]

Then,

[TABLE]

where in $4-2\epsilon$ dimensions,

[TABLE]

(We we always use the massless target approximation – see discussion in Moffat et al. (2019).) It is useful to express transverse momentum in terms of

[TABLE]

In a frame where the incoming and outgoing hadrons are back-to-back, $\bm{q}_{\text{T}}$ is the transverse momentum of the virtual photon (assuming, as always in this paper, that external hadron masses are negligible).

Our treatment of factorization will follow the general style of Collins (2011). The factorization theorem that relates the hadronic and partonic differential cross sections in SIDIS is

[TABLE]

The $1/\xi$ is from the partonic flux factor, and the 1$/\zeta^{2}$ is from the conversion between ${\bf k}_{1}$ and ${\bf P}_{\rm H}$. The indices $i$ and $j$ denote, respectively, the flavors of the parton in the proton (with a momentum fraction $\xi$) and of the outgoing parton that fragments into hadron $H$, whose momentum is a fraction $\zeta$ of parton $j$ momentum. The incoming and outgoing parton momenta $p$ and $k_{1}$ satisfy $p=\xi P$ and $k_{1}=P_{H}/\zeta$. (Indices $i$ and $j$ for incoming and outgoing partons $p_{i}$ and $k_{1,j}$ are not shown explicitly but are understood). $f_{i/P}(\xi)$ and $d_{H/j}(\zeta)$ are the collinear parton distribution and fragmentation functions respectively. It is also useful to define partonic variables

[TABLE]

The differential partonic hard part in Eq. (11) is finite and well-behaved, and at large transverse momentum it starts at $O\left(\alpha_{s}\right)$.

The unpolarized partonic structure tensor for scattering off parton $i$ into parton $j$ is defined in exact analogy with the hadronic tensor:

[TABLE]

with,

[TABLE]

Thus, from Eq. (11),

[TABLE]

The partonic structure function decomposition is

[TABLE]

Then,

[TABLE]

III Organization

III.1 Basic Setup

For a process with $N$ final state partons, ${\rm P}_{g}^{\mu\nu}W_{\mu\nu}$ and ${\rm P}_{PP}^{\mu\nu}W_{\mu\nu}$ are squared amplitudes integrated over the $N$-particle phase space for the outgoing partons,

[TABLE]

Here $\lvert M_{g}\rvert^{2}$ or $\lvert M_{PP}\rvert^{2}$ is the squared amplitude for the process

[TABLE]

with the polarization sum of the virtual photon replaced by $g^{\mu\nu}$ or $P^{\mu}P^{\nu}$. The subtraction terms in Eq. (19) are needed to remove double counting with lower orders of perturbation theory and cancel singularities in the first term. The form of the subtraction term will be discussed in Sec. III.2. It is assumed in Eq. (19) that all integrals that allow kinematical $\delta$-functions to be evaluated have been performed. Also, the phase space factors associated with $k_{1}$ are excluded from the partonic phase space since they give the $\hat{z}$ and $k^{2}_{1\text{T}}$ dependence of the differential partonic cross section. The $1/(2\pi)^{4}$ on the right hand side of Eq. (19) is the same factor in Eq. (4). The $\mathop{\mathrm{d}\mskip-0.59999mu\Pi}\nolimits^{(N)}$ represents a generic phase space factor for $2\to N$ scattering (see Eq. (35) and Eq. (39)).

Thus it is convenient to express the structure functions in the form

[TABLE]

where the $\sum_{N}$ is a sum over all possible partonic final states.

We will express the phase space in terms of Mandelstam variables:

[TABLE]

where $i,j=1,2,3$, with the labeling in Fig. 1. For simplicity, $u_{1}$ and $t_{1}$ will be abbreviated as $u$ and $t$ from here on. All Mandelstam variables will refer to the partonic cross sections. Occasionally it will be useful to change kinematical variables, For example,

[TABLE]

The $O(\alpha_{s}^{0})$ contribution to Eq. (19) is kinematically constrained to $k_{1\text{T}}=0$. At $O(\alpha_{s}^{1})$, only tree level processes contribute, and no singularities appear. At $O(\alpha_{s}^{2})$, soft, collinear and UV singularities arise as $1/\epsilon$ and $1/\epsilon^{2}$ poles in dimensional regularization with space-time dimension $n\equiv 4-2\epsilon$.

These singularities cancel in the sum of real, virtual and counterterm graphs and after applying collinear factorization. More discussion of the singularity structure at $O(\alpha_{S}^{2})$ is available in Daleo et al. (2005); Arnold and Reno (1989). The $O(\alpha_{s})$ partonic cross section is described in details in Nadolsky et al. (1999), and our results match that calculation. At $O(\alpha_{s}^{2})$, details below can also be found in, e.g.,  Aversa et al. (1989); Gordon and Vogelsang (1993); Ellis et al. (1980, 1983); Arnold and Reno (1989). The change of variables in Eq. (20) from momentum fractions to Mandelstam variables is

[TABLE]

where $\zeta$ is replaced by $s_{23}$, the virtuality of the spectator parton system. ($s_{23}$ is $k_{X}^{2}$ in the notation of Gonzalez-Hernandez et al. (2018a)). Overall momentum conservation $s+u+t=-Q^{2}+s_{23}$ ($s\equiv(p+q)^{2}=Q^{2}(1-\hat{x})/\hat{x}$) gives

[TABLE]

The graphical structures needed at order $\alpha_{s}(Q)^{2}$ can be classified as in Fig. 2. Each diagram corresponds to a contribution to the squared amplitude, and the blobs include all possible attachments. For the $2\to 2$ subprocesses, the blobs include up to one loop corrections in one side of the cut to give $O(\alpha_{S}^{2})$ accuracy to the hard parts. The black dot marks the fragmenting parton and the other lines are integrated over all phase space.

The following are assumed and not shown explicitly in the graphs:

• •

Subtractions, consistent with factorization, for all collinear divergences.

• •

One QCD loop at all possible positions inside the blob in virtual processes.

• •

Associated UV counter term graphs for virtual processes.

Define the hard parts $\mathcal{H}$ for individual graphs as

[TABLE]

where $\Gamma\in\{g,pp\}$ and $i,j,k,l\in\{q,\bar{q},g\}$. As before, $i$ and $j$ label incoming and fragmenting parton flavors respectively while $k$ and $l$ label the unobserved parton flavors. The graph is virtual when only one flavor index appears after the “;”. Figures 2(V-A) through 2(V-F) correspond to the $O\left(\alpha_{s}\right)$ graphs if all virtual loops in the blobs are removed.

So, for example, $\mathcal{H}^{g}_{qq;gg}$ is represented by Fig. 2(R-A), when the partonic tensor is contracted with $g^{\mu\nu}$. Fig. 2(R-E) includes both $\mathcal{H}^{g}_{qq;q^{\prime}\bar{q}^{\prime}}$ and $\mathcal{H}^{g}_{qq^{\prime};q\bar{q}^{\prime}}$ contributions where the prime indicates a different flavor. Note that graphs like Fig. 2(R-G) give contributions like $\mathcal{H}^{g}_{q\bar{q};qq}$ or $\mathcal{H}^{g}_{q\bar{q}^{\prime};q^{\prime}q}$.

We work in the approximation that all particles are massless, except for the photon which is highly virtual, $q^{2}=-Q^{2}$. In Eq. (19), for $2\to 2$ scattering,

[TABLE]

Thus the $2\to 2$ phase space factor is

[TABLE]

For the $2\to 3$ case, we follow the strategy in Ellis et al. (1980); Gordon and Vogelsang (1993) where we work in the rest frame of the $k_{2}+k_{3}$ system, in which

[TABLE]

where $\bm{\hat{k}}_{n-3}$ denotes unit vectors for the first $n-3$ components of the $n-1$ dimensional unit spatial vector in spherical coordinates. See Fig. 3.

In the center-of-mass of $p$ and $q$, the spatial vectors $\bm{p}$, $\bm{q}$, $\bm{k_{1}}$, and $\bm{k_{2}}+\bm{k_{3}}$ are in the same plane. We then boost to the $k_{2}+k_{3}$ rest frame $\bm{p}$, where $\bm{q}$, and $\bm{k_{1}}$ are still in the same plane – see Fig. 3. We choose the spatial orientations of this coordinate system so that $\bm{p}$, $\bm{q}$, and $\bm{k_{1}}$ are in the plane created by the last two spatial axes. For these vectors, the first $n-3$ spatial components are zero. Then the scattering amplitudes do not depend on the first $n-3$ spatial components of $k_{2}$ or $k_{3}$, and the 3-body phase space simplifies (see Appendix A for useful identities relating to the 3-body final state phase space). In this frame,

[TABLE]

where in the third equality we use that the scattering amplitudes are independent of the first $n-3$ spatial components of $k_{2}$. That is,

[TABLE]

From overall momentum conservation

[TABLE]

Since $s_{23}$ can become zero in certain regions of $\hat{x}$ and $\hat{z}$, dimensional regularization is needed for the $s_{23}\to 0$ behavior. Further details of the $2\to 3$ phase space are provided in the next subsection.

To summarize the basic steps, integrals for the unobserved final state phase space are first performed over angular variables. Then, integrals over momentum fractions are expressed via Eq. (29) as integrals over $\xi$ and $s_{23}$. In the $s_{23}\to 0$ limit there are soft and collinear singularities, proportional to $\delta(s_{23})/\epsilon$ or $\delta(s_{23})/\epsilon^{2}$, that cancel after combining $2\to 3$ and the corresponding virtual $2\to 2$ processes. There are additional collinear singularities in $2\to 3$ and virtual $2\to 2$ processes which are ultimately canceled by subtraction terms in the factorization algorithm. (These contain also $\delta(s_{23})$ and plus distributions.) All contributions are integrated over the physical range of $s_{23}$ in Eqs. (29)–(32).

III.2 Factorization Order-By-Order

Once all real and virtual Feynman graphs have been computed in dimensional regularization, they must be combined in an algorithm consistent with a factorization theorem.

According to the collinear factorization theorem, hard factors are independent of the species of target and final state particles, so hard scattering calculations can be performed for free and massless partons as the target and final state without loss of generality. (In other words, we will consider perturbatively calculable parton-in-parton PDFs and and parton-to-parton FFs.) This simplifies calculations, so we start our calculations by rewriting Eq. (15) as a factorization theorem for a partonic initial and final state:

[TABLE]

All manipulations will be done in $n=4-2\epsilon$ dimensions until the very end. The left side is now an argument of $p^{\prime}$, $q$, and $k^{\prime}_{1}$, and the PDFs and FFs have subscripts $i/i^{\prime}$ and $j^{\prime}/j$. The momenta $p^{\prime}$ and $k^{\prime}_{1}$ play the roles that $P$ and $P_{\text{H,T}}$ played earlier. The initial $p^{\prime}$ has a flavor $i^{\prime}$ and the final $k^{\prime}_{1}$ has a flavor $j^{\prime}$. The second line in Eq. (41) defines the “$\otimes$” notation as the usual shorthand for the convolution integrals in $\xi$ and $\zeta$. On the last line, $W^{(\text{LO})}_{\mu\nu,i^{\prime}j^{\prime}}$ and $W^{(\text{NLO})}_{\mu\nu,i^{\prime}j^{\prime}}$ are hadronic tensors calculated to LO and NLO in the hard parts, i.e. what one normally means when one speaks of a “leading order” or “next-to-leading order” calculation in collinear pQCD.

Expanding Eq. (41) to the relevant order in $\alpha_{s}$ gives

[TABLE]

The aim is to get the hard parts, $\hat{W}^{(\text{LO})}_{\mu\nu,ij}(q,x/\xi,z/\zeta)$ and $\hat{W}^{(\text{NLO})}_{\mu\nu,ij}(q,x/\xi,z/\zeta)$. The leading order is just

[TABLE]

where the expression for $\hat{W}^{(\text{LO})}_{\mu\nu,ij}(q,x/\xi,z/\zeta)$ is simple and well-known – it is just the sum of tree level $2\to 2$ graphs contributing to $\hat{W}_{\mu\nu,i^{\prime}j^{\prime}}(q,x/\xi,z/\zeta)$.

To get $\hat{W}^{(\text{NLO})}_{\mu\nu,ij}(q,x/\xi,z/\zeta)$, first write

[TABLE]

The term in braces is the correction to the leading order so it vanishes by construction at order $\alpha_{s}$ in the hard part. Also, it contains subtractions for the overlap with the LO, so it is infrared safe through order $\alpha_{s}^{2}$ in hard scattering. Thus, we may rewrite it as:

[TABLE]

The term in braces on the second line is equal to the first line evaluated with hatted partonic variables and expanded only to order $\alpha_{s}^{2}$ in the coupling. Comparing with Eq. (42) and Eq. (44) shows that the factor in braces in Eq. (45) is $\hat{W}^{(\text{NLO})}_{\mu\nu,ij}(q,x/\xi,z/\zeta)$. We use Eq. (43) to evaluate it explicitly in terms of low order PDFs and FFs:

[TABLE]

Here, the $W^{(\text{NLO})}_{\mu\nu}(q,\hat{x},\hat{z})_{\rm unsub}$ is the sum of order $O(\alpha_{s}^{2})$ graphs contributing to $W_{\mu\nu}(q,\hat{x},\hat{z})$ without any subtractions (i.e., graphs in Fig. 2(R-A) through Fig. 2(V-F) integrated over phase space, but absent the subtractions). The second two terms are obtained from Eq. (43), expanded to order $\alpha_{s}^{2}$. The last line in Eq. (46) uses the known results for PDFs and FFs calculated to $O\left(\alpha_{s}\right)$ for massless quarks and gluons:

[TABLE]

Note that, consistent with the $\overline{\text{MS}}$ scheme, the renormalized order $\alpha_{s}$ PDFs and FFs in dimensional regularization for massless pQCD are just $S_{\epsilon}/\epsilon$ poles. There are no $\mu^{\epsilon}$ factors or logarithms of $\mu$. Thus the last line of Eq. (46) involves no factors of $\mu^{\epsilon}$. For completeness, the one-loop splitting functions Altarelli and Parisi (1977) are

[TABLE]

Our calculations need Eq. (46), contracted with specific ${\rm P}_{\Gamma}^{\mu\nu}$ tensors. The substraction scheme has been carry out without contracting the Lorentz indices of the partonic tensor. However the method is equally applicable for any kind of contraction with external momenta. In our calculations we carry out the substractions separately for the two extraction tensors in Eq.(7) and verified analytically the cancellation of all infrared and collinear singularities.

IV Results

IV.1 Combining real and virtual contributions

After the hard real and virtual contributions are calculated, they must be combined into infrared safe squared amplitudes. Table 1 shows graphs from corresponding real and virtual processes. Note that the correspondence between real and virtual processes is not one-to-one. In particular, infrared (soft and collinear) singularities in $\mathcal{H}_{qq;g}$ are canceled by three real processes. A subtlety arises when a quark loop appears on the gluon leg of $\mathcal{H}_{qq;g}$. This creates a collinear pole term proportional to $N_{f}$, the number of massless quark flavors in the loop. This pole term is canceled after adding the real processes $\mathcal{H}_{qq;q\bar{q}}$ and $\mathcal{H}_{qq;q^{\prime}\bar{q}^{\prime}}$ with all massless flavors $q^{\prime}$ (other than $q$) Gordon and Vogelsang (1993). For processes with a non-spectator gluon leg, the $N_{f}$ dependence of the collinear pole is removed by factorization. Also notice that some real processes have no corresponding virtual ones, and in these cases factorization subtractions are sufficient to remove all infrared poles. Since many graphs need to be combined, the results are presented as the six scattering hard parts listed in the last column of Table 1. We also need the hard parts for processes with reversed quark number flow in one or two open fermion lines in the graph. That is, in one or two open quark lines, quarks and their corresponding anti-quarks are interchanged. These are easily related to the hard parts already listed in Tabel 1 by using the fact that QED vertices acquire a minus sign under charge conjugation. The results are summerized in Tabel 2. Note that when a quark line links an outgoing spectator quark anti-quark pair, there is no need for the reversed quark flow. This is because the spectator momenta are integrated over and interchanging the quark and anti-quark will double count the contribution.

IV.2 Comparison with existing results

After combining the graphs in each of the six scattering channels in Table 1, we have verified explicitly that all single and double poles cancel. The terms left are finite in the limit $\epsilon\rightarrow 0$, and constitute the infrared safe hard parts in the last column of Table 1. As a check, we have compared with a (privately obtained) computer calculation that appears to reproduce the results of Ref. Daleo et al. (2005) but is modified to be consistent with the kinematics of the current experimental data. For most kinematics, the calculations agree within experimental uncertainties, but we found several possible differences.

• •

In the last row of Table 1, our calculation of terms from the ${\rm P}_{g}^{\mu\nu}$ projection with the charge structure $e_{q}^{2}$ and $e_{q}e_{q^{\prime}}$ differs from the corresponding terms calculated in the previously existing code by a minus sign. (By contrast, terms with $e_{q^{\prime}}^{2}$ agree.)

• •

As is discussed in subsection IV.1, adding the real processes $\mathcal{H}_{qq;q\bar{q}}$ and $\mathcal{H}_{qq;q^{\prime}\bar{q}^{\prime}}$ with all massless flavors $q^{\prime}$ gives terms proportional to $N_{f}$. The pole parts of these terms are canceled by the corresponding poles from $\mathcal{H}_{qq;g}$. However, a finite part remains since the real processes above with $N_{f}$ final state quark pairs give identical contribution even when the quark pairs are not collinear. These terms contribute to the third row of Table 1. In the previously existing code, we did not find an explicit dependence on $N_{f}$ in this channel.

The numerical result of our calculation can in some cases differ from the previously existing code by as much as $\sim 100\%$ for certain individual subchannels. In the kinematics of the COMPASS experiment, this translates into a discrepancy of up to $\sim 20\%$ for the overall cross section.

IV.3 Phenomenological Results

We examine the impact of the $O(\alpha_{S}^{2})$ corrections in SIDIS by plotting the NLO to LO ratio $K$-factor for the $F_{1}$ structure function. We will consider values of $z$ between $0.2$ and $0.8$ since this is a region where factorization theorems based on current fragmentation in SIDIS is conventionally expected to apply.

For $x$, $Q$, $q_{\rm T}$, we choose values of $Q=2,20{\rm~{}GeV}$, $q_{\rm T}=Q,2Q$ and $0.001\leq x\leq 1$, which correspond to kinematics ranges accessible by existing experimental facilities such as COMPASS and JLab 12. Fig. 4 shows the $K$-factor across the aforementioned kinematical range, and it shows a clear concave up shape in its dependence on $x$ with values that decrease as $Q$ varies from $Q=2$ to $20$ GeV.

It is notable that even for $Q$ as large as $2$ GeV, the $K$-factor is greater than $2$, even at its minimum value. Note also that this minimum is approximately in the region of $0.01\leq x\leq 0.1$, close to the valence region relevant to many hadron structure studies. The $K$-factor increases at smaller values of $x$, indicating that the NLO corrections becomes increasingly important for describing regions with large $(P+q)^{2}$. It also increases at large $x$ until it reaches the kinematic boundary where the phase space for the hadron production vanishes.

The enhancements at large $x$ can be traced to the logarithms of $B$ (see Eq. 32 and appendix D) which are associated with soft gluon effects near the kinematical threshold. Thus it is likely that threshold resummation becomes important in this region. Point-by-point in $q_{\text{T}}$, effects near the kinematical boundaries are more apparent than in $q_{\text{T}}$-inclusive cross sections. It is interesting to compare with observations found in Anderle et al. (2013) where threshold resummation effects for $q_{\rm T}$ integrated SIDIS were found to give sizable corrections, but not the factor of $\gtrsim 5$ enhancement needed to explain the largest $q_{\text{T}}$ present in Fig. 4.

The minimum in the $K$-factor is helpful for identifying regions in $x$ where ordinary fixed order treatments are most likely to be sufficient. Given a set of values for $z$, $q_{\text{T}}$, and $Q$, let $x_{0}$ denote the $x$ value corresponding to the minimum $K$-factor. Then bins in $x$, $Q$, $z$ ,$q_{\text{T}}$ can be classified according to $x<x_{0}$ and $x>x_{0}$. We show this in Fig. 5 for kinematical bins corresponding to recent COMPASS $q_{\text{T}}$-dependent SIDIS multiplicities $h^{\pm}$ Aghasyan et al. (2018) for $q_{\text{T}}/Q>1$. The values of $x$ the for which the $K$-factor reaches its minimum for $q_{\text{T}}\approx Q$ can be read from the plots. Away from these regions, additional resummation techniques may be needed.

It may also be necessary to update collinear PDFs outside these regions in order to fully describe the large $q_{\text{T}}$ behavior. As mentioned in Gonzalez-Hernandez et al. (2018a), the large $q_{\text{T}}$ tails of SIDIS are sensitive to the large $\xi$ PDFs and large $\zeta$ FFs. These could potentially open up new opportunities to constrain collinear PDFs at large momentum fractions.

At the same time, factorization theorems for the full $q_{\text{T}}$-dependent SIDIS spectrum, including $q_{\text{T}}\approx 0$ with non-perturbative TMD PDFs and TMD FFs, need the large $q_{\text{T}}\sim Q$ component in order to have a completely reliable tests of the factorization formalism. Ideally, such validation will take place when $q_{\text{T}}$-dependent SIDIS data are included in the simultaneous extraction of collinear PDFs and FFs as well as non-perturbative TMD PDFs and FFs in QCD global analysis.

V Conclusion

The computational tools necessary to reproduce Fig. 4, Fig. 5 and other similar calculations are available at Gonzalez-Hernandez et al. (2018b), along with documentation. As indicated in the introduction, these generally confirm earlier calculations (e.g., Daleo et al. (2005); Kniehl et al. (2005)) of the overall cross section to within about $20\%$, well within present experimental uncertainties, and thus strengthens our earlier position Gonzalez-Hernandez et al. (2018a) that significant tension exists between existing SIDIS data and collinear factorization calculations at large $q_{\text{T}}\gtrsim Q$.

We have identified region where collinear factorization appears most reliable as the regions with minimal $K$-factors, with current sets of collinear PDFs and FFs. As discussed in Sec. IV.3, minimal values for the $K$-factor lie in the region $0.01<x<0.1$, that is, approximately in valence kinematics. The increase of the $K$-factor at both smaller and larger values of $x$, may point at the importance of resummation effects outside the valence region, or the need to refit collinear PDFs and FFs. This last interpretation is consistent with that of our previous study Gonzalez-Hernandez et al. (2018a). It is also important to note that the $2\to 2$ kinematics of the $O\left(\alpha_{s}\right)$ contribution places severe kinematical constraints on the relationship between the initial and final state partons. Therefore, it is likely that the generally large $K$-factors are at least partly due simply to a kinematical suppression of the $O\left(\alpha_{s}\right)$ contribution, and are not a fundamental problem with the convergence of the perturbation series. We note that Reference Balitsky and Tarasov (2018) has addressed somewhat similar issues, but in the region where $q_{\text{T}}$ is still small enough that a $q_{\text{T}}/Q$ power expansion is still meaningful, and also in the limit of $Q^{2}\gg s$.

As a next step, we plan to refit collinear functions in the large $q_{\text{T}}$ region, using SIDIS and $e^{+}e^{-}$-annihilation to back-to-back, and Drell-Yan scattering to explore the possibility that this allows the region of minimum $K$-factor in Fig. 4 to be fully accommodated with minimal modification to existing fits.

Since the pioneering work in presented in Meng et al. (1996); Graudenz (1997); Nadolsky et al. (1999), there has been a large number of studies on unpolarized SIDIS cross sections Koike et al. (2006); Anselmino et al. (2007, 2014); Signori et al. (2013); Sun and Yuan (2013); Sun et al. (2014); Bacchetta et al. (2017). Unpolarized SIDIS is however, only one component in a broad program of phenomenolical studies where the universality of parton correlation functions plays a central role in testing pictures of nucleon structure Anselmino et al. (2017, 2015); Bacchetta et al. (2015); Scimemi and Vladimirov (2018); Kang et al. (2017, 2015a); Landry et al. (2003); Sun et al. (2017, 2011); Echevarria et al. (2015); Guzzi et al. (2014); Echevarria et al. (2016, 2014); Boglione et al. (2018, 2017); Kang et al. (2015b); Lin et al. (2018); Ye et al. (2017); Bacchetta et al. (2019); Bastami et al. (2018); Bertone et al. (2019). In order for this program to be successful, the precision in the determination of parton correlation functions needs to meet the demands of the increasingly precision of current and future experimental programs  Airapetian et al. (2013); Adolph et al. (2013, 2017); Seidl et al. (2019); Aschenauer et al. (2019); Bradamante (2018). This demands a satisfactory resolution to problems in the $q_{\text{T}}\sim Q$ region for which the pQCD calculation at large $q_{\text{T}}$ is crucial.

Acknowledgements.

We thank J. Owens, J.-W. Qiu, and W. Vogelsang for very useful discussions. We thank R. Sassot for explanations regarding the code in Daleo et al. (2005). T. Rogers’s work was supported by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics, under Award Number DE-SC0018106. This work was also supported by the DOE Contract No. DE- AC05-06OR23177, under which Jefferson Science Associates, LLC operates Jefferson Lab. N. Sato was partially supported by DE-FG-04ER41309, DE-SC0018106 and the DOE Contract No. DE- AC05-06OR23177, under which Jefferson Science Associates, LLC operates Jefferson Lab. B. Wang was supported in part by the National Science Foundation of China (11135006, 11275168, 11422544, 11375151, 11535002) and the Zhejiang University Fundamental Research Funds for the Central Universities (2017QNA3007). J. O. Gonzalez-Hernandez work was partially supported by Jefferson Science Associates, LLC under U.S. DOE Contract #DE-AC05-06OR23177 and by the U.S. DOE Grant #DE-FG02-97ER41028.

Appendix A Identities for three-body partonic phase space

The phase space of three massless particles has the form

[TABLE]

The $k_{1}$ integral can be performed in the center-of-mass frame of $p$ and $q$

[TABLE]

where $\bm{k}_{n-2}$ denotes the first $n-2$ components of the $n-1$ dimensional unit spatial vector in spherical coordinates, and

[TABLE]

In $k_{1}$ phase space, $E$ and $\theta$ are related with observables. The “azimuthal” angles in $\bm{k}_{n-2}$ can be integrated to give

[TABLE]

where $\Omega_{m}$ is the $m$ dimensional angular volume,

[TABLE]

$E$ and $\theta$ are related to Lorentz invariants $u$ and $t$ via

[TABLE]

In terms of $t$ and $u$, the $k_{1}$ phase space becomes

[TABLE]

To simplify the $k_{2}$ and $k_{3}$ phase space, we work in the $k_{2}$ and $k_{3}$ center-of-mass frame where $k_{2}$ and $k_{3}$ take the form of Eqs. (36) and (37). The choice of the spatial orientation of axes also follows the discussion below Eq. (36). In this frame

[TABLE]

Since the scattering amplitudes are independent of the first $n-3$ spatial components of $k_{2}$, the “azimuthal” part of $\mathop{\mathrm{d}\mskip-0.59999mu^{n-1}k_{2}}\nolimits$ can be integrated to give

[TABLE]

Appendix B Phase space integration for $2\to 3$ processes in SIDIS

The main effort of the calculation is in the 2 $\rightarrow$ 3 process $p+q\rightarrow k_{1}+k_{2}+k_{3}$. After partial fractioning (to be discussed Appendix D), each term from the squared amplitudes $\lvert M^{2\to 3}_{g}\rvert^{2}$,$\lvert M^{2\to 3}_{PP}\rvert^{2}$ in Eq. (38) contains at most two Mandelstam variables depending on the angles $\beta_{1}$ and $\beta_{2}$. So, up to overall factors, Eq. (38) takes the form

[TABLE]

where $j$ and $l$ are integers. The coefficients $a$, $b$, $A$, $B$, and $C$ are specific to the squared amplitude and do not depend on $\beta_{1}$ and $\beta_{2}$. They are determined by the spatial orientations of axes in the $k_{2}+k_{3}$ rest frame. Note that Eqs. (36)–(37) do not determine uniquely the components of $p$, $q$, and $k_{1}$. There can still be a rotation in the plane determined by the last two spatial axes. Taking advantage of this freedom, we specify three frames.

• Frame 1:

[TABLE]

• Frame 2:

[TABLE]

• Frame 3:

[TABLE]

where

[TABLE]

There are three cases (see Appendix C for a proof) for the coefficients in the integral in Eq. (67): (1) $a^{2}=b^{2}$ and $A^{2}=B^{2}+C^{2}$; (2) $a^{2}>b^{2}$ and $A^{2}=B^{2}+C^{2}$; (3) $a^{2}=b^{2}$ and $A^{2}>B^{2}+C^{2}$. Integrals in case (3) can be transformed to case (2) by switching to another frame. With multiple frames, the number of integrals to compute can be reduced. In the following we show how the integrals in case (1) and (2) are obtained.

• Case (1): In this case the integral in Eq. (67) has a closed form van Neerven (1986)

[TABLE]

where for simplicity an overall factor is not shown and we choose $a=A=-b=1$, $B=\cos\chi$, and $C=\sin\chi$. In the frames we choose above, $\chi$ can be $\alpha_{1}$ or $\theta_{2}$. Here the signs of the trigonometric function terms in the denominator are chosen to be negative. All other choices can be transformed to this one by applying substitutions $\beta_{1,2}\rightarrow\pi-\beta_{1,2}$, and/or $\chi\rightarrow\pi+\chi$. The derivation of this integral can be found in appendix A of van Neerven (1986).

• Case (2): In this case the integral no longer has a closed form. Following steps similar to appendix A of van Neerven (1986), we arrive at

[TABLE]

where $D>1$. The integral representation of hypergeometric function is used. Note the $I_{j,l}$ notation for integrals with different denominator structures. The result of evaluating Eq. (85) for specific $j$ and $l$ is given in Appendix F. Here $\chi$ can be $\alpha_{3}$ or $\theta_{3}$. The remaining $z$ integral has to be performed by expanding the integrand into series of $\epsilon$, and computing the integral order by order. The subtlety comes from the treatment of the factor $(1-z)^{n/2-l-2}$, which can produce poles for $l\geq 1$. The integral has the form

[TABLE]

where $f(z)$ is a regular function at $z=1$. When $l=1$, the pole term can be made explicit using the identity

[TABLE]

where the “plus” functions are defined, for a function $p(\xi)$, through

[TABLE]

When $l>1$, the integral in Eq. (87) is divergent. Nonetheless, it can be analytically continued. To do this, we write

[TABLE]

Note that, in the first term on the right-hand side of Eq. (90), the power of $1-z$ is increased by 1, with a regular function $f_{1}(z)$. The second term is divergent, but can be analytically continued using $\Gamma$ functions:

[TABLE]

The manipulation in Eq. (90) can be repeated until the power of $1-z$ in the integral becomes $-\epsilon-1$, in which case Eq. (88) can be used. To expand the hypergeometric function in the integrand in Eq. (85), we note that if either $j$ or $l$ is less than 1, the hypergeometric series terminates and the function reduces to a polynomial. For $j=1$ and $l=1$ we use the expansion

[TABLE]

For the case that one of $j$ and $l$ is larger than 1 and the other is at least 1, we use Gauss’s contiguous relations to reduce $j$ or $l$. For example, for $F(1,2,1-\epsilon)$ and $F(2,2,1-\epsilon)$ (with the shorthand notation $F(j,l,k)\equiv\,_{2}F_{1}(j,l,k,w)$), we use

[TABLE]

to get

[TABLE]

Other well-known non-trivial aspects of the computation are the conversion of the expressions in the squared amplitude into forms that allow Eq. (84) or Eq. (85) to be used, which are reviewed in Appendix D, and the algorithm for calculating virtual corrections, which are reviewed in Appendix E. Finally the procedure for combining all Feynman graph calculations consistently into a a factorized cross section is reviewed in Appendix LABEL:a.subtraction.

Appendix C Proof that $a^{2}\geq b^{2}$ and $A^{2}\geq B^{2}+C^{2}$ in Eq. (67)

First, we show that $t_{i}<0$. In the center-of-mass frame of $p$ and $q$,

[TABLE]

where $W=\sqrt{s}$. $k$ can be any one of $k_{i}$, $i=1,2,3$, with the subscript suppressed above. In this frame, $E<\frac{1}{2}\sqrt{s}=\frac{1}{2}W$ and

[TABLE]

The $t_{i}$ are Lorentz invariants and thus $t_{i}<0$ holds in any frame. Since the signs of $t_{i}$ do not change and

[TABLE]

we conclude that $a^{2}\geq b^{2}$ and $A^{2}\geq B^{2}+C^{2}$, where the “$=$” holds when the corresponding variable is $u_{i}$ or $s_{ij}$, and the “$>$” holds for $t_{i}$.

Appendix D Partial fraction algorithm

To perform the three-body phase space integrations as described Sec. B, it is necessary to first convert expressions in the squared amplitude into forms that allow Eq. (84) or Eq. (85) to be used. In this appendix, we review the way this can be automized.

To simplify the discussion, we will refer to angle dependent Mandelstam variables (ADMV):

[TABLE]

The angle independent Mandelstam variables (AIMVs) are:

[TABLE]

A general term in the squared amplitude starts as simple products of AIMVs and ADMVs. A conversion to a useful form is accomplished with a set partial fractions such that each term in the squared amplitude can be expressed as

[TABLE]

with $\alpha_{1}$ and $\alpha_{2}$ being real numbers that could be positive, negative or zero. The quantity in the brackets is some combination of the AIMVs while $X$ and $Y$ are two different types of ADMV that belong to to either $t_{i}$, $u_{i}$ or $s_{ij}$. A term in such a form is called a “reduced” term. Otherwise, we will call it “unreduced”. Example are:

• •

reduced** terms**: (constant), ($t_{2}^{-2}u_{3}^{-1}$), ($s_{12}^{-2}$),($u_{2}$)

• •

unreduced** terms**: ($t_{2}^{-2}t_{3}^{-1}$), ($u_{2}u_{3}$), ($t_{2}^{-2}u_{3}^{-1}s_{12}^{-1}$)

We need the following further sets of definitions:

• •

t-type, u-type, s-type Mandelstam variables: An ADMV is referred to as t-type if it is $t_{2}$ or $t_{3}$, and similarly for u-type and s-type.

• •

The ADMVs are not independent variables; it is possible to construct linear relations relating them. This is achieved by squaring different rearrangements of the momentum conservation equation:

[TABLE]

• Two ADMV relations (2ARs): One such type of relation is between two ADMVs of the same type. These are obtained by moving one of $p$, $q$, and $k_{1}$ in Eq. (110) to the right-hand side and squaring the equation. For instance, moving $p$ to the right and squaring gives

[TABLE]

relating two u-type ADMVs. Similarly, moving $q$ or $k_{1}$ to the right gives the relations for $t$-type or $s$-type ADMVs respectively. There are three 2ARs. The other two are:

[TABLE]

• Three ADMV relations (3ARs): Another type of relation is between three ADMVs, each from a different ADMV type. We can get such relations for any three ADMVs of different types. To do this notice that exchanging the positions of $k_{2}$ or $k_{3}$ with $k_{1}$ in Eq. (110) and squaring will give one of such relation. For instance, squaring $p+q-k_{2}=k_{1}+k_{3}$:

[TABLE]

a relations between $s_{13}$, $t_{2}$, and $u_{2}$. After one such relation is obtained, the others are obtained by replacing the ADMVs in that relation using three 2ARs.

Now the partial fraction is implemented in two steps:

• Step 1: Consider a term $T_{1}$ in the squared amplitude. The first step is to reduce the number of ADMVs in the denominator to two or fewer, and to ensure that no two ADMVs of the same type appear in the denominator. This is done in two substeps:

• Substep 1: First, separate any ADMVs of the same type in the denominator. If two ADMVs of the same type appear in the denominator of $T_{1}$, use a 2AR for the two ADMVs to write a factor of unity in the form

[TABLE]

and write

[TABLE]

with $c_{j}=\pm 1$. Expanding the right side of Eq. (116) produces a sum of new terms. Each will have one power lower of these two ADMVs in its denominator. Further multiplication by unity factors like Eq. (115) can be repeated until in all terms only one of the two ADMVs appears in any denominator. As an example, consider

[TABLE]

$t_{2}$ and $t_{3}$ both appear in the denominator, so we will use Eq. (113) to write Eq. (115):

[TABLE]

Then

[TABLE]

• Substep 2: After Substep 1 is done for all three types, the leftover terms have denominators with at most three ADMVs, all from different types. For terms with two or fewer ADMVs in the denominator, nothing more needs to be done in this step. For terms with three ADMVs in the denominator, one may write a 3AR in the form of Eq. (115) and reduce the number of ADMVs in the denominator to two or fewer in a way similar to Substep 1. For example, say that

[TABLE]

Then it is possible to eliminate $s_{13}$, $u_{2}$ and $t_{2}$ by using Eq. (114) to construct Eq. (115). Then,

[TABLE]

• Step 2: After Step 1 is completed, all terms have two or fewer ADMVs in the denominator. To get to the final form in Eq. (109), we need to make sure no third ADMV appears in the numerator. This is done by noticing that the ADMVs in the numerator can be written in terms of the ADMVs of the denominator using 2ARs and 3ARs. For instance, if one term has $t_{2}u_{3}$ in the denominator, then any t-type ADMV in the numerator can be replaced by $t_{2}$, and u-type by $u_{3}$, using 2ARs. Any s-type ADMV can be replaced by a linear combination of $t_{2}$ and $u_{3}$ using the corresponding 3AR involving $t_{2}$ and $u_{3}$. This puts $T_{1}$ in the form of Eq. (109). (Recall that $\alpha_{1}$ and $\alpha_{2}$ can be negative or zero.) This example had two ADMVs in the denominator. If the number of ADMVs in the denominator is 1, one can express the ADMVs in the numerator in terms of two ADMVs of different types, where one is the ADMV in the denominator and the other is chosen randomly. If the number of ADMVs in the denominator is 0, then the ADMVs in the numerator can be expressed by any two randomly chosen ADMVs of different types. There are six ADMVs, but the four relations from Eqs. (111)–(114) eliminate four in terms of the other two.

Appendix E Virtual Contributions

At order $O(\alpha_{s}^{2})$, virtual corrections to the partonic cross section involve only one-loop integrals. Comprehensive reviews of the methods for this type of calculation can be found in Berger and Forde (2010); Britto (2011); Ellis et al. (2012). In this article, we use the traditional Passarino-Veltman (PV) approach Passarino and Veltman (1979), for which we closely follow the notation of Ellis et al. (2012).

Our $O(\alpha_{s}^{2})$ virtual contributions are the interference terms between tree-level $2\to 2$ amplitudes and $2\to 2$ amplitudes with an additional virtual loop, averaged(summed) over initial(final) states. The structure of each such interference term is

[TABLE]

The denominators contain massless propagators of the form $d_{i}=(l+\Delta_{i})^{2}+i0$, where $\Delta_{i}$ depends only on external momenta, which we will call $v_{i}$. That is, $v_{i}\in\{p,q,k_{1}\}$ or combinations thereof. The trace in (123) contains Dirac gamma matrices contracted with either the external momenta $v_{i}$, or with the loop momentum $l$. Both IR and UV divergences are handled by standard dimensional regularization techniques. After these steps, the numerator in Eq (123) is a collection of terms with Lorentz invariant products of loop and external momenta. All integrals needed to calculate the virtual corrections in Figure 2 can be written in one of the following ways:

[TABLE]

where higher rank integrals cannot appear since, in every case, virtual correction diagrams at one-loop involve at least one gluon propagator. We have left the Feynman prescription implicit in the denominators in Eq. (124), and partonic cross section calculations are done in massless QCD. These tensor integrals can be written in terms the scalar box ($D_{0}$), triangle ($C_{0}$) and bubble ($B_{0}$) integrals using the Passarino-Veltman reduction procedurePassarino and Veltman (1979).

As an example, consider the steps to reduce $C^{\alpha\beta}$. Symmetry properties allow $C^{\alpha\beta}$, $C^{\beta}$ and $B^{\beta}$ to be written in terms of form factors ${\rm B}_{1}(v_{1})$, ${\rm C}_{1}(v_{1},v_{2})$, ${\rm C}_{2}(v_{1},v_{2})$, and ${\rm C}_{ij}(v_{1},v_{2})$:

[TABLE]

Defining the Graham matrix $\mathbf{G}_{i,j}=v_{i}\cdot v_{j}$ and the row vectors $\vec{{\rm C}}^{(1)}_{i}={\rm C}_{i1}$, $\vec{{\rm C}}^{(2)}_{i}={\rm C}_{i2}$, the contraction of both sides of Eq. (125c) with $v_{i}{}_{\alpha}$ and $g_{\alpha\beta}$ respectively gives

[TABLE]

where $\mathbf{R}$ is defined by the relations

[TABLE]

Here we have used the notation $u=v_{1}+v_{2}$. The results for $R_{ij}$ above are obtained by contracting $v_{i}{}_{\alpha}$ with the right-hand side of Eq. (125c) and using the property that the products $v_{i}\cdot l$ can always be written in terms of the denominators $d_{i}=(l+v_{1}+...+v_{i+1})^{2}$ as

[TABLE]

This allows $v_{i}{}_{\alpha}\,C^{\alpha\beta}$ to be reduced to rank-1 integrals, whereupon expansion according to Eqs. (125b), lead to Eqs. (E). Solving Eq. (126) gives ${\rm C}_{ij}$ and ${\rm C}_{00}$ in terms of form factors corresponding to lower rank tensor integrals. In a similar way, one may reduce residual ${\rm C}_{i}$ and ${\rm B}_{1}$ form factors so that only the scalar integrals $C_{0}$ and $B_{0}$ appear at the end of the reduction procedure. Equations analogous to (126) and (E) for all cases in Eq. (124) are provided in Appendix A of Ellis et al. (2012).

After the reduction procedure, all virtual contributions from Figure 2 are in terms of scalar integrals, whose values depend only on Lorentz invariants constructed from momenta that appear in the denominators. The complete set of scalar integrals are in  Ellis and Zanderighi (2008). In our calculation, all of these contain single or double pole singularities.

The singular behavior of $D_{0}$ and $C_{0}$ corresponds to soft and collinear divergences that exactly cancel the soft and collinear singularities in the $2\to 3$ calculation of the corresponding channel.

For our computations, we only need expressions for a reduced number of cases for $C_{0}$ and $D_{0}$, following the notation of  Ellis and Zanderighi (2008): $C_{0}$ divergent type 1 and type 2, $D_{0}$ divergent type 2

The usual UV singularities introduced by virtual corrections ultimately are all produced by $B_{0}$-type integrals. UV divergent tadpole integrals $A_{0}$ also appear but are zero in a massless theory. Thus, keeping track of $B_{0}$ integrals accounts for all UV singularities, which cancel in the sum of virtual graphs. Self-energy diagrams on external legs enter for each channel via the corresponding factors of field strength normalization $Z_{i}^{1/2}$, as prescribed by the LSZ theorem.

At order $O(\alpha_{s}^{2})$, $Z_{i}$ deviates from unity only for $2\to 2$ partonic scattering. Denoting the sum of amputated diagrams for leading order, virtual, and real emission by ${\cal M}_{LO}^{(a)}$, ${\cal M}_{V}^{(a)}$ and ${\cal M}_{R}^{(a)}$, the square-modulus amplitude for partonic scattering is

[TABLE]

Explicit solutions for scalar integrals given in Ellis and Zanderighi (2008) hold when all of the relevant Lorentz invariants are space-like. When needed, we perform analytic continuation on the logarithms, by restoring the Feynman prescription $\log\left(\frac{a}{b}\right)\longrightarrow\log\left(a+i0\right)+\log\left(\frac{1}{b+i0}\right)$. For the dilogarithms in our calculation, we find the following relation useful

[TABLE]

where both $a$ and $b$ are positive.

Appendix F Angular integrals with a virtual photon

For completeness we provide a complete list of the relevant angular integrals. The method to compute these integrals is described in detail in Sec Sec. B. The alternatives for some of the integrals here can be found from appendix C of Ref. Beenakker et al. (1989). But there the list of integrals is computed for heavy quark production process and is not adequate for SIDIS. To compactify notation, we define

[TABLE]

The integrals are computed up to the needed powers in $\epsilon$. The expressions for $I_{j,l}$ from Eq. 85 needed for order $\alpha_{s}^{2}$ SIDIS are then

[TABLE]

Certain order $\epsilon$ terms are not needed and may be dropped in calculations at $O\left(\alpha_{s}^{2}\right)$.