Massive quarks at one loop in the dipole picture of Deep Inelastic Scattering

We calculate the light cone wave functions for a virtual photon to split into quark-antiquark states, including for the first time quark masses at one loop accuracy. These wave functions can be used to calculate cross sections for several precision probes of perturbative gluon saturation at the Electron-Ion Collider. Using these wave functions we derive, for the first time, the dipole picture DIS cross sections at one loop for longitudinal and transverse virtual photons including quark masses. The quark masses are renormalized in the pole mass scheme, satisfying constraints from the requirement of Lorentz invariance of the quark Dirac and Pauli form factors.


INTRODUCTION
It is believed that in very high energy hadronic collisions, the partonic constituents of hadrons and nuclei exhibit a qualitatively new kind of gluon saturation behavior, characterized by strong nonlinear interactions even at short distance scales where the coupling is weak. An experimentally clean way to study this regime are high energy Deep Inelastic Scattering (DIS) experiments. Studying gluon saturation is a key science goal of the future Electron-Ion Collider (EIC) [1, 2], which will address it with a broad program of precision measurements. The EIC can reach further into the saturation regime than previous measurements at HERA, because it also collides heavy nuclei, where saturation phenomena are enhanced [3].
One could search for signals of gluon saturation in the renormalization group evolution of cross sections as functions of the kinematical variables Q 2 and x [4,5]. With the EIC collision energy, however, the kinematical leverarm to distinguish fine details or asymptotic features of evolution is limited, since evolution is only logarithmic in Q 2 or x. Instead, one must most likely look for evidence of saturation in a combination of high precision measurements of different processes. A key part in such a comprehensive set of studies is played by heavy quark observables, both inclusive and exclusive. Of particular interest are processes involving charm quarks, where the quark mass is heavy enough to justify a weak coupling treatment, but light enough to be sensitive to saturation effects. In a collinear factorization picture, the charm cross section is one of the most sensitive probes of smallx gluons at the EIC [6]. To access gluon saturation it is better to use instead the coordinate space dipole picture [7][8][9][10][11][12] of DIS, where the virtual photon emitted by the electron first splits to partonic constituents, which then eikonally interact with the target. The dipole picture naturally involves the eikonal scattering amplitudes, Wilson lines, used to quantify gluon saturation in the CGC picture [13][14][15]. In the dipole picture light quarks are affected by contributions of nonperturbatively large dipoles in the "aligned jet" configurations [16,17], but the mass of the heavy quark makes them safe from this problematic part of phase space.
The theoretical framework of choice to understand saturation and the dipole picture in high energy DIS is QCD light cone perturbation theory [18][19][20][21] (LCPT). Here one first calculates the photon light cone wave function (LCWF) describing the probability amplitude of the photon to split into a partonic state. The LCWF is a universal quantity in perturbative field theory. It is a necessary ingredient in cross section calculations for different inclusive and exclusive scattering processes [17,[22][23][24][25][26][27][28][29][30][31]. Recently, the photon LCWF has been calculated to one loop accuracy in QCD perturbation theory [32][33][34] leading to a description of the HERA inclusive cross section [35] with massless quarks (see also [36,37]). In this letter we report the result of the calculation of the so far unknown NLO γ * → qq wavefunction with massive quarks. This Letter will be accompanied by a longer paper [38] with full technical details on the calculation for the transverse photon, the longitudinal photon having already been presented in Ref [39] (see also [40]). In a separate follow-up paper we will discuss the issue of quark mass renormalization in LCPT in more detail.

CALCULATIONAL SETUP
Loop calculations in LCPT We use the Hamiltonian LCPT formulation of perturbative QCD [18][19][20][21]. This approach is an ideal one for high energy scattering, where particles move on lightlike trajectories and their dynamics is thus naturally quantized with a lightlike time coordinate. The additional advantage of the LCPT formulation in light cone gauge is that only physical degrees of freedom are present in the calculation, which comes at the expense of having additional "instantaneous" 4particle interactions resulting from the exchanges of nonpropagating degrees of freedom. An unfortunate disadvantage is that because of the separate treatment of longitudinal and transverse coordinates, the theory is not manifestly Lorentz-invariant at the quantum level.
In the LCPT approach, one develops the full quantum state of the incoming particle, in this case the virtual photon, in a Fock state expansion of bare states. At small-x, the partons interact with the color fields of the target, thus only Fock states consisting of quarks and gluons are of relevance here. The leading such component in the photon state is the quark-antiquark dipole, depicted in Fig. 1. At NLO one also needs to include corrections from gluon loops, and gluon emission diagrams, i.e. qqg Fock states.
The coefficients of the expansion of the interacting (photon) state in terms of bare states are known as light cone wave functions. Perturbatively they are obtained in terms of a set of diagrammatical rules [21,41]. For every vertex one includes a matrix element depending on the helicities, polarizations and momenta of the participating particles. Instantaneous vertices are denoted by vertical crossed lines (time propagates from left to right). For every intermediate state (including the final state), one includes a light cone energy denominator which, in a covariant perturbation theory language, originates from integrating over the light cone energy k − and setting it on shell using the pole of a propagator. One then integrates over loop momenta and sums over internal helicities.
The leading order γ * → qq wavefunction (see e.g. Refs. [7,9,10,22]) is obtained by evaluating the diagram of Fig. 1, with one gauge boson-fermion vertex and one energy denominator. The diagrams needed for the NLO calculation are the same as in the massless case [32][33][34], as are most other calculational details. For our calculation we need first the self-energy corrections for the fermions, Fig. 2. For transverse photons, there is also a self-energy correction from an instantaneous interaction, Fig. 3. Longitudinal photons are properly speaking parts of an instantaneous interaction with the lepton [42], and thus do not have instantaneous vertices. The photon-quark-antiquark vertex gets corrections from normal physical gluons, Fig. 4, and also from instantaneous interactions, Fig. 5. Gluon emission can happen via normal, Fig. 6, and, again only for transverse photons, instantaneous interactions, Fig. 7. We have evaluated all these diagrams.
The loop momenta are integrated over in 2 − 2ε transverse dimensions, with a cutoff α regularizing any soft divergences arising from longitudinal momentum integrals in the k + → 0 limit. After integrating over the loop momenta one sums over internal helicities and gluon polarizations. We have performed the helicity sums both in the conventional dimensional regularization (CDR) scheme as in Ref. [32,33], and in the four-dimensional helicity (FDH) scheme as in Ref. [34], with equal results for the cross sections.
Issue of mass renormalization At this order in perturbation theory also the quark mass is renormalized. We work in the Hamiltonian LCPT approach, where from the Lorentz-invariant Lagrangian one first derives a Hamiltonian, and only then sets out to canonically quantize the latter in light cone gauge A + = 0. In the Hamiltonian the fermion mass appears in two separate terms [43]. The free part of the Hamiltonian has a "kinetic mass", which determines the relation between light cone energy k − = (k 2 + m 2 )/(2k + ) and 3-momentum (k, k + ). There is also the "vertex mass," the coefficient of the light cone helicity flip term of the gauge boson-qq vertex (see Eq. (1)). The latter did not need to be renormalized in our earlier calculation [39] of the longitudinal photon state.
Lorentz-invariance at the original Lagrangian level guarantees that the kinetic and vertex masses are equal in nature. Regularization methods that break Lorentzinvariance, such as the transverse dimensional regularization combined with longitudinal cutoffs used in our previous calculations for massless quarks Refs. [32][33][34], require the restoration of this invariance at the loop level by separate renormalization conditions for the kinetic and vertex masses. This was already known from the pioneering LCPT calculations of Refs. [44][45][46][47]. One can, however, slightly modify the regularization procedure by including, in addition to the diagrams appearing here, also the "self-induced inertia" or "seagull" diagrams [21,48,49] before the integrations. In the latter case, it becomes possible to maintain the equality of the vertex and kinetic masses.
Spinor structure Our calculation is organized in terms of possible independent spinor structures of the wave function. For a transversally polarized photon at leading order, the spinor structure of the leading order light-cone gauge γ * (q) → q(k 0 )q(k 1 ) matrix element can be decomposed (see e.g. [39]) in terms of three independent spinor structures as where i, j are transverse indices and P = (k + 1 /q + )k 0 − (k + 0 /q + )k 1 is the qq relative transverse momentum. The result for the γ * → qq wavefunction after evaluating all the loop diagrams, Figs. 2, 3, 4 and 5, can be decomposed in terms of four structures where α s = g 2 /4π is the QCD coupling constant, C F = (N c 2 − 1)/(2N c ) and N c is the number of colors. We obtain the 4 separate scalar functions V T , N T , S T and M T by evaluating the loop diagrams. For the longitudinal photon, one can perform a similar, but simpler de-composition. Comparing (1) and (2) we can see that the vertex mass is related to M T . On shell renormalization scheme For mass renormalization in the on-shell scheme we must look at the wave function in a specific kinematical configuration that we refer to as the on-shell point, corresponding to a timelike virtual photon with q − = (q + /(2k + 0 k + 1 ))(P 2 + m 2 ) (note that the physical region for DIS is spacelike, q − < 0). From Lorentz-invariance we know that at the on-shell point the whole γ * qq vertex function can be expressed in terms of two known scalar functions, the Dirac and Pauli form factors It is a straightforward exercise to express F D (q 2 /m 2 ) and One mass renormalization condition is given by the requirement that the self-energy diagrams in Figs. 2 and 3 do not have a pole at the on-shell point, as discussed explicitly in Ref. [39]. For a Lorentz-invariant regularization including the self-induced inertia diagrams, no other conditions are needed and the four conditions for V T , N T , S T and M T at the on-shell point are additional nontrivial checks of our result. On the other hand, with the regularization scheme of Refs. [32][33][34], the condition on M T becomes an additional vertex mass renormalization condition and one is left with three consistency checks for V T , N T , S T . In both cases our result for the massrenormalized wave function is the same.
From wave function to cross section To calculate the inclusive DIS cross section, one additionally needs to specify the interaction of the state with the target proton or nucleus. In the CGC formalism [15] this is described by an eikonal interaction with a nonperturbatively strong color field. The field is parametrized in terms of Wilson lines as functions of the transverse coordinate. Thus one must, after performing the mass renormalization, transform the LCWF's into mixed transverse coordinate-longitudinal momentum space. The interactions of the mixed space states with the target bring in Wilson line correlators that are the same as in the massless case. Also similarly to the massless case, there are cancellations of divergences (appearing as 1/ε poles) and scheme-dependent terms between the qq and qqg contributions. In order to obtain a manifestly finite expression for the cross section these must be subtracted from the qqg terms and added to the qq terms. We are performing this step within the same subtraction scheme as in Ref. [34].

RESULT AND DISCUSSION
For high energy QCD calculations one needs the wavefunction in mixed transverse coordinate-longitudinal mo-mentum space. Some of the spinor matrix elements in Eq. (2) depend on the relative qq transverse momentum P. Thus, what is needed are the scalar functions V T , N T , S T and M T multiplied by specific powers of the transverse momentum and by the leading order energy denominator, Fourier-transformed to coordinate space. We denote this multiplication and transformation by F . The NLO γ * → qq LCWF, the main result of this paper, can be written as The Fourier transformed scalar function V T reads where we have defined Here +[z ↔ 1 − z] adds a term corresponding to the whole preceding expression with the replacement. The Fourier transformed scalar function N T reads where the functions Ω T N and I T N are given by: The Fourier transformed scalar function S T reads Finally, the Fourier transformed scalar function combination V T + M T − S T /2 reads where the function I T VMS is given by The above expressions use the function L, which is defined as where Li 2 is the standard dilogarithm function and the parameter γ = 1 + 4m 2 /Q 2 . To our knowledge this is a completely new fundamental result in perturbative QCD.
We have also calculated the total DIS cross section to one loop order, using our results for the qq LCWF and the more straightforward, but algebraically complicated gluon emission wave functions. After the cancellation of UV divergences between the qq and qqg contributions, the cross section has a similar structure as for massless quarks [33,34]: where α em is the QED coupling constant. The "dipole" contribution σ γ * L,T subt qq corresponds to just the quarkantiquark pair crossing the shockwave color field of the target. The qqg-term σ γ * L,T subt qqg corresponds to a quarkantiquark-gluon system crossing the shockwave. The integration in the limit k + → 0 for this gluon develops a logarithmically large contribution, which must be resummed into the B/JIMWLK evolution of the target color fields in the same way as for massless quarks [35,50]. The transverse coordinate and gluon momentum fraction integrals cannot be performed analytically in the general case, since they depend on the properties of the Wilson line correlators describing the target. The quark momentum fraction integrals are also best left for numerical evaluation, similarly as in the case of massless quarks.
In addition to integrals that are similar to the massless case, the mass-dependent parts include additional integrals over Schwinger parameters that we have not been able to perform analytically. These integrals are generalizations of Bessel K 0,1 -function integral representations appearing in the massless case. They are very well convergent, and we do not expect their numerical evaluation to be significantly more complicated than a numerical evaluation of a normal Bessel K 0,1 -function. All the explicit expressions of the cross sections are written out in the Supplemental material of this Letter.
In conclusion, after a lengthy calculation, we have obtained the one loop LCWF's for the process γ * → qq. These are new results in field theory by themselves, expressing the full one-loop structure of the photon-quarkantiquark vertex in light cone gauge. We believe our result will be an important element in many future calculations. For example, the LCWF's will enable several calculations of exclusive processes in high energy DIS, such as diffractive structure functions, diffractive dijets, and exclusive vector meson production, at NLO accuracy and including massive quarks. As a first important application, we have computed the full NLO cross section for DIS in the dipole picture with quark masses. The cross section expressions obtained in this work will pave the way for simultaneous global fits of total and heavy quark cross sections measured at HERA, following the massless quark case [35]. These cross sections will be crucial for obtaining more precise predictions for EIC cross sections including the effects of gluon saturation.

Supplemental material
Here we present the results for the transverse and longitudinal photon cross sections, where the latter was already published in [39] but is here shown in the notations of this paper. These require also the contribution of the radiative diagrams shown in Figs. 6 and 7.  Explicit expressions for transverse photon cross section Using the generic notation (15) the (subtracted) quark-antiquark (qq) contribution to the transverse photon total DIS cross section at NLO in α s with massive quarks reads where S 01 is the quark-antiquark dipole amplitude (see the definition in Ref. [39]), Q is the virtuality of the photon, z = k + 0 /q + is longitudinal momentum fraction, m is the quark mass and the short-hand notation for i = 0, 1 denotes a 2-dimensional transverse coordinate integral. The functions f T i which encode the one loop QCD corrections with massive quarks are given by the following expressions: Here the function L is defined in Eq. (14) and the functions Ω T V , I T V , Ω T N , I T N and I T VMS are given by Eqs. (6), (7), (9), (10) and (13), respectively.
The (subtracted) quark-antiquark-gluon (qqg) contribution to the transverse photon total DIS cross section at NLO in α s with massive quarks can be written as a sum of four contribution where the first and the second term in Eq. (18) are the UV subtracted parts and the last two terms are the UV finite contribution.
It is useful in the following to introduce the compact notation for coordinate differences with x nm = x n − x m . It is also useful to introduce the relative coordinates corresponding to diagrams (j) and (k) as x 3;(j) = x 0+2;1 , x 3;(k) = x 0;1+2 , x 2;(j) = x 20 and x 2;(k) = x 21 . Using this notation, the UV subtracted contributions can be written as and where S 012 is the quark-antiquark-gluon amplitude (see the definition in Ref. [39]). We have also introduced the generalized Bessel function integral with (x) being either (j) or (k) and and We note that the integral (22) could be seen as a generalization of the integral representation of Bessel functions that appear in the massless case [33,34]. These integrals are very rapidly converging at both small and large values of the integration variable, and hence they should be well suited for numerical evaluation as is.
The UV finite contributions can be written into the following form and Here we have also introduced the function H (x) defined as with (x) being (j) or (k).