Massive quarks in NLO dipole factorization for DIS: Transverse photon

We calculate the light cone wave functions for the QCD Fock components in a transverse virtual photon necessary for applications at next-to-leading order in the QCD coupling, including quark masses. We present a detailed calculation of both the one loop wave function for the quark-antiquark Fock component and the tree level wave function for the quark-antiquark-gluon Fock component. The quark masses are renormalized in the pole mass scheme, satisfying constraints from the requirement of Lorentz invariance. In particular the quark Pauli form factor at NLO is recovered from the on-shell limit of the quark-antiquark Fock component. We use our result to calculate the next-to-leading order correction to the high energy deep inelastic scattering (DIS) transverse structure function on a dense target in the dipole factorization framework. Together with our earlier result for longitudinal photons, this completes the calculation of the total deep inelastic scattering cross section in the dipole picture with massive quarks at next-to-leading order, enabling a comparison with experimental data.


I. INTRODUCTION
In the limit of high scattering energies, Quantum Chromodynamics (QCD) is believed to exhibit the phenomenon of gluon saturation. This means that partial wave scattering amplitudes can become of the order of unity, i.e. sensitive to unitarity requirements, even for processes at weak coupling (transverse) momentum scales. Understanding the behavior of QCD in this limit has been the object of much attention both experimentally and theoretically. In addition to high energy hadronic and nuclear collision experiments, this small-x regime of QCD can be probed in high energy deep inelastic scattering experiments, at the HERA collider and in a future Electron-Ion collider (EIC) [1,2].
A convenient theoretical tool to understand the behavior of QCD in this limit is provided by the Color Glass Condensate (CGC) effective theory description [3][4][5]. Here one starts from the fact that in the small-x limit the dominant degrees of freedom of the scattering are gluonic states in the Fock state of the high energy target hadron or nucleus. One then formulates the high energy scattering process as an eikonal scattering off a classical color field [6]. The advantage of this eikonal approximation is that one works in transverse coordinate space, where the unitarity of the scattering matrix is manifest. The relevant physical degrees of freedom describing the dense gluonic target are eikonal scattering amplitudes of a dilute probe in the target gluon field, light-like Wilson lines. The Wilson lines resum nonlinear interactions involving any number of gluons in the target, and are thus very well adapted for describing the nonlinear physics of gluon saturation. For the deep inelastic scattering (DIS) process at high energy the eikonal scattering approach leads to the dipole picture [7][8][9][10][11], where one factorizes the perturbative partonic structure of the virtual photon from the eikonal scattering of the partonic states from the possibly dense target.
The formalism of light cone perturbation theory (LCPT) [6,[55][56][57] provides a calculational and conceptional tool to develop the picture of high energy scattering in a systematical perturbative expansion. We use LCPT to formulate a weak coupling Fock-space expansion of the quantum state of the probe, expressed in terms of light cone wave functions (LCWF)'s . When Fourier transformed to transverse coordinate space, the LCWF's are naturally combined with the Wilson lines describing the target to obtain scattering cross sections.
In this work we address the DIS process with massive quarks at one loop order in QCD perturbation theory. This calculation is highly relevant both for its phenomenological applications and for the development of the LCPT framework. On the phenomenological side, consistently describing HERA results for F 2 and F c 2 simultaneously has been challenging in leading order fits with BK (or JIMWLK) evolution [58,59]. For understanding data from HERA and the EIC in terms of high energy QCD and it will be important to include F c 2 in the description at NLO accuracy. This paper, taken together with our earlier results for longitudinal photons [60], provides the full expressions needed to calculate the total DIS cross section including quark masses.
On the more theoretical level, there are not many LCPT calculations at higher orders in perturbation theory. While the LCPT formulation provides a more physical explicit picture of the scattering process, the technology for loop calculations has been less developed than in covariant theory. Both the Hamiltonian formalism and the light cone gauge break explicit Lorentz-invariance, which can make intermediate expressions more cumbersome. Also renormalization has been less well understood in LCPT than in covariant perturbation theory, both because of the small number of loop calculations, and the complications related to Lorentz-invariance and the gauge choice. In this paper we will be faced with the problem of quark mass renormalization in LCPT. To our knowledge, our calculation here is the first practical LCPT calculation of a physical observable at NLO where not just the divergent part of the mass counterterm is extracted [61][62][63][64], but one-loop mass renormalization is carried out in full in the pole mass scheme (see also discussion in [65]), and all the finite leftover NLO terms are calculated with full mass dependence.
This paper is a part of a series of papers on extending the CGC (or more generally dipole picture) calculations of small-x DIS at NLO accuracy to massive quarks using Light Cone Perturbation Theory. In our first paper [60] we calculated the photon wave function and the total DIS cross section for longitudinally polarized photons. The final result for the wave function, and the cross sections for both polarizations, are presented in an accompanying shorter paper [66], but all the details of the calculation for the transverse photon will be described here. Since much of the calculation is relatively similar to the previous calculations with massless quarks [43,44,46] and to the longitudinal polarization case [60], we will be relatively brief with the introduction in this paper and move straight to the point. We refer the reader to these earlier papers for more background on the physics and the calculational methods.
In addition to the longitudinal photon case being algebraically simpler than our present calculation, the calculation in Ref. [60] also did not encounter the full complexity of quark mass renormalization in LCPT in the same way as in this paper (although the issue is discussed in Ref. [60]). We will in this paper need to fully renormalize the quark mass in both the propagator correction (so called "kinetic mass") and vertex correction ("vertex mass"). However, we will treat the issue in a concise way in this paper, ensuring that we get the correct result that maintains Lorentzinvariance of the quark form factor at one loop. We will return to a more detailed discussion of mass renormalization, regularization and the related issue of the self-induced inertias [67] in a separate future paper.
The paper is structured as follows. We will first, in Sec. II, write down the expressions for the total NLO DIS cross section, defining our notations and normalization for the LCWF's. We will then briefly describe the calculation of the leading order cross section in Sec. III. At NLO, before moving to the individual diagrams, we will first discuss the overall spinor structure, kinematics and mass renormalization procedure of the wavefunctions in Sec. IV. We then calculate the loop diagrams in Sec. V and combine them to get the momentum space mass renormalized γ * → qq wavefunction in Sec. VI. This is then transformed to mixed longitudinal momentum and transverse coordinate space in Sec. VII. We then write down and Fourier-transform the tree-level γ * → qqg wavefunctions in Sec. VIII. In Sec. IX we input the wave functions into the expressions for the cross section and effectuate the cancellation of the remaining UV divergences between the contributions of the qq and qqg Fock states, to arrive at our result for the total DIS cross section which is summarized in Sec. X. We provide brief conclusions and an outlook for the future in Sec. XI. Several more technical parts of the calculation are presented in the Appendices.

II. DIPOLE FACTORIZATION FOR DIS: CROSS SECTION AT NLO
The DIS cross section can be expressed in terms of the cross section of a virtual photon scattering on the hadronic target. Following the discussion presented in [60], the total NLO cross section for a virtual transverse photon scattering from a classical gluon field takes the form Here, the quark (antiquark) helicity and color indices are denoted as h 0 (h 1 ) and α 0 (α 1 ). In the vertex the light cone spatial three-momentumq ≡ (q + , q) =k 0 +k 1 is conserved, i.e. q + = k + 0 + k + 1 and q = k 0 + k 1 .
where the color factor C F = (N c 2 − 1)/(2N c ) and N c is the number of colors. The phase space sums over the mixed space quark-antiquark (qq) and quark-antiquark-gluon (qqg) Fock states are given by: where x 0 is the transverse coordinate of the quark, x 1 that of the antiquark and x 2 that of the gluon. The reduced LCWF's ψ γ * T →qq and ψ γ * T →qqg (see the discussion in [60]) are independent of the photon transverse momentum q and cannot depend on the absolute transverse positions of the Fock state partons, just on their differences. Finally, the quark-antiquark (S 01 ) and quark-antiquark-gluon (S 012 ) amplitudes are defined as for the quark-antiquark (S 01 ) and quark-antiquark-gluon (S 012 ) amplitudes. Here the fundamental (F) and the adjoint (A) Wilson lines are defined as light-like path ordered exponentials for a classical gluon target where t a and T a are the generators of the fundamental and adjoint representations, respectively. What we will do in this paper is to calculate the LCWF ψ γ * T →qq to one loop order (up to α s ) and the gluon emission LCWF ψ γ * T →qqg to tree level, including quark masses. We will then insert the result into Eq. (1) and subtract and add a term to make explicit a cancellation of a transverse UV divergence between the two contributions.

III. LEADING ORDER WAVE FUNCTION
As in the longitudinal case [60], let us first write down the leading order LCWF contribution to the transverse photon splitting into a massive quark-antiquark dipole. Following the notation shown in Fig. 1, the leading order γ * T → qq LCWF can be written as to coordinate space, and for the longitudinal momentum in the phase space integration measure (2). The leading order light cone energy denominator ED LO simplifies to Here we have introduced the notations with the momentum fraction z = k + 0 /q + and 1 − z = k + 1 /q + with z ∈ [0, 1]. The interpretation of P is that it is the relative transverse momentum of the quark-antiquark pair on the light cone, where the plus-momentum plays the role of a mass in a 2-dimensional nonrelativistic system. After Fourier-transforming to transverse coordinate space, κ z will be the typical inverse size of the dipole.
The transverse photon splitting vertex into massive quark-antiquark (qq) pair is given by where the parameter e f is the fractional charge of quark flavor f , e is the QED coupling constant, and the compact notation for the spinorsū(0) =ū(k 0 , h 0 ) and v(1) = v(k 1 , h 1 ) has been introduced. The transverse momentum dependence in the massive spinors u and v and in the photon polarization vector can be extracted out by using the decomposition derived in [60]. For the leading order quark spinor structure in Eq. (10), this decomposition yields Indeed, in each term in the right-hand side of Eq. (11), the γ + Dirac matrix projects out the components that depend on the transverse momenta (and masses) from the spinors u and v, so that only the components depending only on the light cone momenta and on the light cone helicities (sometimes called good components) survive (see e.g. [43]). For the DIS cross section, we need to Fourier transform the momentum space expression of the LCWF, supplemented with the delta function for overall transverse momentum conservation (see [60]), into mixed space. The reduced wave function ψ γ * T →qq LO is extracted from the two-dimensional Fourier-transform of the wavefunction Ψ γ * T →qq LO by separating an overall kinematical factor related to the total transverse momentum, and a color factor, as: Thus the reduced LCWF in mixed space is given by the following expression Here we have defined a Fourier transform operator F that corresponds to the Fourier-transform of its argument, divided by the leading order energy denominator and an overall factor, as We have also introduced the compact notation x 01 = x 0 − x 1 for the difference of transverse coordinates x 0 and x 1 . The shorthand notation F will be very convenient in the NLO calculation. The evaluation of such Fourier transforms is discussed in Appendix C.
It is then straightforward to perform the remaining (D − 2)-dimensional Fourier integrals in Eq. (14). All in all, the final result for the reduced leading order LCWF in mixed space can be written in the following form where the function K ν (z) is the modified Bessel function of the second kind. Setting D = 4, it is easy to check that one recovers the conventional result for the wavefunction [7].

IV. STRUCTURE OF THE NLO CORRECTIONS
A. Spinor structure and mass renormalization Before proceeding to calculate these contributions, let us first discuss the general Lorentz and gauge invariance constraints for the wave functions, and its relation to mass renormalization. As discussed already in our previous paper [60], the quark mass appears in two separate terms in the QCD light front Hamiltonian. The mass in the free fermion term that determines the relation between light cone energy k − and light cone momentum k 2 + m 2 is referred to as the "kinetic mass". On the other hand, as can be explicitly seen from the decomposition (11), also a part of the interaction term of the Hamiltonian is linearly proportional to the mass, which we call the "vertex mass". Lorentzinvariance at the level of the Lagrangian guarantees that the two masses stay equal. However, the regularization procedure of transverse dimensional regularization and a longitudinal cutoff, which has proven convenient in loop calculations [43,44,46], breaks Lorentz-invariance. One is left with two options. One is to modify the regularization procedure to take this into account. It turns out that such a regularization procedure in fact can be found, by only a slight extension of the one that has been used so far, affecting only the treatment of the "self-induced inertia" or "seagull" diagrams that vanish in our earlier scheme [57,67,68]. We will discuss this option in more detail in a future work. The other option is to accept that the two masses will become different after loop contributions, and must be renormalized by separate renormalization conditions at each order in perturbation theory. We have verified explicitly that, at least at one loop, both procedures lead to exactly equivalent results.
What happens in a practical calculation is the following. As discussed in the case of the longitudinal polarization [60], the kinetic mass renormalization is determined by the "quark propagator correction diagrams". Here the renormalization condition in the pole mass scheme is determined by the requirement that the leading divergence resulting from the energy denominator approaching zero in the on-shell limit is absorbed into mass renormalization. This procedure, discussed explicitly in [60], is identical in the case of the transverse polarization, since the propagator corrections completely factorize from the photon vertex. If the propagator correction diagrams are evaluated with a Lorentz-invariance conserving regularization scheme, the ensuing counterterm can directly be used to renormalize the vertex mass. If, on the other hand, one uses the straightforward k + -cutoff of [43,44,46], one needs a separate renormalization condition for the vertex mass, to restore Lorentz-invariance. In the case of the longitudinal photon polarization, the leading order vertex does not have a light cone helicity flip term proportional to the mass. There is thus no vertex mass parameter to renormalize and, and vertex mass renormalization is not needed.
Since the difficult issues with mass renormalization are related to the interplay of the regularization method and Lorentz-invariance, it is useful to start by checking the constraints of Lorentz-invariance on the wavefunction. In the transverse photon case, the NLO correction to the initial-state LCWF for γ * T → qq must be proportional to the polarization vector of the photon. It can also only depend on one transverse momentum vector P, and involve matrix elements of different Dirac matrix structures between the quark and antiquark spinors, which can be related to each other by Dirac matrix commutation relations, and by using the Dirac equation. Based on these facts it can be seen that the LCWF can be written as a linear combination of four independent spinor structures. In fact, for different stages of the calculation it is convenient to use different bases for these spinor structures. We start by introducing a basis that is convenient for the loop calculations and first choose the independent structures as: In terms of these four structures we can write the wave function up to NLO in terms of four form factors V T , N T , S T and M T which are defined, extracting some constants for future convenience, by Here the form factor V T multiplies the leading-order-like structure, which has both light cone helicity nonflip contributions, and flip contributions proportional to the quark mass. Out of the new structures that only appear at one-loop order, N T is a nonflip contribution and S T and M T are light cone helicity flip contributions. In general, the parametrization of the γqq vertex function based on Lorentz and gauge invariances, relevant for the γ * T → qq amplitude, can be written as where F D and F P are the Dirac and Pauli form factors, respectively. As explained in Appendix F, one can show that the form factors V T , N T , S T and M T have to obey the following constraints: All of these constraints are evaluated at what we call the on-shell point Q 2 = (P 2 + m 2 )/(z(1 − z)). The on-shell point corresponds to the kinematical configuration of a time-like virtual photon decaying to a quark-antiquark pair with four-momentum conservation. For the DIS process the virtual photon is spacelike, and the on-shell point is outside of the physical region. We have, however, analytical expressions for the form factors, and they can easily be extended up to the on-shell point. The on-shell point is, also, the renormalization point in the on-shell, i.e. pole mass, quark mass renormalization scheme. Depending on the details of the regularization procedure, we can use these four constraints in two ways. If the regularization procedure preserves Lorentz-invariance, all four conditions serve as cross-checks of the result of the loop calculation. If, on the other hand, Lorentz-invariance needs to be restored, one of the constraint equations becomes a renormalization condition for the vertex mass. To see explicitly how this happens, let us see how a vertex mass counterterm appears in the wave function (17). For the purpose of mass renormalization, it is convenient to change from the basis in (17) to a different one, where the longitudinal component of the polarization vector is eliminated using ε + (q) = 0, q · ε(q) = 0. This leads to the following structure where clearly the first two lines are light cone helicity conserving, and the second two flip terms.
Following a bare perturbation theory approach, the mass appearing in the unrenormalized wave function (17) is in fact the bare mass m 0 . One then replaces m 0 = m − δm with δm ∼ α s , and inserts this relation into the equation. In a renormalized perturbation approach, on the other hand, one works all the time with the physical quark mass m, and inserts an additional 3-point vertex mass counterterm in the Hamiltonian. In both cases, it is obvious that the mass counterterm at order α s is associated with the last line of (23), which is the only one where the mass appears at LO, and gives a contribution of the form Since the form factor M T is the only one that is only associated with the same spinor structure, it is clear that the vertex mass renormalization will affect only the form factor M T . Thus we can separate M T into the loop contribution M T loop (P, z, Q 2 ), which depends on the kinematics, and the constant counterterm It is now clear how to use one of the form factor constraints, namely the last one (22), as a mass renormalization condition. One simply determines the counterterm contribution to M T by requiring (22) to be satisfied, i.e.
This also points to a consistency condition that the loop calculation must, and will, satisfy. Namely, when P 2 is set to the on-shell value P 2 = −Q 2 − m 2 the form factor M T loop (P, z, Q 2 ) must no longer depend on z or Q 2 . Also, if a Lorentz-invariance preserving regularization has been used, (26) must be automatically satisfied with the mass counterterm extracted from the kinetic mass renormalization condition. Thus we have established the effect of the vertex mass renormalization. Once the loop contribution M T loop (P, z, Q 2 ) is calculated, we will simply subtract from it its value at the on-shell point, corresponding to the vertex mass counterterm. The additional check of our result obtained from comparing the first three constraint equations (19), (20) and (21) to the well known result for the Pauli form factor at one loop in QCD (identical to the QED result [69,70] up to the replacement α s C F ↔ α em e 2 f ) is done in detail in Appendix F.
In practical calculations of loop diagrams, we will want to use a symmetry under exchanging the quark and the antiquark to restore some contributions, without calculating the corresponding diagrams explicitly. The relevant symmetry here is that the light cone wave function should stay invariant if, for a fixed photon light front helicity, one exchanges both the momenta and the helicities of the quark and the antiquark. In terms of the invariant momenta the quark and antiquark transverse momenta are given by k 0 = P + zq and k 1 = −P + (1 − z)q. Thus the exchange of three-momentak 0 ↔k 1 is achieved by the substitution z → 1 − z and P → −P. The spinor matrix elementsū(0)γ + v(1),ū(0)γ + [γ i , γ j ]v(1) andū(0)γ + γ i v(1) in Eq. (23) are independent of transverse momenta, and symmetric under the exchange z → 1 − z. Under exchanging the light cone helicities h 0 ↔ h 1 , the first light cone helicity conserving matrix elementū(0)γ + v(1) ∼ δ h 0 ,−h 1 and the light cone helicity flip oneū(0)γ + γ i v(1) ∼ δ h 0 ,h 1 are symmetric, and the second heliclity conserving oneū(0)γ + [γ i , γ j ]v(1) ∼ h 0 δ h 0 ,−h 1 antisymmetric. The scalar form factors only depend on the square of the transverse momentum P, but have a nontrivial dependence on the momentum fraction z. Requiring the invariance of the wave function under the simultaneous exchangek 0 ↔k 1 , h 0 ↔ h 1 we can then deduce the following symmetrization requirements for the scalar form factors: We will use these conditions to deduce the correct sign for the contributions of diagrams related to each other by the exchange of the quark and antiquark. Before moving on, let us note that for the purposes of Fourier-transforming the wavefunction to mixed space (which is done after mass renormalization) and for calculating the total DIS cross section, it is convenient to change the basis yet again. In fact, the Fourier-transform of both the S T contribution and the M T + V T one are more easily calculated if one splits the symmetric tensor P i P j into its trace and a traceless part. This additional reorganization changes the decomposition (23) into the form It will also turn out that the expression for the combination M T + V T − S T /2 will be simpler than the terms separately. The spinor matrix elements of this basisū(0)γ do not depend on the transverse momentum P and can be factorized from the Fourier-transform. The result of the Fourier-transform of this decomposition is the one already given in the accompanying shorter paper [66].

B. Diagrams, energy denominators and kinematics
The LCPT diagrams relevant for the calculation of the γ * T → qq LCWF at NLO with massive quarks are shown in Figs Fig. 6 there is an instantaneous gluon exchange (i) between the quark and the antiquark, in which the gluon momentum can either flow upwards into the quark or downwards into the antiquark. It is convenient to split up the contribution from this diagram into terms that are combined with one or the other of the diagrams in Fig. 4 according to the direction of the momentum flow. Due to the symmetry of the kinematics by exchange of the quark and the antiquark between the two classes of graphs, only the calculation of half of the diagrams is necessary; in this case we will calculate the ones labeled (a), (c), (e), (g) and the part of (i) where the momentum flows to the quark as in (e). The energy denominators appearing in the diagrams are the same as in the longitudinal photon case [60]. Following the notation introduced in Figs. 2 and 4, and imposing the plus and the transverse momentum conservation we find 3. Time ordered one-loop instantaneous self-energy diagrams (c) and (d) contributing to the transverse virtual photon LCWF at NLO. Diagram (c): Imposing plus and transverse momentum conservation at each vertex gives:q =k 0 +k 1 and k 0 +k =k 0 .
where L = −(k + 0 − k + )P/k + 0 . It is also convenient to introduce the change of variables (k, k + ) → (K, ξ) by parameterizing the transverse momentum integration in the loop by the relative momentum K of the gluon with respect to the quark after the loop as In these notations the energy denominators in Eqs. (32) and (33) can be cast into the following form where the momentum variable L and the coefficients ∆ 1 , ∆ 2 are given by: In the following subsections, we present the detailed computation of the NLO form factors corresponding to the one-loop self-energy diagrams, one-loop instantaneous diagrams and one-loop vertex diagrams.

V. CALCULATION OF THE LOOP DIAGRAMS
A. One-loop quark self-energy contribution In this section we compute the contribution to the γ * T → qq LCWF in Eq. (17) from the massive quark self-energy diagrams (a) and (b) shown in Fig. 2. Applying the diagrammatic LCPT rules formulated in momentum space yields the expression where the energy denominators are written down in Eqs. (6) and (32), and the spinor structure coming from the vertices (note that the summation is implicit over the internal helicities, gluon polarization and color) is given by the numerator By making the change of variables (k, k + ) → (K, ξ) (see Eq. (34)) and regulating the small k + → 0 (or ξ → 0) divergences by an explicit cutoff k + > αq + (or ξ > α/z) with the dimensionless parameter α > 0, we can simplify the expression in Eq. (37) to where the detailed calculation of the numerator in Eq. (38), performed in Appendix A, gives Similarly as in the longitudinal case, we note that from this expression one obtains both the FDH scheme result by taking the limit D s → 4, and the CDR one by setting D s = D.
Following the same steps as presented in [60], we find where the elementary scalar integral A 0 is defined in Ref. [43] (see Appendices D and E there): We now perform the renormalization of the kinetic mass precisely as in our earlier calculation for the longitudinal photon [60] (see Sec. V C). In the end, this amounts to a mass subtraction that makes the leading singularity at the on shell point P 2 + Q 2 + m 2 → 0, i.e. ∆ 2 1 → ξ 2 m 2 (see Eq. (36)) cancel. Note that the denominator P 2 + Q 2 + m 2 in the second term of Eq. (41) is precisely the kind of doubling of the energy denominator in Ψ γ * T →qq LO that results from a Taylor series development in powers of the mass counterterm. Thus the effect of mass renormalization is to replace A 0 (∆ 1 ) in the second term of Eq. (41) by A 0 (∆ 1 ) − A 0 (ξ 2 m 2 ). With some abuse of notation we will continue denoting the mass-renormalized wavefunction with the same notation Ψ . Together with the explicit expression [43] we can now evaluate the wavefunction as where the form factor V T (a) is given by The factor (D s − 4)/2(D − 4) is the regularization scheme dependent coefficient coming from the following integral and we have defined the function I V (a) as Note that this integral over ξ ∈ [0, 1] is finite. The mass renormalized LCWF for the quark self-energy diagram in diagram (b) shown in Fig. 2 can be now easily obtained by using the symmetry between the diagrams (a) and (b). As discussed in Sec. IV A, this implies that the contribution to the form factor V T from the antiquark propagator correction diagram (b) is obtained from the diagram (a) by the substitution z → 1 − z. Using this symmetry we get where from Eq. (45) one obtains with Summing the expressions in Eqs. (44) and (48), we obtain for the full contribution of the one-loop quark self-energy to the γ * where the NLO form factor V T (a)+(b) can be written as

B. One-loop instantaneous contribution
Let us then calculate the one-loop instantaneous diagrams shown in Figs. 3, 5 and 6. We start by writing down the expression for the diagram (c) in Fig. 3. Applying the diagrammatic LCPT rules formulated in momentum space yields the expression where the spinor structure coming from the vertices is given by the numerator After making the change of variables (k, k + ) → (K, ξ), we can simplify the expression in Eq. (53) to where the numerator Eq. (54) in these variables simplifies (see Appendix A) to Noting that the term proportional to K i vanishes upon integration over K, we then obtain Next, we compute the diagram (g) in Fig. 5. The expression in momentum space reads where the spinor structure coming from the vertices is given by the numerator In the (K, ξ) variables, we can simplify the expression in Eq. (58) to where the numerator in Eq. (59) (see Appendix A) simplifies to Again, the term proportional to K i vanishes, and we are left with the following expression Finally, we focus on the diagram (i) in Fig. 6. We calculate this in the kinematics k + 0 > k + 0 and denote this contribution as (i) 1 . The corresponding expression in momentum space reads where the energy denominator ED v is defined in Eq. (33) and the spinor structure coming from the vertices is given by the numerator By the change of variables (k 0 , k + 0 ) → (K, ξ) we can simplify the expression in Eq. (63) to where the numerator in Eq. (64) (see Appendix A) simplifies to Now, since the term proportional to (K + L) i vanishes by symmetry due to the integration in K, we obtain By comparing the expressions in Eqs. (57), (62) and (67) to Eq. (17), we observe that the instantaneous diagrams contribute only to the form factor M T as: and To simplify the expressions in Eqs. (68) and (69) we utilized the relation − (D−4) . The corresponding contributions from diagrams (d), (h) and (i) 2 can be obtained from Eqs. (68), (69) and (70) by symmetry, as z ↔ 1 − z.

C. One-loop vertex contribution
The only diagrams left to calculate for the full γ * T → qq LCWF at NLO are the non-trivial vertex correction diagrams (e) and (f) shown in Fig. 4. These two diagrams give contributions to all four form factors introduced in the decomposition Eq. (17). Since the contributions coming from the diagrams (e) and (f) are symmetric to each other by exchange of the quark and antiquark, we need to explicitly calculate only one of them, for example the diagram (e).
For diagram (e), the momentum space expression of the one-loop LCWF can be written as where the numerator (again the summation is implicit over the internal helicities, gluon polarization and color) is given by Applying again the change of variables k → K, we obtain Before starting to calculate this contribution, it is useful to make a small power counting argument for the transverse momentum integrals in relation to the quark masses and light cone helicities. One is performing a (D − 2)-dimensional integral with a total denominator ∼ K 4 at large K from the energy denominators. Therefore one can have UVdivergences in the integral if the numerator, coming from a product of three elementary quark-gauge boson vertices, is proportional to at least K 2 . From the decomposition of the elementary vertex (11) we see that light cone helicity conserving terms are proportional to the momentum K i , whereas light cone helicity flip terms are proportional to one power of the mass ∼ m. We will thus have two kinds of UV-divergent contributions, proportional to m 0 or to m 1 . The first ones behave just like the corresponding diagrams in the massless case. The second ones, proportional to the mass, will contribute to the renormalization of the quark (vertex) mass, which appears linearly in the leading order vertex. There will also be 4 kinds of finite contributions, proportional to m n , n = 0, 1, 2, 3. The number of powers of mass corresponds to the number of light cone helicity flips following the quark line. Thus n = 0 and n = 2 contributions are nonflip contributions for the overall diagram and can interfere with each other. Similarly n = 1 and n = 3 have an overall helicity flip between the outgoing quark and antiquark, and can interfere with each other. Only even powers of the mass up to m 4 will appear in the NLO cross section. Quark masses regulate the energy denominators, so there are no IR divergences in the transverse integrals.
From a technical point of view, the calculation of the numerator and the corresponding tensor integrals in Eq. (73) turns out to be a quite tedious task. Therefore, the details of the intermediate steps are documented in Appendix A 3, keeping the main text as readable as possible. Collecting the results from Appendix A 3, we find the following result where the form factors can be cast into the form: and In the above expressions the following compact notation has been introduced where the elementary integrals B 0 and B j are the ones defined in Ref. [43] (see Appendices D and E there). Diagram (f) is obtained from (e) by exchanging the quark and the antiquark. Thus, as discussed in IV A, we can obtain its contribution to the form factors with the substitution z → 1 − z and changing the sign for N T , so that the symmetry requirements Eqs. (27), (28), (29) and (30) are recovered.

VI. FULL MOMENTUM SPACE LCWF
A. Form factor V T The full form factor V T gets contributions from the self-energy diagrams (a) and (b) and from the non-trivial vertex correction diagrams (e) and (f). Consequently, by symmetry, one has where the mass renormalized expression for the V T (a) contribution is given in Eq. (45). To obtain the full V T contribution in momentum space, we shall next evaluate further the individual terms appearing in Eq. (75) for V T (e) . We follow the same strategy that was used in Ref. [60], and study the UV divergent and UV finite parts separately.

UV divergent part
To evaluate the UV divergent part in Eq. (75), we first remind that the UV divergent function A 0 (∆ 2 ) can be written as Before performing the integration over the ξ, we rewrite the coefficient ∆ 2 as where the zeroes in ξ are given by Substituting Eq. (82) into Eq. (81) and then performing the remaining ξ-integral analytically gives the following result where the functions I ξ;i for i = 1, 2, 3 are independent of transverse momenta P and are given explicitly in Appendix B.
Putting everything together, we get

UV finite parts
Let us then focus on computation of the UV finite part in Eq. (75), which turns out to be quite tricky. It is possible to directly compute the B 0 and B i integrals, but the result would be too complicated for further analytical integration, in particular for the required Fourier transform. Instead, we first Feynman parametrize the denominator appearing in B 0 and B i as where L = −(1 − ξ)P and the denominator can be rewritten as From Eqs. (86) and (87) we then find that a combination appearing in the definition of T + in Eq. (79) can be expressed as where the function I + is given by the simple expression The relation in Eq. (88) turns out to be very useful in what follows. The difficulty in this part of the calculation is that we want to perform the divergent integral 1 α/z dξ/ξ analytically. Being able to do this requires further manipulations of the somewhat complicated functions T ± in the integrand. However, the relations in Eqs. (86) -(89) make it possible to split the 1 α/z dξ/ξ integral into a divergent part that can be integrated analytically, and a remaining more complicated term that is finite in the limit α → 0. To this end, we first substitute the expression (88) into the definition of T + in Eq. (79) yielding the following result Furthermore, we note that one can rewrite the function T − from Eq. (79) as 1 Finally, combining the expressions in Eqs. (90) and (91), we arrive at the following expression With these intermediate results at hand, we obtain after the direct integration the term proportional to T − in the expression for V T (e) given in Eq.
where the functions J ξ;i with i = 1, 2, 3 which depend explicitly on the transverse momentum P are given in Appendix B. The point of these manipulations is that we have now managed to split the integral (93) that had a divergence in the momentum fraction regulated by α into, on the one hand, a part where both the momentum fraction and Feynman parameter integrals have been evaluated and, on the other hand, a part proportional to m 2 where the integrals remain, but the regulator has been taken to α = 0. The contribution (93) is, as it appears in the expression for V T (e) in Eq. (75), multiplied with 1/P 2 . However, as can both be seen in the general decomposition in Eq. (17) and deduced from the fact that the original expressions before transverse momentum integrals never have energy denominators ∼ P 2 , such a transverse IR-divergent behavior is completely spurious and must in fact cancel. Thus we know that in fact the expression (93) must in fact vanish in the limit where The "double" denominator DD P 2 →0 appears after combining the original expression ∼ 1/D and the subtraction term ∼ 1/D P 2 →0 into a common fraction. This expression manifestly vanishes at P 2 → 0 and thus cancels the factor 1/P 2 multiplying it in Eq. (75), as it should. As a separate point, we also note that the remaining ξ and x integrals in Eq. (94) are finite and could be done analytically. However, it turns out that they, unlike the parts that diverge as α → 0, need to be kept in an unintegrated form in order to Fourier transform our expressions into mixed space later. The form (94) is the one we then use to obtain the final result. Finally, for the T + term appearing in Eq. (75) we obtain after the direct integration Collecting all the UV divergent and UV finite terms together, we can cast the one-loop form factor V T (e) into the following form where we have defined the function Ω V (e) as and I V (e) is given by the following integral expression The double integrals appearing in Eq. (99) can be further simplified by performing a set of change of variables. Since we will use the same variable change in several parts of our calculation, let us introduce it here for a more general case. We consider the following double integral of the type appearing e.g. in Eq. (99): where the denominator D is defined in Eq. (87) and f is some smooth function depending on ξ and x. First, introducing the following change of variable where the denominator D now becomes Next, introducing the following change of variable In this expression the denominator is now linear in y. Finally, with the change of variable one obtains Interestingly, we will find that all the double integrals that are written as an our example integral (100), can be greatly simplified by performing the above chain of three changes of variables. Now we can further simplify our expression for I V (e) in Eq. (99) by applying the change of variables to y and χ. After some algebra, we obtain (108)

The full form factor V T in momentum space
At this stage, it is convenient to combine the results from the propagator correction diagram (a) from Eq. (45) and from the vertex correction (e) that we have just calculated in Eq. (97). This gives We note that the log 2 (α/z)-term and a part of the log(α/z) term from the vertex corrections canceled against the self-energy contribution, as in the massless case. Finally, we are ready to write down the full result for the form factor V T in Eq. (80). By reconstructing the contributions of diagrams (b) and (f) using the z ↔ 1 − z symmetry, as discussed in Sec. IV A, we obtain the full leading order-like form factor V T as where the functions Ω V and I V are defined as in terms of Ω V (e) from Eq. (98) and I V (a) , I V (e) from Eqs. (47) and (108) and the function L, which also appears in the longitudinal photon case [60], is given by It is now straightforward to check that in the massless quark limit (γ → 1) the functions defined above satisfy and correspondingly the expression in Eq. (110) reduces to the one obtained in Refs. [44,46].
The form factor N T gets contributions only from the two non-trivial vertex correction diagrams (e) and (f). Again, by using the symmetry discussed in Sec. IV A, we have The expression above is UV finite and free from the ξ → 0 singularities. Since the calculation follows the same steps as before, we can immediately write down the final result where the function Ω N is given by with and I N is given by the following integral expression It is straightforward to check that this contribution vanishes in the massless quark limit. This is an agreement with the findings in the massless case in Ref. [43]. The double integrals in Eq. (119) can simplified by performing the change of variables to (y, χ) introduced in Sec. VI A 3. This yields the following result C. Form factor S T The form factor S T gets contributions only from the two non-trivial vertex correction diagrams (e) and (f). Again, by using the symmetry discussed in Sec. IV A, we have The fully finite contribution S T (e) from Eq. (77) can be easily rewritten into the following form Collecting the expressions obtained in Eqs. (68), (69), (70) and (78) together, we find with T − from Eq. (92). In the on-shell scheme for mass renormalization, the vertex-mass counterterm M T c.t. is determined by the renormalization condition M T → 0 for P 2 → −(Q 2 + m 2 ). Consequently, the expression in Eq. (125) yields the following result We discuss the value of the mass counterterm a bit more in Appendix F 3, and continue here with our primary objective, the mass-renormalized LCWF. The full mass renormalized expression for the form factor M T in the general kinematics now reads Finally, we note that the first term can be further evaluated by using the following result and that the third term containing B 0 can be rewritten by using the Feynman parametrized form introduced in Eq. (86). Combining these final remarks, we obtain the following result with T − from Eq. (92).
We are now almost ready to pass on to the next stage of our calculation and Fourier-transform the LCWF to transverse coordinate space. As discussed in Sec. IV A, at this stage it is convenient to switch from our original basis of four scalar form factors V T , M T , S T and M T to one where the most rank-2 tensor, most complicated to transform, is traceless. This makes appear, as the coefficient of the "leading-order-mass-like" spinor structure, the combination V T + M T − S T /2, as seen in Eq. (31). It turns out that a significant simplification occurs in this particular linear combination. As the last step in evaluating the momentum space LCWF let us now write it down. This combination has quite a simple form where I VMS is given by the following integral expression We can further simplify the expression in Eq. (131) by first performing the change of variables to (y, χ). This yields

VII. COORDINATE SPACE LCWF
We now Fourier transform the full NLO result of γ * T → qq in Eq. (31) into mixed space. We first factor out the exponential dependence on the center-of-mass coordinate of the dipole and the momentum of the photon as described in [60]. This yields where the reduced LCWF ψ γ * T →qq up to NLO accuracy in coupling α s reads Here the leading order result ψ γ * T →qq LO is given by Eq. (13) and the NLO piece Fourier transformed to mixed space reads where the Fourier operator F is defined in Eq. (14). We remind the reader that in this basis, our expression for the form factor S T is written in a symmetric and traceless form, which turns out to be very convenient for the Fourier transformation into mixed space. In addition, as discussed in previous section, the combination M T + V T − S T /2 is quite straightforward to compute both in momentum and mixed space due to the simple analytical structure. Before computing the Fourier transforms in Eq. (135), we would like to clarify some important points. Firstly, all the UV finite terms appearing in this expression can be Fourier transformed in four dimensions and only the UV regularization dependent terms (including the D → 4 pole term) need to be Fourier transformed in D dimensions. Secondly, the terms which are independent of the transverse momenta P trivially factor out from the Fourier transform. Thirdly, the terms exhibit a number of different types of P dependence. All the corresponding transverse Fourier integrals needed for these terms are given in Appendix C.
Putting all these steps together, we first Fourier transform the momentum space expression of V T in Eq. (110). The final result can be cast into the following form where the function I V is given by To obtain Eq. (137), we have additionally performed the following change of variable y → u = y/(1 − y) in Eq. (108).
In addition, the above expressions use again the function κ v = v(1 − v)Q 2 + m 2 for the inverse transverse length scale associated with momentum fraction v ∈ {z, χ}. Next, we Fourier transform the momentum space expression of N T in Eq. (114). The final result can be cast into the following form where the function I N is given by Similarly, Fourier transforming the momentum space expression of S T in Eq. (121) yields Finally, the result for the Fourier transformed form factor combination V T + M T − S T /2 in Eq. (130) reads where the function I VMS is given by We have now obtained the full mass-renormalized one-loop mixed space LCWF for γ * → qq with massive quarks, the most important result of our paper. This result, without any details of it derivation, is also shown in the shorter companion paper [66].

A. Momentum space wavefunction
We then move to the tree level wave functions for gluon emission from a transverse photon state, which are needed for the full cross section at NLO. As in the massless case [44,46], we need to calculate four gluon emission diagrams, shown in Figs. 7 and 8.
Applying the diagrammatic LCPT rules and following the notations shown in Figs. 7 and 8, we obtain for the diagrams (j) and (l) in momentum space and where the energy denominators appearing in Eq. (143) and (144) can be written as and the coefficients introduced in the denominators are defined as Similarly, for the diagrams (k) and (m) we get and where the energy denominators appearing in Eqs. (148) and (149) can be written as and the coefficients introduced in the denominators are defined as Note that the energy denominators (146) and (151) are in fact equal, but it is convenient to write them in terms of different variables in the context of different diagrams. The transverse momentum dependence in the massive spinors u and v and in the gluon or photon polarization vectors can be extracted by using the decompositions where the remaining spinor matrix elements are independent of transverse momentum. Inserting the decomposition in Eqs. (153) and (154) into Eqs. (143) and (144), and noting that in the light-cone gauge ε / * σ (k 2 )γ + ε / λ (q) = −ε l λ ε * j σ γ + γ j γ l we find for the sum of the diagrams (j) + (l) the expression where the Dirac structure M jl (j) is defined as Similarly for the sum (k) + (m) inserting the decomposition in Eqs. (155) and (156) where the Dirac structure M jl (k) is defined as B. Fourier transform to coordinate space The Fourier transform to mixed space of the γ * T → qqg LCWF is given by the expression To simplify the Fourier transformation, we first make the change of variables (k 1 , k 2 ) → (P, K) in Eq. (157) with and correspondingly (k 0 , k 2 ) → (P, K) in Eq. (159) with Here we also introduce the following compact notation for the coordinate difference between a quark/antiquark p and the center of mass of the two-particle system m, n that it is recoiling against: with x nm = x n − x m . Using these we introduce the relative coordinates corresponding to the gluon emission diagrams (j) and (k). For the diagram (j) with emission from the quark, the natural coordinates are the separation between the gluon (k 2 ) and the quark (k 0 ) and the separation between the center of mass of the quark-gluon system and the antiquark (k 1 ): For the other diagram (k) with an emission from the antiquark, the natural coordinates that appear are the separation between the gluon and the antiquark and the separation between the quark and the center of mass of the antiquark-gluon system Next, performing the integration in Eq. (161) over the delta function of transverse momenta, we obtain the expression where the LCWF's in Eqs. (157) and (159) are written in terms of the new variables P and K and the relative coordinate separations. Adding the contributions in Eqs. (169) and (170) together, we obtain the following result for the full qqg LCWF with massive quarks in mixed space where we have factorized out the dependence on the center of mass coordinate of the partonic system 2 and the color structure from the reduced LCWF ψ γ * T →qqg . The reduced LCWF can be split into two parts: the part without light cone helicity flips (not proportional to m), and the part with one or two flips (proportional to m, m 2 ) as Here the light cone helicity nonflip contribution Σ lj is given by the following expression and correspondingly the contribution with one or two light cone helicity flips, proportional to powers of the mass, Σ lj m is given by It is straightforward to check that in the massless quark limit, the expression in Eq. (172) reduces to the result obtained in [44,46]. The (D − 2)-dimensional Fourier integrals of the type I,Î and J are defined and calculated in Appendix E, for which the following compact notation has been introduced: with the coefficients 2 See also the discussion in [60] the other momentum fraction ratios λ (x) from Eqs. (147), (152), and the coordinate differences defined in Eqs. (165), (166), (167) and (168). We now have the full expressions for the gluon emission wave functions in coordinate space.

IX. CALCULATING THE DIS CROSS SECTION
Let us now compute the transverse virtual photon cross section at NLO with massive quarks in the dipole factorization framework. As in the massless case [43,44,46], this consists of squaring separately the qq and qqg-contributions, and then effectuating a cancellation of an UV-divergent contribution between them. These steps happen in a similar way here, since the additional algebraic complexity for massive quarks resides more in the UV-finite pieces that are not affected by this subtraction. Thus we will be relatively brief here and refer the reader to the massless calculation for more details. We follow here the exponential subtraction procedure of [46].

A. Quark-antiquark contribution
Let us first write down the qq contribution to the DIS cross section at NLO in α s . Applying the formula for the cross section in Eq. (1), we find the following expression where the γ * T → qq LCWF squared is summed over the quark spins (h 0 , h 1 ), transverse photon polarization (λ) and the factor 1/(D s − 2) comes from averaging over the transverse photon polarization states.
We can now compute the LCWF squared by using the compact expressions for the LO and NLO wave functions in Eqs. (13) and (135), respectively. This leads to the following expression After performing the straightforward Dirac algebra, we obtain where the UV finite terms has been simplified by taking the D s = D = 4 limit.
Using the results obtained in Eqs. (136), (138) and (141), the qq contribution to the total cross section can be written as where the short-hand notation [x i ] = d D−2 x i /(2π) for i = 0, 1 denotes a (D − 2)-dimensional transverse coordinate integral. The function f UV contains the UV divergent term plus the regularization scheme dependent part, and it is given by Similarly, the UV finite functions f V , f N and f VMS are given by: Here we recall that Eq. (112), Ω N by Eqs. (116) and (117)

B. Quark-antiquark-gluon contribution
The qqg contribution to the DIS cross-section at NLO in α s is given by the second term in Eq. (1) where the γ * T → qqg LCWF squared is also summed over the emitted gluon polarization (σ). Following the straightforward but lengthy derivation shown in Appendix D, we find the result Here the functions K T qqg | UV and K T qqg | UV m are the UV divergent parts of the squares of the contributions without a light cone helicity flip and with one flip, respectively. They are given by the following expressions: and Note that a flip and nonflip contribution cannot interfere and thus there are no cross terms.
The functions K T qqg | F and K T qqg | F m are the remaining UV finite terms, without and with light cone helicity flips respectively, i.e. without or with explicit powers of the quark mass. The UV finite light cone helicity nonflip terms can be cast into the following form Similarly, the UV finite contributions with one or more light cone helicity flips are given by Note that there are contributions proportional to m 2 , coming from either two light cone helicity flips in the amplitude and none in the conjugate or one flip in both, as well as contributions ∼ m 4 with three flips in the amplitude and one in the conjugate. Finally, inserting the expression in Eq. (184) into Eq. (183), we obtain the gluon emission contribution to the cross section (189)

C. UV subtraction
Like in the longitudinal photon case [60], the UV renormalization of the coupling g is not relevant at the accuracy of the present calculation. Therefore, the remaining UV divergences have to cancel between the qq and qqg Fock states contributions on the cross section level. Due to the complicate analytical structure of the gluon emission contribution in Eq. (189), the UV divergent phase-space integrals cannot be performed analytically for arbitrary dimension D. Hence, it is desirable to understand the cancellation of UV divergences at the integrand level.
In the expression (189), the first two contributions inside the curly brackets contain two terms which are UV divergent when x 2 → x 0 and x 2 → x 1 , respectively. In addition, the remaining two contributions are UV finite and hence we have immediately taken the limit D s = D = 4 at the integrand level, which simplifies the further calculation considerably. In order to subtract the UV divergences, we will follow the same steps as presented in [60] and use the following property of Wilson lines at coincident transverse coordinate points which implies that lim Thus, the UV divergences in Eq. (189) are subtracted by replacing: and where the subtraction terms are given in terms of a two functions I ik UV and I i UV : Now, for the function I ik UV be a good UV approximation of the full integral, it must satisfy lim from which it follows that lim An identical analysis holds for the function I i UV . It is important to note that there is no unique choice for the UV divergent subtraction in Eq. (195). The only requirement for the subtraction is that the UV divergence needs to cancel between the qq and qqg contributions. Thus, it is sufficient for the subtraction to approximate the original integrals by any function that has the same value in the UV limits (for any D). Because of this cancellation, the integrals of the expressions inside the curly brackets in Eqs. (192) and (193) are finite, and one can safely take the limit D s = D = 4 under the x 2 integral.
In an arbitrary dimension D, the integrals I ik and I i (see Eqs. (E6) and (E7)) are given by It is straightforward to see that to get the leading behavior in the limit |r| 2 → 0 we can set λ = 0. This leads us to and where we have suppressed the dependence on the variable λ in the notation. Now there are several possible ways of performing the UV subtraction. Using the exponential subtraction procedure introduced in [46], we approximate the incomplete gamma function with where the exponential is independent of u allowing for an analytical calculation of the u-integral. This replacement has the correct behavior in the UV limit |r| 2 → 0, but also is regular in the IR limit of large |r| 2 → ∞. Here we note that another option would be to follow the polynomial subtraction scheme used in [44]. Here the subtraction function is polynomial in r. This however, introduces a new IR divergence, which must be compensated with another subtraction 3 . Finally, substituting Eq. (201) into Eqs. (199) and (200) gives and The following master integral will be used to compute the UV divergent r-integral:

D. UV subtracted results
For the UV subtraction terms coming from the diagram (j), we find and Similarly, for the diagram (k), we find and In order to simplify our subtraction terms further, we introduce the same parametrization as in the qq contribution: k + 0 = zq + , k + 1 = (1 − z)q + and k + 2,min = αq + and change of variables from (k + 1 , k + 0 ) → z. Hence, the sum of Eqs. (205) and (207) yields the following expression for the UV divergent light cone helicity nonflip subtraction contribution and for the light cone helicity flip contribution 4). We can now gather here the main result of our paper, which is the transverse virtual photon total cross section at NLO with massive quarks. The cross section can be written as a sum of two UV finite terms where the first term in Eq. (211) is the mass renormalized and the UV subtracted qq contribution, which is obtained by adding the UV subtraction terms in Eqs. (209) and (210) into Eq. (180). This gives where the short-hand notation x i = d 2 x i /(2π) for i = 0, 1 denotes a 2-dimensional transverse coordinate integral and we recall that  212), not proportional to α s , is explicitly the known leading order cross section for massive quarks and the functions f i encode the order α s NLO corrections. We have checked that in the limit of zero quark mass these expressions reduce to the known results in Refs. [43,46]. The second term in Eq. (211) is the UV subtracted qqg contribution. It contains both the UV-subtracted terms corresponding to the parts inside the curly brackets of Eqs. (192) and (193), and the parts of the qqg-contribution that were finite to begin with: where the last two terms correspond a contribution coming from the finite terms without a light cone helicity flip Eqs. (187) and the ones with one or more helicity flip terms, proportional to powers of the mass, (188). The final result for the qqg-part contains a significant number of terms. In order to clarify the structure we use the relative coordinates corresponding to the gluon emission diagrams (j) and (k): for the diagram (j) with emission from the quark the q(k 0 ) − g(k 2 ) relative separation x 2;(j) ≡ x 20 defined in Eq. (165) and the separation of the qg system to the recoiling antiquark (q(k 1 )) x 3;(j) ≡ x 0+2;1 from Eq. (166), and conversely for the diagram (k) with emission from the antiquark theq(k 1 ) − g(k 2 ) relative separation x 2;(k) ≡ x 21 defined in Eq. (167) and the separation of the recoiling quark to theqg system (q(k 0 )) x 3;(k) ≡ x 0;1+2 from Eq. (168).
In the case of the exponential subtraction scheme developed [46], i.e. using Eq. (201), the sum of two UV subtracted terms in can be simplified to and Here we have introduced the compact notation (216) and .
(217) The notations Q 2 (j) , Q 2 (k) , λ (j) , λ (k) are defined in Eqs. (147), (152) and ω (j) , ω (k) in Eq. (176). As discussed also in the longitudinal case [60], these integrals could be seen as generalizations of the integral representation of Bessel functions that appear in the massless case. In addition, 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.
Let us finally turn to the last two terms in Eq.
Similarly, the last term in Eq. (213) reads Here we have also introduced the functions H (j) and H (k) , which are defined as To summarize the abbreviated notations, we recall that the definitions needed for an evaluation of the qqg results can be found in • The momentum scales Q 2 (j) , Q 2 (k) and longitudinal momentum fraction ratios λ (j) , λ (k) in Eqs. (147), (152) and the momentum fraction ratios ω (j) , ω (k) in Eq. (176).

