The spin-dependent quark beam function at NNLO

TL;DR

This paper computes the polarized quark beam function at NNLO in QCD, completing the ingredients for NNLO analysis of polarized collision observables and providing a cross-check for unpolarized results.

Contribution

It presents the first calculation of the polarized quark beam function at NNLO, filling a key gap for precise theoretical predictions in polarized high-energy collisions.

Findings

First independent check of unpolarized quark beam function

Provides detailed derivation of polarized beam function at NNLO

Enables NNLO analysis of polarized N-jettiness observables

Abstract

We calculate the beam function for longitudinally-polarized quarks through next-to-next-to-leading order (NNLO) in QCD perturbation theory. This is the last missing ingredient needed to apply the factorization theorem for the $N$-jettiness event-shape variable in polarized collisions through the NNLO level. We present all technical details of our derivation. As a by-product of our calculation we provide the first independent check of the previously-obtained unpolarized quark beam function. We anticipate that our result will have phenomenological applications in describing data from polarized collisions.

Equations

$${\mathcal{T}}_N=\sum _k{\text{min}}_i\left\{\frac{2p_i\cdot q_k}{Q_i}\right\}.$$

$$\frac{d\sigma }{d{\mathcal{T}}_N}=H\otimes B\otimes S\otimes \left[\prod _n^N{J}_n\right]+\cdots .$$

$$\frac{d{\sigma }_{LL}}{d{\mathcal{T}}_N}=\Delta H\otimes \Delta B\otimes S\otimes \left[\prod _n^N{J}_n\right]+\cdots .$$

$$\Delta B_i=\sum _j\Delta {\mathcal{I}}_{ij}\otimes \Delta f_j\left[1+\mathcal{O}\left(\frac{{\Lambda }_{\mathrm{QCD}}^2}{{\mathcal{T}}_N}\right)\right].$$

$$\Delta B_q(t,x,\mu )=$$

$$\Delta B_{\bar{q}}(t,x,\mu )=$$

$${\chi }_n(y)=W_n^{†}(y){\xi }_n(y),$$

$$W_n(y)=\left[\sum _{\mathrm{perms}}\exp \left(-\frac{g}{{\overline{\mathcal{P}}}_n}\bar{n}\cdot A_n(y)\right)\right].$$

$$\Delta B_{ij}(t,z,\mu )=\sum _k\Delta {\mathcal{I}}_{ik}(t,z,\mu )\otimes \Delta f_{kj}\left(z\right)\equiv \sum _k{\int }_z^1\frac{d{z}^{\prime }}{z^{\prime }}\Delta {\mathcal{I}}_{ik}(t,z^{\prime },\mu )\Delta f_{kj}\left(\frac{z}{z^{\prime }}\right).$$

$$\Delta B_{qq}^{(0)}(t,z,\mu )=⟨q_n(p),+|\theta (\omega ){\bar{\chi }}_n(0)\delta (t-\omega {\hat{p}}^{+})\frac{/\bar{n}{\gamma }_5}{2}[\delta (\omega -{\overline{\mathcal{P}}}_n){\chi }_n(0)]|q_n(p),+⟩=\delta (t)\delta (1-\omega /p^{-}).$$

$${\mathcal{I}}_{qq}^{(0)}(t,z,\mu )={\mathcal{I}}_{\bar{q}\bar{q}}^{(0)}(t,z,\mu )=\delta (t)\delta (1-z).$$

$${\mathcal{I}}_{qg}^{(0)}(t,z,\mu )={\mathcal{I}}_{gq}^{(0)}(t,z,\mu )=0.$$

$$\Delta B_{ij}^{bare}(t,z)=\int d{t}^{\prime }{Z}_i(t-t^{\prime },\mu )\Delta B_{ij}(t^{\prime },z,\mu ),$$

$$\mu \frac{d}{d\mu }\Delta B_{ij}(t^{\prime },z,\mu )=\int d{t}^{\prime }{\gamma }_B^i(t-t^{\prime },\mu )\Delta B_{ij}(t^{\prime },z,\mu ),$$

$${\gamma }_B^i(t,\mu )=-\int d{t}^{\prime }(Z_i{)}^{-1}(t-t^{\prime },\mu )\mu \frac{d}{d\mu }{Z}_i(t^{\prime },\mu ).$$

$$\int d{t}^{\prime }(Z_i{)}^{-1}(t-t^{\prime },\mu )Z_i(t^{\prime },\mu )=\delta (t).$$

$$\Delta B_{ij}=$$

$$Z_i=$$

$$\Delta {\mathcal{I}}_{ij}=$$

$$\Delta f_{ij}=$$

$$\Delta B_{ij}^{(0)}(t^{\prime },z,\mu )={\delta }_{ij}\delta (t^{\prime })\delta (1-z),Z_i^{(0)}(t-t^{\prime },\mu )=\delta (t-t^{\prime }),$$

$$\Delta B_{ij}^{bare(1)}(t,z)=\Delta B_{ij}^{(1)}(t,z,\mu )+Z_i^{(1)}(t,\mu ){\delta }_{ij}\delta (1-z).$$