XI. CONCLUSION
In this paper, we have calculated, for the first time in the literature, the one-loop light cone wave function for a transverse photon splitting into a quark-antiquark pair including quark masses. After obtaining the one loop LCWF we also Fourier-transformed our result to mixed transverse coordinate-longitudinal momentum space. Combined with the tree level wave function for a quark-antiquark-gluon state this enabled us to obtain an explicit expression for the total transverse photon cross section in the dipole picture, suited for cross section calculations in the nonlinear gluon saturation regime.
We believe that our results are a significant in several ways. Firstly, we have provided the first calculation of the total DIS cross section for massive quarks in the dipole picture at NLO accuracy. This paves the way for a description of existing and future total DIS cross sections in the small-x saturation regime, following the recent first successful NLO fit for massless quarks in [47]. We provide in this paper the explicit expressions for the cross section which, together with a BK-factorization procedure, e.g. one of the three used in [47], can directly be used in similar fits. More generally, the light cone wave function is a central ingredient in any NLO calculation in the small-x dipole factorization formulation for DIS or photoproduction processes involving heavy quarks. Thus it is an essential ingredient in bringing the phenomenology of gluon saturation in high energy deep inelastic scattering to NLO accuracy. The light cone wave function is needed for calculations of, e.g. exclusive vector meson production, dijet production or diffractive structure functions in DIS. All of these processes can now be accessed at NLO accuracy in the saturation regime, also for massive quarks. On yet a more fundamental level, the light cone wavefunction, independently of the needs of a specific cross section calculation, is a fundamental universal quantity in perturbative QCD. It encodes the structure of a photon in light cone gauge, which is required for a proper partonic interpretation. In calculating the wave function to one loop accuracy with massive quarks, we have also had to address explicitly the issue of quark mass renormalization in light cone perturbation theory. This is an issue that had previously been discussed in the literature at the level of the divergent terms. To our knowledge our calculation is, however, the first LCPT loop calculation to perform a full quark mass renormalization in a specific scheme, including the finite terms. We will discuss the details of different regularization procedures and their relation to mass renormalization in more detail in a separate publication. However, already here we have demonstrated in practice the relation between mass renormalization in LCPT and Lorentz-invariance, encoded in the structure of the quark form factor. The final result for the wavefunction and cross section obtained in this paper are also summarized in a shorter companion paper [66].
1. Numerator for the quark self-energy diagram (a) The numerator for the quark self-energy diagram (a) can be written as where the four-vectors k 0 and k 0 are the same since the plus and the transverse momentum are conserved and both k 0 and k 0 are on-shell. The spinor structures inside the first two square brackets can be written in the helicity basis as:ū where the variable K is defined in Eq. (34). Using the completeness relation for the spinors and noting that γ + (k / 0 + m)γ + = 2(k + 0 − k + )γ + , we obtain Here the terms linear in the transverse integration variable K vanish due to rotational symmetry. Summing over the gluon helicity states yields where the definition of the leading order QED photon splitting vertex V γ * T →qq h 0 ;h 1 is given in Eq. (10). Finally, using the parametrization k + /k + 0 = ξ, we find the following expression 2. Numerator for the instantaneous diagrams (c), (g) and (i) 1 The numerator for the instantaneous diagram (c) can be written as To simplify this expression, we first use the completeness relation (see e.g. Appendix A in Ref. [43]) where n µ a µ = a + and Eq. (A4). This gives where the standard γ-matrix relations have been used, as well as the plus and transverse momentum conservation relationk 0 =k 0 −k. Note that this implies k / 0 γ + = (k / 0 − k /)γ + . Finally, using the Dirac equationū(0)k / 0 = mū(0), we get Note that in the shift k / − (k + /k + 0 )k / 0 the γ − term vanishes, leaving only the terms involving γ i . Thus, we can rewrite Using the variables K and ξ introduced in Eq. (34) the expression in Eq. (A11) can be rewritten as The numerator for the instantaneous diagram (g) can be written as Following the same steps as above we find Finally, the numerator for the instantaneous diagram (i) 1 can be written as The simplifications here again follow the same steps as above leading to 3. Numerator for the quark vertex correction diagram (e) The numerator for the diagram (e) can be written as With the kinematical variables as in Fig. 4(e), the independent spinor structures can be decomposed as where we have introduced the variables with L = ω 2 P and Substituting Eqs. (A19), (A20) and (A21) into (A18), and following the same steps as before, we find the expression where we have introduced the following compact notation and defined the coefficients a, b, e, f, g, h as and the coefficient a m , c m as To perform the contractions in Eq. (A24) requires tedious but straightforward algebraic manipulations. To this end, we have automatized the computation of the Lorentz contractions with FORM [71,72], and we have also cross-checked our results manually by pen and paper. The one-loop vertex correction diagram (e) can be written in the spinor basis of Eq. (17), and the result reads Here we have used the following compact notation γ = 1 + 4m 2 /Q 2 introduced in Eq. (83). Similarly, the functions J ξ;i which depend on the transverse momenta P are defined as follows: and (B6) Next we provide two important intermediate results used for the computation shown in Section V C. After integration over the ξ, we obtain the following result which turns the integral in Eq. (C1) into a simple Gaussian form in P. The resulting expression can be evaluated by using the standard Gaussian integral: where y = (y 1 , y 2 , . . . , y n ), b = (b 1 , b 2 , . . . , b n ) and a ij is a symmetric, non-singular and positive defined n × n matrix. These steps lead to the result The remaining one-dimensional integral can be performed by using the general formula Here K −β (x) is the modified Bessel function. Notice that for positive values of x K β (x) is an analytic function of β and that K β (x) is even in β, i.e. K β (x) = K −β (x). Hence, the integral in Eq. (C4) yields the result We also note that if the expression in Eq. (C1) contains a logarithmic function, one can first use the relation log(A) = lim α→0 ∂ α A α and then apply the Schwinger parametrization formula in Eq. (C2). It is then easy to show that for the rank-one tensor integral In addition, we need to compute the following symmetric and traceless rank-two tensor integral in D = 4 plus delta functions in x 01 , which we ignore. Furthermore, Fourier transforms of a product can be treated by a partial fraction decomposition before integrating: which leads to the Fourier-transform being a difference of two Bessel K functions with different arguments. Doing a double Schwinger parametrization of the denominators, which one would also be tempted to do, ends up expressing the result in the same form. Here the two terms come from the integration limits of the innermost Schwinger parameter integral after the usual shift, with then the outermost Schwinger parameter integration transforming the terms to Bessel K functions. One can also first combine the denominators with a Feynman parametrization and then perform the integral with a single Schwinger parameter. This gives a result as a Feynman parameter integral of a single Bessel K function. However, this Feynman parameter integral can be performed using the Bessel equation satisfied by the K functions. This again results in the same difference of two terms as the partial fraction method, coming from the upper and lower integration limits of the Feynman parameter integral. In some cases we have been able to express the differences of Bessel K functions as a single Bessel function, using a partial integration in some other Feynman parameter. However, this is only possible on a case-by-case basis, not in general.
Furthermore, for the terms which contain instantaneous diagrams, we obtain the following expressions: and Next, we concentrate on the contributions with one or two light cone helicity flips in the amplitude. We start by computing the second term |Σ lj m | 2 in Eq. (D1). Again, to simplify the algebra, we divide the function Σ lj m in Eq. (174) into the two parts where we have defined the functions A lj m and B lj m as: and Here the coefficients a i and b i for i = 1, 2, 3 are defined as , a 2 = (k + 2 ) 2 (k + 0 + k + 2 ) 2 q + , a 3 = This yields the following expression By taking the complex conjugate of the function A lj m , we find for the first term on the right-hand side of Eq. (D17) the following expression Let us then perform further simplifications for the terms ∝ a 2 1 from the above expression. This part is given by Performing the sum over h 1 with γ + (k / 1 − m)γ + = 2k + 1 γ + , and then using the relations we find that Eqs. (D19) can be simplified to Next, utilizing the fact that I i (j) ∝ x i 20 (see Appendix E) so that only the term symmetric in i, i in the above expression survives, and finally performing a simple trace h 0ū (0)γ + u(0) = 4k 0,µ g µ+ = 4k + 0 (D22) then gives for Eq. (D19) the following result (D19) = (D s − 2)a 2 1 8k + 0 k + 1 4k + 0 (k + 0 + k + 2 ) + (k + 2 ) 2 (D s − 2) |I i (j) | 2 . × +a 2v (1) (2k + 1 − q + )δ k l (D s ) + q + 2 [γ k , γ l ] γ j γ + u(0)(Î k (j) ) * .