$$\Delta B_{ij}^{bare(2)}(t,z)=\Delta B_{ij}^{(2)}(t,z,\mu )+Z_i^{(2)}(t,\mu ){\delta }_{ij}\delta (1-z)+\int d{t}^{\prime }{Z}_i^{(1)}(t-t^{\prime },\mu )\Delta B_{ij}^{(1)}(t^{\prime },z,\mu ).$$

$$\{{\gamma }_5,{\tilde{\gamma }}_{\mu }\}=0,[{\gamma }_5,{\hat{\gamma }}_{\mu }]=0.$$

$$\Delta B=\Delta \tilde{I}\otimes \Delta \tilde{f},$$

$$\Delta {\mathcal{I}}_{q_i{q}_j}(z)=$$

$$\Delta {\mathcal{I}}_{q_i{\bar{q}}_j}(z)=$$

$$\Delta I_{qq}^{(V)}=I_{qq}^{(V)},\Delta I_{q\bar{q}}^{(V)}=-I_{q\bar{q}}^{(V)}.$$

$$\Delta B=\left(\Delta \tilde{I}\otimes {\bar{Z}}^5\right)\otimes \left(Z^5\otimes \Delta \tilde{f}\right)=\Delta I\otimes \Delta f$$

$$\Delta I=\Delta \tilde{I}\otimes {\bar{Z}}^5,\Delta f=Z^5\otimes \Delta \tilde{f}$$

$$\Delta {\tilde{f}}_{ij}^{(1)}(z)=$$

$$\Delta {\tilde{f}}_{ij}^{(2)}(z)=$$

$${\beta }_0=\frac{11C_A-4T_R{N}_F}{3}.$$

Full text

I Introduction

The precision of modern hadron collider experiments has reached an impressive level. At the Large Hadron Collider (LHC), sub-percent experimental errors have been obtained in numerous measurements that span orders of magnitude in energy. Theoretical advances have matched the pace of this experimental progress. New calculations and simulation tools have been developed over the past decade that have made possible precision comparisons between theory and data that reach the percent level for some measurements. These achievements together test the high-energy predictions of the Standard Model more stringently than ever before.

There is a desire to extend such rigorous quantitative comparisons to polarized collisions as well. New data is becoming available for polarized collisions with direct implications for measurements of fundamental importance, such as the determination of the polarized gluon content of the proton at the Relativistic Heavy Ion Collider (RHIC) Adamczyk:2014ozi . Several other experiments have released measurements of the double-spin asymmetry in polarized collisions Anthony:1999ac ; Airapetian:2010ac . The theoretical description of this data using perturbative QCD predictions at next-to-leading order (NLO) in the strong coupling constant is not entirely satisfactory Hinderer:2017ntk . More precise data is expected at a future electron-ion collider (EIC). The EIC is expected to be a unique laboratory in which to discover novel QCD phenomena and to determine the spin structure of the proton with an accuracy currently not achievable. In order to achieve these goals perturbative QCD effects will need to be disentangled from new non-perturbative phenomena and proton structure. The situation will be very similar to the current one at the LHC, where higher-order perturbative QCD effects must be distinguished from potential signals of beyond-the-Standard Model phenomena.

In order to quantitatively describe and analyze data from polarized collisions, more theoretical advances are necessary. As mentioned, NLO perturbative QCD fails to properly describe double-spin asymmetries appearing in deeply-inelastic scattering (DIS) in kinematic regions where non-perturbative effects are expected to be suppressed Hinderer:2017ntk . The situation may be compared to the case of inclusive jet production in unpolarized DIS that has received attention recently. The NLO corrections were found to be large, increasing the leading-order result by nearly a factor of two Hinderer:2015hra . The perturbative series stabilized only after the inclusion of next-to-next-to-leading order (NNLO) corrections Abelof:2016pby . This suggests that calculations at this order in perturbation theory may be necessary for a proper description of current polarized hadron-collider data. As the experimental errors are expected to decrease at an EIC, it is expected that NNLO theory will become an integral part of the EIC physics program.

One very successful approach to higher-order calculations for unpolarized collisions at the LHC has been $N$-jettiness subtraction Boughezal:2015dva ; Gaunt:2015pea . This method has been used to compute the NNLO corrections for a host of jet production processes Boughezal:2015dva ; Boughezal:2015dva ; Boughezal:2015aha ; Boughezal:2015ded ; Boughezal:2016yfp ; Boughezal:2016dtm ; Boughezal:2016isb ; Campbell:2016lzl , as well as numerous color-singlet production processes Boughezal:2016wmq . It has been extended to handle inclusive jet production in electron-nucleon collisions as well Abelof:2016pby . $N$-jettiness subtraction is based on the event-shape variable ${\cal T}_{N}$ Stewart:2010tn ; Stewart:2009yx ; Stewart:2010qs :

[TABLE]

The subscript $N$ denotes the number of jets desired in the final state, and is an input to the measurement. The $q_{k}$ denote the four-momenta of any final-state radiation, while the $p_{i}$ denote the momenta of the initial-state hard partons and any final-state jets. The $Q_{i}$ are dimensionful variables that characterize the hardness of the beam-jets and final-state jets. The cross section factorizes in the limit where ${\cal T}_{N}$ is less than any other hard scale in the problem Stewart:2010tn ; Stewart:2009yx ; Stewart:2010qs . Schematically, the form of the cross section in this limit becomes

[TABLE]

Here, $H$ describes the effect of hard radiation, $S$ describes the soft radiation, and $J_{n}$ contains the radiation collinear to a final-state jet. We note that the $J_{n}$ can be straightforwardly replaced by a fragmentation function if a final-state hadron is instead measured. $B$ encodes the effect of radiation collinear to an initial beam direction. We have assumed a single hadronic beam as for DIS; this formula would contain two beam functions for proton-proton collisions, or none for $e^{+}e^{-}$ collisions. Depending on the observable and process under consideration, only a subset of the other terms may be present. The ellipsis denote power-suppressed terms which become negligible for ${\cal T}_{N}\ll Q_{i}$. There has been recent progress in understanding these terms for color-singlet production Moult:2016fqy ; Boughezal:2016zws ; Moult:2017rpl . The derivation of this factorization theorem relies heavily upon the machinery of Soft-Collinear Effective Theory (SCET) scet . As a NNLO calculational technique, $N$-jettiness subtraction works by partitioning the phase space using a cut ${\cal T}_{N}^{cut}$ with ${\cal T}_{N}^{cut}\ll Q_{i}$, using the factorization theorem of Eq. (2) below this cut, and noting that the cross section above the cut becomes a simpler NLO cross section. We refer the reader to the original papers for more details Boughezal:2015dva ; Gaunt:2015pea . We note that the study of $N$-jettiness has intrinsic interest besides its use as a NNLO subtraction scheme. It is widely used in study of jet substructure through its $N$-subjettiness incarnation Thaler:2010tr , and has been proposed as a measure of the nuclear medium in electron-nucleus collisions Kang:2012zr ; Kang:2013wca . A key aspect of the usefulness of $N$-jettiness at the LHC is our ability to calculate to high orders in the QCD perturbative expansion the objects that appear in the factorization theorem of Eq. (2). The beam function $B$ Gaunt:2014xga ; Gaunt:2014cfa , the jet function $J_{n}$ Becher:2006qw ; Becher:2010pd and the soft function ${\cal S}_{N}$ for jets Boughezal:2015eha and for the massive case Li:2016tvb are all known to the NNLO level, as are the hard functions for many processes of interest.

Given the success of this framework in describing a host of data at the LHC and at other unpolarized colliders, and the need for higher-order corrections to better describe current data from polarized collisions as well as expected future data from an EIC, it is interesting to extend the $N$-jettiness framework to cover polarized collisions. As a concrete example we consider the double spin asymmetry in lepton-proton collisions, ${\cal A}_{LL}$. This observable begins at the leading order in the twist expansion. It is straightforward to write down the analogous factorization theorem in the low-${\cal T}_{N}$ limit for the relevant polarized cross section that enters the numerator of this asymmetry111The possible contribution of perturbative Glauber modes to this factorization theorem is expected to occur at higher orders in the strong coupling constant than the NNLO level of interest here Zeng:2015iba .:

[TABLE]

Here, the $\Delta$ denotes the polarization dependence caused by taking the appropriate difference of helicities needed to obtain the double-spin asymmetry. We note that the soft function $S$ and the jet/fragmentation function $J_{n}$ are unchanged upon considering polarized collisions. Furthermore, the hard function is obtained from the virtual corrections to the scattering process under consideration. Since these are generally computed for the separate helicity states, $\Delta H$ is known for most processes of interest. Only the polarized beam function $\Delta B$ is not known at the NNLO level. In analogy to the beam function for unpolarized collisions, the polarized beam function is a non-perturbative object that can be matched to the polarized parton distribution functions Stewart:2010tn ; Stewart:2009yx ; Stewart:2010qs :

[TABLE]

The $\Delta{\cal I}_{ij}$ are perturbatively calculable matching coefficients, while the $\Delta f_{j}$ are the standard polarized PDFs. The $i$ and $j$ are parton labels. From this expression it is also apparent that studying the polarized beam function may reveal aspects of the polarized PDFs.

It is our goal in this manuscript to calculate the polarized quark beam function matching coefficients through the NNLO level. This is the last missing ingredient needed to bring the theoretical status of polarized collisions at leading twist to the same level as for unpolarized collisions. It will allow for NNLO calculations in polarized collisions, and will be a necessary ingredient in extending the global fit of polarized proton structure deFlorian:2014yva to the NNLO level. The beam function can be thought of as a generalized PDF, where the operators appearing in the definition are separated along both light-cone directions. In this sense our calculation of $\Delta B$ represents the first NNLO calculation of such an object for polarized collisions. We discuss all relevant details of our calculation, including our treatment of $\gamma_{5}$ in dimensional regularization. As a by-product of our calculation we provide the first independent check of the NNLO unpolarized quark beam function. We find complete agreement with the previous result for this quantity Gaunt:2014xga . However, we do not find any need for an additional regulator beyond dimensional regularization in our calculation, in apparent contrast to this previous result.