(D24)
This term is fully UV finite, and therefore we can immediately set D s = D = 4 everywhere. After performing the spinor and gamma-matrix algebra, we obtain the following expression h 0 ,h 1 2 e[A lj m (B lj m ) * ] = 4m 2 8k + 0 k + 1 − a 1 b 1 (2k + 0 + k + 2 )(2k + 1 + k + 2 ) + (k + 2 ) 2 e[I i (j) (I i (k) ) * ] Finally, we concentrate on the third contribution in Eq. (D1), which is the cross term between the helicity nonflip and the helicity one and two flip contributions. It is straightforward to show that terms coming from the helicity nonflip and single helicity flip contributions yield a vanishing result, since they involve a trace over an odd number of gamma matrices. Therefore, we are left with the terms which involve only a product of a helicity nonflip contribution and a two helicity flip one. These terms give the following result K i e iP·b e iK·r P 2 + Q 2 + m 2 K 2 + ω P 2 + Q 2 + m 2 + λm 2 (E2) P i e iP·b e iK·r P 2 + Q 2 + m 2 K 2 + ω P 2 + Q 2 + m 2 + λm 2 (E3) The remaining u-integral can be performed for generic dimension D, yielding the following results

Appendix F: Form factors and Lorentz-invariance
In this appendix we discuss some of the constraints from the Lorentz-invariance on the momentum space scalar form factors S T , N T and V T , and the longitudinal photon V L from our earlier paper [60], at the on-shell point, as discussed in Sec. IV A.
Our starting point is the usual parametrization of the γqq vertex function (based on Lorentz and gauge invariance), Γ µ (q) =F D (q 2 /m 2 ) γ µ + F P (q 2 /m 2 ) q ν 2m iσ µν , which introduces the Dirac and Pauli form factors F D (q 2 /m 2 ) and F P (q 2 /m 2 ). This form should now be compared with the general parametrization of the spinor structure of the wavefunction in Eq. (17). The difference between vertex corrections as encoded in the form factors F D and F P and the LCWF is that in the former the "outgoing" qq is an asymptotic outgoing state in the scattering process, but in the latter it is an intermediate state that will scatter off the target shockwave. This means that in the form factors the quarks are on mass shell, which in our case means P 2 = −Q 2 − m 2 . Note that for a quark-antiquark pair in the final state, which is the situation here, this means that the virtual photon must be timelike, q 2 = −Q 2 > 4m 2 . Relatedly, the qq state in the LCWF has its own energy denominator ED LO . In the comparing the LCWF with the form factors this energy denominator (which becomes zero at the on-shell point) must be left out. Leaving out the additional color and electric charge factors and contracting with the physical polarization vector for the photon, the condition that should be satisfied is u(0)ε / λ (q)v(1)F D (q 2 /m 2 ) + q ν ε λ,µ (q) 2m iū(0)σ µν v(1)F P (q 2 /m 2 ) = ū(0)ε / λ (q)v(1) 1 + α s C F 2π V T + q + 2k + 0 k + 1 (P · ε λ )ū(0)γ + v(1) α s C F 2π N T + q + 2k + 0 k + 1 (P · ε λ ) P 2 P j mū(0)γ + γ j v(1) This same constraint with a longitudinal photon polarization (the Γ + component of the constraint) was already considered in our earlier paper [60]. For a transverse photon polarization one can use the relation ε λ,µ (q)q ν u(0) iσ µν v(1) = 2m u(0) ε / λ v(1) − 2m (k + 0 −k + 1 ) q + q + 2k + 0 k + 1 (P · ε λ ) u(0) γ + v(1) (F3) − 2P 2 q + 2k + 0 k + 1 (P · ε λ ) P 2 P j u(0) γ + γ j v(1) − P 2 +Q 2 +m 2 q + 2k + 0 k + 1 u(0) γ + ε / λ (q) v(1) (F4) to arrive at the constraints (19), (20), (21) and (22). In the remainder of this appendix, we will verify that our results for the form factors satisfy the constraints (19), (20) and show how to use (22) for quark mass renormalization. We will also check an additional relation that relates the transverse and longitudinal polarization states. It would also be interesting to check what happens with the Γ − component, which does not couple to photons in the light cone gauge, but we leave such a discussion for another day.
where one uses the symmetry of the χ integrand under χ ↔ 1 − χ. Thus the l.h.s of the condition Eq. (19) is which is indeed the one-loop Pauli form factor. We have now verified that our LCWF's satisfy the constraints for the scalar functions S T and N T , and used the condition on M T for mass renormalization. There remains the condition (21) for V T . This is, however, a a bit more problematic. Indeed, the Dirac form factor and V T involve quantities which have soft divergences and UV divergences (related to quark wave function renormalization, not mass). Hence, there should be some regularization scheme dependence, making it more difficult to use these relations to compare our results with earlier ones in the literature. Nevertheless, such soft and UV divergences have some degree of universality. This allows us to obtain a constraint between two finite quantities by taking the difference between the constraints for the longitudinal and transverse photons. By subtracting Eq. (H7) in Ref. [60] for the longitudinal photon from its analogue (21) for the transverse one, we find that Note that the quark self-energy diagrams cancel in the difference. Then, one has, using the notations of Ref. [60] We will be comparing to results from our earlier paper [60], so we must start from intermediate results that are in a slightly less integrated form. Using Eqs. (73) and (78) from [60] one finds