Our paper is organized as follows. We define the spin-dependent beam function and establish our notation in Section II. In Section III we discuss the renormalization and matching needed to convert the bare results into the renormalized matching coefficients needed in Eq. (4). We present the details of our NLO calculation in Section IV. The details of our NNLO derivation are presented in Section V. Finally, we conclude in Section VI.

II Setup and definition of the spin-dependent beam function

We begin by establishing our notation and defining the spin-dependent beam function. We will use the standard light-cone vectors $n^{\mu},\bar{n}^{\mu}$ with $n^{2}=\bar{n}^{2}=0$ and $n\cdot\bar{n}=2$. Any four vector can be written in terms of these directions as $p^{\mu}=(p^{+},~{}p^{-},~{}p_{\perp}^{\mu})=(n\cdot p,~{}\bar{n}\cdot p,~{}p_{\perp}^{\mu})$.

The longitudinal spin-dependent beam function represents the difference of beam functions with positive and negative helicity in a parent proton with positive helicity. Focusing on the quark and anti-quark beam functions, we can define them as the proton matrix element of quark operators:

[TABLE]

Here, the “$+$” represents the positive helicity of the proton, $x=\omega/P^{-}$ is the momentum fraction of the proton carried by the parton that enters the hard scattering, the $\delta(\omega-\overline{\mathcal{P}}_{n})$ operator constrains the total minus momentum of the composite quark/gluon field to $\omega$, and $\delta(t-\omega\hat{p}^{+})$ sets the total plus momentum of all initial state radiations to $t/\omega$. $\chi$ is the composite quark operator

[TABLE]

where $\xi_{n}$ is the n-collinear quark field, and $W_{n}$ is the Wilson line

[TABLE]

We refer the reader to the SCET literature for more details on the operators that appear in these definitions scet . Our definition of the polarized beam function follows the definition of the polarized PDF Manohar:1990jx with the appearance of the additional plus momentum component $t/\omega$. We note that $t$ represents the beam-sector contribution to the measured $N$-jettiness ${\cal T}_{N}$.

As discussed in the introduction the polarized beam function is a non-perturbative quantity that can be matched to polarized PDFs, in analogy to the unpolarized beam function Stewart:2010tn ; Stewart:2009yx ; Stewart:2010qs . In order to calculate the matching coefficients, we replace the proton state by $n$-collinear quark and gluon states with momentum $p=(0,p^{-},0)$. The desired matching coefficients are unchanged upon making this replacement. With this substitution the matching equation takes the form

[TABLE]

The quantity $\Delta f_{kj}$ is the distribution function for a parton of flavor $k$ within another parton of flavor $j$ (we have replaced the proton state by $j$). The $\Delta B_{ij}$ are the polarized beam functions with this replacement for the PDFs. The tree level diagram for quark beam function with an external quark is

[TABLE]

The tree-level matching coefficient is therefore

[TABLE]

Inserting instead a gluon in place of the the initial-state proton leads to a vanishing result, allowing us to conclude that

[TABLE]

III Renormalization and matching

At higher orders in the strong coupling constant we must renormalize the beam function, and also perform the matching to the PDFs. The matching equation has already been presented in Eq. (9). The bare and renormalized beam functions are related through the renormalization constants $Z_{i}$:

[TABLE]

where the bare beam function depends on the renormalized $\overline{\rm MS}$ coupling $g$, and the renormalization constants $Z_{i}$ are defined to remove UV divergences in the bare beam function. Taking derivatives on both sides of Eq. (13) with respect to $\mu$, one can derive the renormalization group equations (RGEs) for the polarized beam functions:

[TABLE]

where the anomalous dimension for the quark beam function in the $\rm\overline{MS}$ scheme is defined as

[TABLE]

The inverse of $Z_{i}$ is defined as

[TABLE]

As we will see later from explicit calculations at NLO and NNLO, the renormalization constant is the same in the polarized and unpolarized cases. This indicates that the RGEs for polarized beam functions follow exactly the same form as in the unpolarized case.

To facilitate the expansion in the strong coupling constant we introduce separate expansions for each of the quantities that appear in our result:

[TABLE]

The different choices of two and four in these expansions match the typical conventions in the literature for the various objects. At order $\alpha_{s}$, using the fact that

[TABLE]

we can derive the following relation between the renormalized and bare beam functions at NLO:

[TABLE]

Similarly, we can expand Eq. (13) to obtain the analogous relation for the NNLO beam functions:

[TABLE]

We will use standard $\overline{{\rm MS}}$ renormalization, so that the $Z_{i}^{(n)}$ renormalization constants will contain only poles in the dimensional regularization parameter $\epsilon=(4-d)/2$, where $d$ is the space-time dimension.

An important technical issue to discuss is the treatment of $\gamma_{5}$ in $d$-dimensions. Several consistent schemes have been proposed for this purpose. We use the HV scheme tHooft:1972tcz ; Breitenlohner:1977hr , in which $\gamma_{5}$ maintains its 4-dimensional definition: $\gamma_{5}=i\gamma^{0}\gamma^{1}\gamma^{2}\gamma^{3}$. Denoting 4-dimensional quantities with a tilde and $\epsilon$-dimensional ones with a hat, this definition leads to the following commutation and anti-commutation rules for $\gamma_{5}$:

[TABLE]

These rules are easy to track in standard algebraic manipulation programs. We use Tracer Jamin:1991dp to implement these rules, together with in-house routines written in Form Kuipers:2012rf as a cross-check.

The use of the HV scheme necessitates an additional transformation in order to obtain the standard $\overline{\text{MS}}$ factorization scheme for the PDFs, as is well known in the literature Gordon:1993qc ; Vogelsang:1995vh ; Vogelsang:1996im ; Vogt:2008yw ; Moch:2008fj ; Moch:2014sna . Switching to a matrix notation in parton flavors, the beam function computing using Eq. (21) can be written as

[TABLE]

where the tilde represents the results before scheme transformation. We decompose the matching coefficients into singlet and non-singlet pieces in analogy to the usual decomposition performed for the splitting functions:

[TABLE]

The scheme transformation can be derived by demanding helicity conservation for massless quarks, which relates the polarized and unpolarized splitting functions. These requirements naturally extend to the matching coefficients, leading to the relations

[TABLE]

The naively-computed $\Delta\tilde{I}$ do not satisfy these constraints. We can restore these relations by the transformations

[TABLE]

where

[TABLE]

with $Z^{5}$ the scheme transformation matrix, and $\bar{Z}^{5}$ is its inverse which satisfies $\bar{Z}^{5}\otimes Z^{5}=1$. More details on the transformation matrix are given later in this manuscript.

Having established the factorization scheme transformation and the expansion of the renormalization condition through the necessary NNLO order, we now consider the matching condition. To obtain the matching coefficients we replace the proton states in Eq. (5) by perturbative quark or gluon states. With this replacement the polarized PDFs can be calculated as an expansion in $\alpha_{s}$, taking on the familiar forms

[TABLE]

where $\beta_{0}$ is the usual leading-order QCD beta function,

[TABLE]

The polarized splitting functions $\Delta\tilde{P}_{ij}$ needed in our calculation are defined in the Appendix. Expanding the matching equation in terms of $\alpha_{s}$, we can derive the NLO matching coefficient in terms of the renormalized beam function and PDFs:

[TABLE]

Similarly, we can derive the matching coefficient at order $\alpha_{s}^{2}$

[TABLE]

The convolutions required for the calculation were computed with the mathematica package MT Hoeschele:2013gga . Notice that the PDFs of Eq. (III) and the matching coefficients of Eqs. (30, 31) in our matching calculations still require the scheme transformation of Eq. (27). Performing this transformation, we obtain the physical results which restore the helicity conservation equations as shown in Eq. (25).

As discussed in Ref. Stewart:2010qs and reviewed above, the beam function satisfies a renormalization group equation that allows the logarithmic dependence of the matching coefficients on the renormalization scale $\mu$ to be derived. Solving this equation allows us to predict all logarithmically-enhanced terms in $t$ in terms of known anomalous dimensions:

[TABLE]

The non-cusp anomalous dimension for the quark beam function in $\overline{\rm MS}$ is the same as that in unpolarized case, and reads

[TABLE]

The cusp anomalous dimension is

[TABLE]

$\mathcal{L}_{n}$ is the standard plus distribution, defined as

[TABLE]

The only terms to determine are the coefficients of the scale-independent $\delta(t)$ contributions, which we label as $\Delta{I}^{(1)}_{ij}(z)$ and $\Delta{I}^{(2)}_{ij}(z)$. In our calculation we derive all the scale-dependent terms as well in order to check that our results satisfy Eqs. (32) and (33).

IV Calculation at next-to-leading order

We begin by discussing the calculation of the matching coefficients at NLO. At this order, there are two matching coefficients to consider, ${\cal I}^{(1)}_{qq}$ and ${\cal I}^{(1)}_{qg}$. The first can be obtained by setting the external proton to a quark, while the second can be obtained by setting the external proton to a gluon. We note that ${\cal I}^{(1)}_{\bar{q}\bar{q}}={\cal I}^{(1)}_{qq}$ and ${\cal I}^{(1)}_{\bar{q}g}={\cal I}^{(1)}_{qg}$, while ${\cal I}^{(1)}_{q\bar{q}}=0$. To perform our calculation we work in light-cone gauge with $n\cdot A=0$. This has the effect of setting all diagrams with gluons emitted from Wilson lines to zero. In this gauge there is a single diagram contributing to each channel. These are shown in Fig. 1. The cut particles in these diagrams are those satisfying the on-shell constraints.

We use the quark channel as an example to show the details of our calculation. In light-cone gauge, only the left diagram of Fig. 1 contributes to the quark beam function. This diagram can be calculated by using the standard QCD Feynman rules,

[TABLE]

where $\mathcal{P}_{R}=(1+\gamma_{5})/2$ is the spin projection operator. As mentioned previously we use light-cone gauge for the gluon to avoid extra diagrams with a gluon radiated from gauge link, leading to the propagator numerator

[TABLE]

The final-state phase space at NLO can be written as

[TABLE]

To obtain the expression in the second line we have separated the final-state phase space into a 4-dimensional part and an $\epsilon$-dimensional part, where $\hat{k}_{\perp}$ denotes the $\epsilon$-dimensional momentum component. We have also used the delta functions to simplify the expression.

At this point the calculation of the trace in Eq. (38) and the integration over the phase space in Eq. (40) are straightforward. We derive the following bare NLO beam functions:

[TABLE]

Expanding this result to $\mathcal{O}(\epsilon^{0})$ (we note that higher-order terms in $\epsilon$ will be needed when using these expressions in the NNLO calculation), we obtain

[TABLE]

Upon rewriting the coefficient of the $1/\epsilon$ pole in terms of the polarized splitting function $\Delta P_{qq}$ that appears in the matching, we can read off the renormalization constant from Eq. (43), which is the same as that in the unpolarized case:

[TABLE]

We note that there is an overall sign difference in our result for $Z_{q}^{(1)}$ as compared to other results in the literature Ritzmann:2014mka . This is simply due to a different relation used here between the bare beam function and the renormalized one. We believe that the choice we make here is the more conventional one. Substituting the renormalization constant into Eq. (19), we can derive the renormalized quark beam function:

[TABLE]

where the polarized splitting functions can be found in the Appendix. The last step is to match to the polarized PDFs using Eqs. (30) and (III). We note that the residual $1/\epsilon$ poles cancel, and we reproduce the logarithmic terms in Eq. (32) upon using the anomalous dimensions presented in Section III and the polarized splitting functions in the Appendix. This serves as a check on our calculation.

At this point we must match to PDFs and perform the scheme transformation to obtain standard $\overline{\text{MS}}$-scheme expressions. Referring to the Appendix, the necessary transformation is

[TABLE]

We can then derive the unknown coefficients of the scale-independent $\delta(t)$ structures in the matching coefficients:

[TABLE]

We note that the last term in Eq. (49) arises from the $z_{qq}^{(1)}$ correction to the violation of helicity conservation when using the HV scheme to deal with $\gamma_{5}$. We note that these matching coefficients are in agreement with a recent NLO calculation of the matching coefficients for transverse-momentum dependent PDFs Gutierrez-Reyes:2017glx .

V Calculation at next-to-next-to-leading order

Having established the matching coefficients at NLO we move onto the NNLO calculation. We begin by decomposing the $\delta(t)$ contributions to the matching coefficients into singlet and non-singlet pieces, in analogy to Eq. (23):

[TABLE]

As we will see in our explicit calculations, $\Delta{I}_{qq}^{(2,S)}(z)$ is determined by the $q^{\prime}\to q$ channel, $\Delta{I}_{q\bar{q}}^{(2,V)}(z)$ is determined in the $\bar{q}\to q$ channel, $\Delta{I}_{qq}^{(2,V)}(z)$ is determined in the $q\to q$ channel, and $\Delta{I}_{qg}^{(2)}(z)$ is determined in the $g\to q$ channel. For completeness, we reproduce the necessary scheme transformations of these quantities from the Appendix:

[TABLE]

We organize our calculation in terms of cut diagrams, which distinguish whether the two additional partons that appear in the NNLO calculation are virtual or real. The double-virtual corrections, in which both additional partons are virtual, are scaleless and vanish in dimensional regularization. This leaves us with real-real and real-virtual diagrams to calculate. Fig. 2 shows the symmetric diagrams contributing to the real-real corrections at NNLO as a representative example of the types of contributions which occur. Interference diagrams are not explicitly shown. Fig. 3 shows the real-virtual corrections at NNLO. Mirror diagrams are not explicitly shown.

V.1 The $q^{\prime}\to q$ and $\bar{q}\to q$ channels

We will begin by showing the details of our treatment of the $q^{\prime}\to q$ transition, since this calculation is sufficiently compact to demonstrate explicitly. In light-cone gauge there is only one diagram that contributes: Fig. 2(e), in which the quark entering the hard-scattering cross section is a different flavor than the one coming from the PDF. This diagram can be calculated by using the standard QCD Feynman rules:

[TABLE]

Here, $k_{1}$ and $k_{2}$ are the momenta of the intermediate particles that pass the cut. The final-state phase space for the NNLO real-real correction can be parameterized as

[TABLE]

It is straightforward to evaluate the trace appearing in Eq. (V.1) in the HV scheme. We are left with integrals containing the momenta $k_{1}$ and $k_{2}$ over the phase space of Eq. (56).

We facilitate our calculation using integration-by-parts identities (IBP) Tkachov:1981wb ; Chetyrkin:1981qh , implemented in the computer code LiteRed Lee:2013mka . The only non-standard aspect of our implementation of the IBP identities is our treatment of $\epsilon$-dimensional momenta. The HV scheme relations of Eq. (21) require us to separate the momenta $k_{1}$ and $k_{2}$ into 4-dimensional and $\epsilon$-dimensional pieces, $k_{\mu}=\tilde{k}_{\mu}+\hat{k}_{\mu}$. Upon doing so we obtain integrals that depend explicitly upon the $\epsilon$-d momenta $\hat{k}_{\mu}$. We introduce auxiliary vectors that parameterize the $\epsilon$-dimensional direction to handle such contributions. As an example, suppose our evaluation of the trace in Eq. (V.1) contains the dot product $\hat{k}_{1}\cdot\hat{k}_{2}$, and we wish to evaluate the corresponding integral $I^{d}[\hat{k}_{1}\cdot\hat{k}_{2}]$ that occurs upon integrating this expression over phase space. We note that $I^{d}$ can depend upon any other manifestly $d$-dimensional dot products in addition to its dependence on $\hat{k}_{1}\cdot\hat{k}_{2}$. A simple form-factor decomposition of this integral reveals that we can write

[TABLE]

where $v_{\perp}$ is a space-like vector with support only along the $\epsilon$-dimensional direction. The integral on the right-hand side is now written in a manifestly $d$-dimensional form, and can be handled using the standard IBP machinery. Similar relations can be derived for all structures appearing in our integrand. To obtain all the integrals needed in our calculation we must introduce two such auxiliary momenta, both with support only in the $\epsilon$-dimensional momenta but with an angular separation in this subspace.

After IBP reduction, $\Delta B_{qq^{\prime}}$ can be expressed in terms of four master integrals:

[TABLE]

The four needed master integrals can be derived through direct phase space integration:

[TABLE]

We note that there is no $z\to 1$ divergence in this channel, but there is a $t\to 0$ divergence that occurs upon combining these integrals. Dimensional analysis tells us that the result of the full calculation must be proportional to $t^{-1-2\epsilon}$. We expand this factor together with the renormalization scale $\mu$ in terms of plus distributions as follows:

[TABLE]

Upon performing this expansion we obtain the bare beam function $\Delta B^{bare(2)}_{qq^{\prime}}$. It is clear from Eq. (20) that the renormalized and bare beam functions are the same at this order for this channel. We can use Eq. (31) to express the beam function in terms of the matching coefficient:

[TABLE]

We note that the one-loop matching coefficient $\Delta\tilde{\mathcal{I}}_{qg}^{(1)}$ must be derived to higher orders in $\epsilon$ when used in this relation. All poles cancel in this expression for $\Delta\tilde{\mathcal{I}}_{qq^{\prime}}^{(2)}$, which is a powerful check of our calculation. We are left with the finite matching coefficient, which reads as follows after scheme transformation:

[TABLE]

Comparing to Eq. (33) we can confirm that the plus distributions in $t$ have the required coefficients, which serves as a second check of our result. We isolate the desired scale-independent matching coefficient:

[TABLE]

We further verify our identification of this term as $\Delta{I}_{qq^{\prime}}^{(2,S)}$ in our calculations of the $\bar{q}\to q$ and $q\to q$ channels. We note that all master integrals appearing in this channel are well-defined in dimensional regularization, without the need for an additional regulator. A nearly identical calculation for the $\bar{q}\to q$ channel gives the result

[TABLE]

As a final check of our result, we note that the master integrals we obtain for the polarized beam function are identical to those that appear for the unpolarized one. It is simple for us to re-evaluate the Feynman diagrams without any polarization in order to obtain the standard unpolarized matching coefficients. We agree exactly with the results of Ref. Gaunt:2014xga for these channels. We note that $\Delta{I}_{q\bar{q}}^{(2,V)}$ satisfies the helicity conservation relation of Eq. (25).

V.2 The $g\to q$ and $q\to q$ channels

We now discuss aspects of the calculation that first appear when we consider $\Delta{I}_{qg}^{(2)}$. This channel has both real-virtual corrections and a non-trivial renormalization. The real-real corrections require nine independent master integrals. We use the method of differential equations Kotikov:1990kg to calculate them, facilitated by the use of the canonical form Henn:2013pwa . This special form of the differential equation was found through the Magnus algorithm Argeri:2014qva . Since they otherwise present no new feature beyond those already presented in the $q^{\prime}\to q$ channel, we do not discuss them further.

After performing the IBP reductions for the real-virtual diagrams we find three master integrals which can be solved by direct integration:

[TABLE]

We have let $q$ be the loop momentum while $k$ is the momentum for the on-shell parton crossing the cut. In all three master integrals, the final-state phase space integration is trivial,

[TABLE]

This fixes the momentum components of $k$ as: $k^{+}=\frac{t}{\omega},~{}k^{-}=p^{-}-\omega=(1-z)p^{-},~{}k_{\perp}^{2}=k^{+}k^{-}=\frac{1-z}{z}t$. The rest of the integration can be done by using standard Feynman parameterization techniques.

We note that all of these master integrals are well-defined in dimensional regularization, as are all the real-real master integrals. No new master integral appears in the calculation of the $q\to q$ channel, meaning that our entire calculation is well-defined without the introduction of an additional regulator. In addition, each cut diagram is well-defined without an additional regulator, since each can be expressed as a linear combination of master integrals with rational coefficients depending on $\epsilon$ and $z$. The calculation of the unpolarized beam function of Ref. Gaunt:2014xga was stated to require a regulator beyond dimensional regularization in intermediate stages. Since it is completely straightforward for us to also apply our master integrals to the calculation of the unpolarized beam function, we have done so and found complete agreement with the matching coefficients of Ref. Gaunt:2014xga , while finding that no additional regulator is needed for that calculation either.

We now detail the matching and renormalization in the $g\to q$ channel, where ultraviolet divergences first appear. The first step is to take into account the renormalization of the strong coupling constant in the $\overline{\rm MS}$ scheme, for which we need the following renormalization constant:

[TABLE]

The two-loop contribution to the bare beam function after $\alpha_{s}$ renormalization is

[TABLE]

From the relation between the bare and renormalized beam functions we can derive the renormalized beam function in the $g\to q$ channel as

[TABLE]

where $Z_{q}^{(1)}$ and $\Delta B_{qg}^{(1)}(t,z,\mu)$ have already been determined in Section IV. We note that the one-loop beam function is needed to ${\cal O}(\epsilon^{2})$ in this equation. With the renormalization finished, we next consider the matching. For the $g\to q$ channel, it is

[TABLE]

where the polarized PDFs are defined in Eq. (III). We have verified that all the divergences vanish in the matching, as required. After scheme transformation, the logarithmic terms all agree with the constraints shown in Eq. (33), providing a check on our result. We present here the desired scale-independent coefficient of the $\delta(t)$ term:

[TABLE]

The only new feature present in the calculation of the $q\to q$ channel is the appearance of $z\to 1$ singularities. They are straightforwardly extracted by a standard expansion in plus distributions:

[TABLE]

We note that the needed renormalization factor turns out to be equivalent to the one found for the unpolarized case Ritzmann:2014mka . The desired scale-independent matching coefficient is

[TABLE]

where

[TABLE]

We have applied our techniques to also calculate the unpolarized matching coefficients, finding complete agreement with the results of Ref. Gaunt:2014xga . We note that $\Delta{I}_{qq}^{(2,V)}$ satisfies the helicity conservation relation of Eq. (25).

VI Conclusions

In this paper we have calculated the longitudinally-polarized quark beam function through NNLO in QCD perturbation theory. This is the last remaining quantity needed to apply the $N$-jettiness factorization theorem through the NNLO level to polarized collisions. We expect that our result for $\Delta B$ will be useful in numerous applications, from enabling NNLO results for polarized collisions at RHIC and future colliders, to the study of event shapes at a future EIC. We have presented in detail the technical features of our derivation, which we believe will be useful in future calculations of similar quantities. As a by-product of our calculation we have provided the first independent check of the unpolarized quark beam function obtained in Ref. Gaunt:2014xga . This quantity has been heavily used in the $N$-jettiness subtraction approach to fixed-order calculations, and an independent derivation of the result was highly desirable. We have shown that the NNLO beam function can be obtained using only dimensional regularization to handle all singularities. We look forward to phenomenological applications of our results.

Acknowledgements

We thank A. Vladimirov for questions regarding the factorization scheme dependence in an initial version of this manuscript, which led to corrections in our NNLO matching coefficients and an expanded discussion of this topic. We thank X. Liu and P. Mastrolia for helpful conversations. R. B. and U. S. are supported by the DOE contract DE-AC02-06CH11357. F. P. and H. X. are supported by the DOE grants DE-FG02-91ER40684 and DE-AC02-06CH11357. This research was also supported in part by the NSF under Grant No. NSF PHY11-25915 to the Kavli Institute of Theoretical Physics in Santa Barbara, which we thank for hospitality during the completion of this manuscript.

Appendix A Polarized splitting functions

We reproduce here the polarized splitting functions needed in our calculation for both before and after scheme transformation. These are consistent with the literature Vogelsang:1995vh ; Vogelsang:1996im . The leading-order results are

[TABLE]

The NLO results are

[TABLE]

with the scheme transformation relations shown below:

[TABLE]

The scheme transformation factors are defined in Appendix B, and the final physical splitting functions are

[TABLE]

with

[TABLE]

We note that contrary to previous definitions of these splitting functions Vogelsang:1995vh ; Vogelsang:1996im we have not included a factor of $2N_{F}$ in $\Delta P_{qg}^{(0)}$, $\Delta P_{qg}^{(1)}$ and $\Delta P_{qq}^{S(1)}$.

Appendix B Factorization scheme transformation

In our calculation, we consider the following scheme transformation matrix

[TABLE]

where $a_{s}=\alpha_{s}/(4\pi)$ and $z_{qq}=z_{qq}^{V}+z_{q\bar{q}}^{V}+z_{qq}^{S}$. In the above equation, we have dropped $z_{qg}^{(n)}$, $z_{gg}^{(n)}$ and $z_{gq}^{(n)}$, although $z_{gq}$ can be reintroduced if desired Moch:2014sna . The helicity-conservation requirements lead to the transformations

[TABLE]

We can derive

[TABLE]

by requiring

[TABLE]

From helicity conservation, we can also derive $z_{qq}^{(2,V)}$ and $z_{q\bar{q}}^{(2,V)}$ from the $q\to q$ and $\bar{q}\to q$ channels, respectively. The quantities we derive are in agreement with previous results in the literature for these quantities Ravindran:2003gi . We cannot derive $z_{qq}^{(2,S)}$ from our calculation. However, given the agreement of $z_{qq}^{(2,V)}$ and $z_{q\bar{q}}^{(2,V)}$ with Ref. Ravindran:2003gi , we take the result for $z_{qq}^{(2,S)}$ from this reference. For completeness we list all needed terms here:

[TABLE]