Massive quarks in NLO dipole factorization for DIS: Longitudinal photon

In this work, we will present the first complete calculation of the one-loop longitudinal photon-to-quark-antiquark light cone wave function, with massive quarks. The quark masses are renormalized in the pole mass scheme. The result is used to calculate the next-to-leading order correction to the high energy Deep Inelastic Scattering longitudinal structure function on a dense target in the dipole factorization framework. For massless quarks the next-to-leading order correction was already known to be sizeable, and our result makes it possible to evaluate it also for massive quarks.

state in terms of the bare Fock states. These coefficients are known as light cone wave functions (LCWF's). Recent work [16,17,19,30,54,55] has led to important technical advances in performing LCPT calculations at loop level. However, the introduction of quark masses introduces some additional issues that must be dealt with. The typical way of regularizing these recent LCPT loop calculations has been to use a cutoff for longitudinal and dimensional regularization for transverse momentum integrals. It has long been known [56][57][58][59] that using such a cutoff procedure causes a problem for the fermion 1 mass renormalization. At a fundamental level the problem is associated with the well-known fact that the regularization procedure should respect the symmetries of the underlying theory. In the case of LCPT, the fermion mass appears in two different places in the Hamiltonian that one quantizes. Firstly there is the free fermion term, where the "kinetic" mass determines the relation between the energy and momentum of a free fermion. There is also a quark mass in the quark-gluon interaction term, where the quark-gauge boson vertex consists of two parts. Out of these two the light-cone helicity conserving part is independent of the quark mass. The light-cone helicity flip term, on the other hand, is proportional to the "vertex mass" of the fermion. The equality of the kinetic and vertex masses is due to the rotational invariance of the underlying theory at the Lagrangian level. In LCPT, however, one first derives from the Lagrangian the Hamiltonian formulation of the theory and only then starts to quantize it. If this quantization uses a regularization procedure that breaks (3-dimensional) rotational invariance, it can happen that one needs to separately renormalize the two masses at each order in perturbation theory with an additional renormalization condition restoring rotational invariance [60]. Both this conceptual issue, and the more complicated structure of the basic quark-gluon vertex due to the light-cone helicity [61] flip-term, make the calculation of the DIS cross section for massive quarks more involved than the corresponding one for massless quarks [16,17,19]. The factorization of high energy divergences into BK evolution of the target, on the other hand, is rather orthogonal to these additional complications from the quark masses. Thus, the high energy factorization aspect of the calculation is quite similar to the massless case, and will only be discussed briefly in this work, although it of course needs to be treated carefully in order to eventually compare the calculations to experimental data. This is the first in a set of papers, where we will fully analyze the calculation of the DIS cross section in the dipole picture to NLO accuracy with massive quarks. In this first paper we will concentrate on the case of a longitudinal virtual photon. For the longitudinal polarization the numerator algebra is simpler, making the calculation slightly less lengthy. More importantly, only the relatively simple "propagator correction" diagrams lead to a renormalization of the quark (kinetic) mass. Therefore, the mass renormalization in a pole scheme can be performed in a relatively straightforward way, without encountering the intricacies discussed above. Transverse photons will be addressed in a separate paper. There the helicity and polarization algebra is somewhat more complicated, and also the renormalization of both the kinetic and vertex mass needs to be addressed. Separately from these calculations, we plan to address more formal aspects of mass renormalization in LCPT, its relation to the regularization procedure and the treatment of so called "self-induced inertia" or "seagull" diagrams [50,62,63] in more detail in separate future work.
The rest of the paper is structured as follows: first, in Section II, in order to keep the paper self-contained, we give some basic background notions of light cone perturbation theory and explain the regularization approach used in this paper, although relatively briefly since this is very similar to the previous calculations in Refs. [16,17,19]. In Section III we recall how the DIS cross section is, in the dipole picture, factorized into light cone wave functions encoding the Fock state of the virtual photon, and Wilson line operators describing the interactions of these states with the target. We then derive in Sec. IV the leading order virtual photon-to-quark-antiquark light cone wavefunction in D spacetime dimensions as a warmup. In Sec. V we calculate the (mass renormalized) one loop corrections to this wavefunction in momentum space, which are then transformed into mixed transverse coordinate-longitudinal momentum space in Sec. VI. The tree level gluon emission diagrams needed for the real corrections to the cross section are calculated in Sec. VII. We then derive in Sec. VIII the cross section from these light cone wave functions, including the subtractions needed to cancel divergences between the real and virtual contributions. Finally we summarize our result for the cross section in Sec. IX and briefly discuss future steps in Sec. X. Many technical details on the calculations are explained in the Appendices.

A. Light cone coordinates and conventions
This section will be rather brief, as our notations are a combination of the conventions used in Refs. [16,17,19]. We refer the reader there for a more thorough explanation. For an arbitrary Minkowskian four-vector x µ = (x 0 , x) with x = (x 1 , x 2 , x 3 ) we define the light cone coordinates as where x + is the light cone time along which the states are evolved, x − is the longitudinal coordinate, and x = (x 1 , x 2 ) is the transverse position with x 2 = |x| 2 = x · x. In this paper, the spatial (in the light cone sense) three-vectors are denoted byx ≡ (x + , x). The components x + and x − are related to the Minkowski coordinates by The inner product of two four-vectors is from which one sees that x ± = x ∓ . The canonical conjugate of the longitudinal coordinate x − is the longitudinal momentum p + , and the evolution in light cone time x + is generated by the light cone energy p − . With the form of the metric in Eq. (3), the on-shell light cone energy becomes

B. Regularization
In loop calculations one needs to integrate over internal (on-shell) momenta. Here (contrary to e.g. the review [50]) we use conventional relativistic field theory conventions for the normalization. Thus, the momentum space integral measure is given by In the evaluation of loop or final state phase space integrals over longitudinal and transverse momenta, we encounter divergences which have to be properly regularized. In this work, following the same procedure as the one described in Refs. [16,17,19], we regularize the ultraviolet (UV) divergent (k → ∞) transverse momentum integrals via dimensional regularization by evaluating them in (D − 2) transverse dimensions and regularize (if needed) the soft divergences (k + → 0) with a cutoff. We consider two variants of dimensional regularization and present our results in both schemes, checking explicitly that the final result for the cross section is scheme-independent. These variants are the conventional dimensional regularization (CDR) scheme [64] that was used in [16,17], and the four dimensional helicity (FDH) scheme [65,66] that was used in Ref. [19]. The precise implementation of these schemes is carefully explained in [19]. Thus, in here, we only give an small overview of these two approaches. Both regularization schemes involve continuing space-time from four to D dimensions, but differ in the way how they treat the momenta and the polarization vectors of unobserved and observed particles. Here, the unobserved particles are either virtual ones which circulate in internal loops or particles which are external but soft or collinear with other external particles. All the rest are observed particles. The following rules listed below will be used to compute the LCWF's in D dimensions.
• In the CDR scheme, the momenta and polarization vectors of the observed and the unobserved particles are continued to D spacetime dimensions.
• In the FDH scheme, the momenta and polarization vectors of observed particles are kept in four dimensions (i.e. observed gluons have 2 helicity states) and the momenta of unobserved particles are continued to D > 4. The helicities (spinor structures) of unobserved internal states are treated as D s -dimensional, where D s > D at all intermediate steps in the calculation. Once all the helicity (Dirac) and Lorentz algebra is done, one sets D s = 4 and analytically continues to D < 4 dimensions.
In order to perform intermediate computations for both schemes at the same time, we will use the following rules: • Any factor of spacetime dimension arising from the Dirac and Lorentz algebra for spin or polarization vectors should be labelled as D s and should be distinct from the dimension of the momentum vectors D.
• The vertices are proportional to (D s − 2)-dimensional gluon polarization vectors ε * i λ , and summing over the helicity states of gluons yields where by δ ij (D s ) we denote a Kronecker delta in (D s − 2)-dimensional transverse space.
• The tensoral structures resulting from transverse momentum integrals are kept in (D − 2) dimensions. For example, if the integrand in the transverse momentum integral is proportional to k i k j , then the value of the dimensionally regulated integral is proportional to a (D − 2)-dimensional Kronecker delta δ ij (D) .
• Since D s > D, we have • Since both D > 4 and D s > 4 when the algebra is done, contractions of Kronecker deltas with fixed momentum p or polarization vectors ε i λ of external particles preserve these vectors • Only after the spin and tensor algebra is done, one can analytically continue the obtained result to D < 4 and take the limit D s → 4 in the FDH scheme or D s → D in the CDR scheme. If the calculation was done only in one of the two schemes, we would not need the notation D s at any intermediate step; this could be replaced by 4 or D. Here we will, however, present the results in both schemes simultaneously, and for this it is necessary to keep D s general.

III. DIPOLE FACTORIZATION FOR DIS: CROSS SECTION AT NLO
We use here the standard procedure where the DIS cross section is expressed in terms of the cross section of a virtual photon scattering on the hadronic target. The virtual photon cross section can be obtained by the optical theorem as twice the real part of the forward elastic scattering amplitude where the forward elastic amplitude M fwd γ * λ →γ * λ is defined in light cone quantization from the S-matrix element as [2] At high energy (or, equivalently, at small Bjorken x) the interactions with the target described by the operatorŜ E are eikonal interactions with a classical color field. These interactions are represented by Wilson lines, which are the scattering matrix elements of bare partonic states in mixed space: transverse position -longitudinal momentum. Thus in order to calculate the amplitude we need to develop the incoming virtual photon state in terms of these bare states.
The state |γ * λ (q, Q 2 ) i is the physical one-photon-state in the interaction picture. We start by its perturbative Fock state decomposition in momentum space, which is given by whereq = (q + , q), Q and λ are the spatial three momentum, virtuality and polarization of the virtual photon, respectively. The explicit phase space sum over the quark-antiquark (qq) and the quark-antiquark-gluon (qqg) Fock states are defined in D dimensions as: From now on we will leave the helicity h 0 , h 1 , polarization σ, λ and color α 1 , α 1 , a indices as well as the quark flavor implicit, summed over when appropriate. While the leading order cross section is of order α em , for NLO accuracy in this context we want to calculate the cross section to order 2 α em α s . This requires the knowledge of the quarkantiquark wave function Ψ γ * λ →qq to one loop order, and that of the quark-antiquark-gluon one Ψ γ * λ →qqg at tree level. The non-QCD Fock states containing electromagnetic interactions via photons or leptons, higher order Fock states (represented by the dots) and the photon wave function renormalization constant Z γ * = 1 + O(α em ) can be ignored since they do not contribute at the order considered in the present calculation.
In the dipole factorization picture with eikonal scattering, we need to switch from the full momentum space representation to mixed space, in which the kinematics of a particle is described by its light cone longitudinal momentum and transverse position. In this case, the Fock state expansion in Eq. (11) reduces to the following form in terms of the mixed space states defined as transverse Fourier transforms and so forth. The phase space sums over the mixed space qq and qqg Fock states are: The Fourier transforms to the mixed space wavefunctions Ψ γ * λ →qq and Ψ γ * λ →qqg are defined as: It is convenient to factorize out the dependence on the "center of mass" 3 coordinate and the color structure of the partonic system out from the LCWFs as: 2 Here α em = e 2 /(4π) and α s = g 2 /(4π) are the QED and QCD coupling constants, respectively. 3 In LCPT the "mass" that serves as a weight in this linear combination is the longitudinal momentum.
Time ordered light cone diagram (momenta flows from left to right) contributing to the longitudinal virtual photon wave function at leading order. Here, the quark(antiquark) helicity and color index are denoted as h 0 (h 1 ) and α 0 (α 1 ), respectively. In the vertex, the momentum is conservedq =k 0 +k 1 , where the spatial three momentum vectors are denoted asq = (q + , q), The reduced LCWFs ψ γ * λ →qq and ψ γ * λ →qqg 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. The color structures δ αβ and t a αβ are the only invariant SU(N c ) color tensors available for the qq and the qqg Fock states, respectively.
Using these notations, the total NLO cross section for a virtual photon scattering from a classical gluon field takes the following general form where we have introduced the notation 4 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.

IV. LEADING ORDER LONGITUDINAL PHOTON WAVE FUNCTION
We now recall the well known leading order light cone wave function for the longitudinal virtual photon splitting into a quark antiquark dipole. The labeling of the kinematics for this process is shown in Fig. 1. Using the light cone perturbation theory (LCPT) rules as presented in [16,19], the leading order γ * L → qq light cone wave function can be written as Here, the function 8 is an effective QED photon splitting vertex into a qq pair. In the physical DIS process with a longitudinal photon this is strictly speaking a part of an instantaneous interaction vertex with the lepton. However, as discussed in Refs. [14,16], it can in practice be treated as a separate vertex. As a remnant of this nature as a part of an instantaneous interaction, the longitudinal photon does not couple to instantanous quarks. Following Ref. [16], we use a condensed notation u(k 0 , h 0 ) =ū(0), v(k 1 , h 1 ) = v(1) etc. to shorten the expressions. In Eq. (23) the parameter e f is the fractional charge of quark flavor f and e is the QED coupling constant. In the LO diagram, there is only one intermediate state, which is the qq state 5 . The corresponding light cone energy denominator ED LO appearing in Eq. (22) is given by the differences of the light cone energies of the states Following the discussion of [16], it is possible to generalize the notation of LCWF to the case of an off-shell (virtual) photon by assigning the off-shell value to the light cone energy of the photon. The virtuality of the photon in LCPT actually corresponds to the light cone energy difference between the incoming electron and the outgoing electron+photon state. The off-shell value of q − is then used in each light cone energy denominator appearing in the perturbative expansion of the LCWF. The quarks energies are given by the mass shell relation where m is the quark mass. Thus, using Eqs. (25) and (26), the energy denominator given in Eq. (24) can be written as At this stage, it is convenient to introduce the relative transverse momentum P and the normalized photon virtuality squared Q 2 as where, in the second line, the variables P and Q 2 are expressed in terms of the longitudinal momentum fraction 1]. Using these notations, the LO energy denominator becomes where we have dropped the factor iδ since we only consider the case Q 2 > 0 in which the energy denominator is strictly negative. For the DIS cross section (given in Eq. (19)), we need to Fourier transform momentum space expression of LCWF into mixed space. Using the leading order expression in Eq. (22), we find for the reduced LCWF (see Eq. (18)) in mixed space the following expression FIG. 2. Time ordered (momenta flows from left to right) one-gluon-loop quark self-energy diagrams (a) and (b) contributing to the longitudinal virtual photon wave function at NLO. In diagram (a) imposing the spatial three momentum conservation at each vertex gives:q =k 0 ′′ +k 1 ,k 0 ′′ =k 0 ′ +k andk 0 ′ +k =k 0 .
Here we have introduced the notation x 01 = x 0 − x 1 . Performing the remaining (D − 2)-dimensional transverse integral by using the result in Eq. (D8) yields where the function K ν (z) is the modified Bessel function of the second kind. Setting D = 4, and calculating explicitly (see e.g. [67]) the matrix elementū(0)γ + v(1) one recovers the conventional result for the wavefunction [3].

A. Spinor structures and energy denominators
In the longitudinal photon case, at NLO accuracy in QCD, one finds that the initial-state LCWF for γ * L → qq can be written as a linear combination of two spinor structures. Using a convenient choice of basis for that space of spinor structures, one can write where the NLO form factors V L and S L are the light cone helicity non-flip (h.nf.) and helicity flip (h.f.) contributions, respectively. Note that while the transverse photon wave function has both light cone helicity flip and nonflip terms already at LO, for the longitudinal photon the flip term only appears at NLO. This term is related to the quark Pauli form factor, or the quark anomalous magnetic moment, and is discussed in more detail in Appendix H. The LCPT diagrams relevant for the calculation of the γ * L → qq LCWF at NLO with massive quarks are shown in Figs. 2, 3 and 4. There are two "propagator correction" diagrams, (a) and (b) in Fig. 2, with a gluon loop attached to the quark or the antiquark. Then, in Fig. 3, there are two "vertex correction" diagrams (c) and (d), corresponding to two different kinematical possibilities, with longitudinal momentum (which is always positive) flowing either up from the antiquark to the quark or vice versa. Finally, in Fig. 4, there is an instantaneous gluon exchange (e) 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 contributing to one or the other of the diagrams in Fig. 3 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) and the part of (e) where the momentum flows to the quark as in (c). Note that since the longitudinal photon is really fundamentally a part of an instantaneous interaction, there is no diagram with an instantaneous quark line.
In order to set the stage for our NLO (one-loop) computation, we start by writing down all the energy denominators appearing in the first class diagrams. In diagram (a), there are two energy denominators ED LO and ED a . Following the notation in Fig. 2(a) and imposing the plus and transverse momentum conservation, it is straightforward to show FIG. 3. Time ordered (momenta flows from left to right) one-gluon-loop vertex diagrams (c) and (d) contributing to the longitudinal virtual photon wave function at NLO. In diagram (c) imposing the spatial three momentum conservation at each vertex gives:q =k 0 ′ +k 1 ′ ,k 0 ′ +k =k 0 andk 1 ′ =k +k 1 .
FIG. 4. Time ordered (momenta flows from left to right) one-gluon-loop instantaneous diagram (e) contributing to the the longitudinal virtual photon wave function at NLO. Diagram (e): Imposing the spatial three momentum conservation gives: q =k 0 ′ +k 1 ′ . For the case, in which the gluon momentum flows to the quark, the spatial three momentum conservation also gives: and the LO energy denominator ED LO is given by Eq. (27). In diagram (c), there are three distinct energy denominators ED v , ED a and ED LO . Again, following the notation in Fig. 3(c) and imposing the momentum conservation, we obtain where the notation L = −(k + 0 − k + )P/k + 0 has been introduced. Finally, for the instantaneous diagram (e), in which the gluon momentum flows to the quark, there are two distinct energy denominators ED v and ED LO .
For all three diagrams, it is convenient to change variables from the momentum of the gluon k and k + to the relative transverse momentum K and the longitudinal momentum fraction ξ of the gluon with respect to the final state quark. These are the natural variables of the gluon emission and absorption vertices in diagram (a) and the gluon absorption vertex in diagram (c). They are defined as In these notations, the energy denominators in Eqs. (33) and (34) can be cast into the following form where the coefficients ∆ 1 and ∆ 2 are given by In the following subsections, we present the detailed computation of the NLO form factors V L and S L coming from the one-loop self-energy, vertex and instantaneous diagrams.

B. One-loop quark self-energy
We start by computing the contribution to the γ * L → qq LCWF in Eq. (32) from the one-loop massive one-loop quark self-energy diagrams shown in Fig. 2(a) and (b). Applying the diagrammatic LCPT rules formulated in momentum space yields the following expression for the diagram Fig. 2 where the energy denominators are written down in Eqs. (27) and (33). The product of light cone vertices (the summation is implicit over the internal helicities, gluon polarization and color) is given by the numerator We make the change of variables (k, k + ) → (K, ξ) as defined in Eq. (35) and regulate the small k + → 0 (or ξ → 0) divergences by an explicit cutoff k + > αq + (or ξ > α/z) with the dimensionless parameter α > 0. Using Eq. (36) together with Eq. (37), we can simplify the expression in Eq. (38) to The detailed calculation of the numerator in Eq. (39) is performed in Appendix B and gives We remind the reader 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.
At this point, it is convenient to define the UV divergent one-loop scalar integral as where µ 2 is the scale introduced by transverse dimensional regularization. Note that in the framework of transverse dimensional regularization Under these simplifications Eq. (40) reduces to where the result for the integral A 0 (∆ 1 ) can be written as [16] A 0 ( We note that this result is the same for any variant of dimensional regularization.

C. Quark mass renormalization
From the expression (45) for the integral A 0 (∆ 1 ), it is clear that the one-loop quark self energy diagram (a) is UV divergent. UV divergences have been already found for that diagram and for the other ones in the massless quark case [16,17,19]. In that case, the UV divergences cancel each other at the cross section level. UV divergences with such a behavior are expected to occur as well in the massive quark case studied in the present paper. However, since the expression for the longitudinal photon wave function (and thus the one for the longitudinal photon cross section) involves the quark mass already at LO, NLO corrections are expected to also include UV divergences associated with quark mass renormalization. Hence, one expect two types of UV divergences, which need to be disentangled.
In the LO expression (22) for the longitudinal photon wave function, the quark mass appears in the energy denominator ED LO , Eq. (29), but not in the numerator (23). In the bare perturbation theory approach, one uses the bare mass m 0 when writing the expression of the diagrams from the Feynman rules in light-front perturbation theory, in particular in the energy denominators. Then, the bare mass is rewritten as m 2 0 = m 2 − δm 2 in the result, and a Taylor expansion at small mass shift δm 2 is performed, assuming δm 2 ∼ α s . Mass renormalization then amounts to imposing an extra condition in order to determine δm 2 . When following such a bare perturbation theory approach for the longitudinal photon wave function, the energy denominator ED LO becomes By contrast, in the renormalized perturbation approach, one uses the renormalized mass m when writing down the diagrams, and instead includes mass counter-terms in the light-front Hamiltonian, corresponding to additional twopoint vertices. In that context, the contributions δm 2 /2k + 0 and δm 2 /2k + 1 in Eq. (46) are reinterpreted as coming from the NLO diagrams obtained by inserting such a mass counterterm on the quark or antiquark line respectively in the LO diagram in Fig. 1. The important point to note is the doubling of the energy denominator for these terms. Hence, these terms are enhanced in the limit ED LO → 0, corresponding to the on-shell limit for the quark-antiquark Fock state. In the bare perturbation theory approach, the doubling of the energy denominator comes naturally from the Taylor expansion. In the renormalized perturbation theory approach, on the other hand, it comes from the fact that one should include an energy denominator both before and after the insertion of the mass counterterm. Only NLO corrections coming with an extra copy of the energy denominator ED LO can thus be absorbed into the quark mass in the energy denominator, via mass renormalization.
In the massless quark case, the extra copy of the energy denominator ED LO disappears in the course of the evaluation of diagram (a) (see for example section IV. in Ref. [16]), showing that no quark mass can be radiatively generated by such a self-energy diagram, and that the UV divergences encountered in the massless case have nothing to do with the renormalization of the quark mass in the energy denominator.
Let us now come back to the expression (44) for the quark self-energy diagram (a). Remember that one copy of the energy denominator ED LO is already included in the LO wavefunction, and that ED LO is proportional to (P 2 + Q 2 + m 2 ) (see Eq. (29)). One finds that only the second term of (44) exhibits a doubling of the energy denominator and can thus be associated with mass renormalization. Note also from Eq. (37) that ∆ 1 depends on ED LO , and that ∆ 1 → ξ 2 m 2 in the on-shell limit (P 2 + Q 2 + m 2 ) → 0. We now add and subtract a term with where now only the first term can contribute to mass renormalization. By contrast, the UV divergence in the second term is one that should cancel at the cross section level, like in the massless case. The terms on the second line in (47) are UV finite.
Other contributions to mass renormalization in the energy denominators come from the "self-induced inertia" or "seagull" diagrams [50,62,63]. These diagrams, which correspond to instantaneous self-energy loops, are highly sensitive to the details of the UV regularization procedure. They bring a pure extra power of the energy denominator ED LO , and can thus be entirely absorbed into the quark mass in the energy denominator, via mass renormalization. We plan to present a complete analysis of mass renormalization in light-front perturbation theory at one loop in a separate future work, with a special emphasis on UV regularization issues and on "self-induced inertia" contributions.
Here, in order to avoid entering further into these issues, we choose the on-shell scheme for mass renormalization. In the context of the present calculation, this scheme is defined by imposing a strict cancellation between all the enhanced contributions in the on-shell limit ED LO → 0, or equivalently (P 2 + Q 2 + m 2 ) → 0. 6 Hence, the quark mass counterterm is chosen in order to cancel exactly the "self-induced inertia" insertions on the quark line and the first term in Eq. (47) from the diagram (a). Adding the "self-induced inertia" contributions and the quark counterterm to the diagram (a) and choosing the on-shell scheme for mass renormalization thus simply leads to a result where the first term of Eq. (47) is absent We keep calling this expression the (mass renormalized) contribution of diagram (a), by a slight abuse of language. The same treatment is done on the antiquark line, with diagram (b).
In the case of the transverse photon wavefunction, which we will study in a future publication, the quark mass also appears in the numerator in the LO expression, in the transverse photon to quark-antiquark vertex. Hence, at NLO, one has to deal with mass renormalization in the numerator as well. This involves vertex correction diagrams analog to diagrams (c) and (d), and another quark mass counterterm which is a three-point vertex. In light-front perturbation theory, mass renormalization corrections in either the denominator or the numerator thus occur in a very different pattern. In principle, the final result for the mass shift should be the same in both case, but this property is typically lost if the UV regularization procedure does not preserve the full Poincaré symmetry [56,[58][59][60]68].
The expression (48) for the mass renormalized contribution from the self-energy diagram (a) can now be written as 6 Note that the total energy of the quark-antiquark state is positive. This means that the renormalization condition is determined at a time-like virtuality for the photon Q 2 < 0, away from the physical (spacelike) region for the DIS process.
Using Eq. (45), the form factor V L (a) can be evaluated as In the above expression, the factor (D s − 4)/2(D − 4) is the regularization scheme dependent coefficient coming from the following integral and the function I V (a) is defined as This integral is ξ → 0 finite and we could in principle perform the integration over ξ analytically. However, it turns out that it is more convenient to Fourier transform first and then perform the remaining integral numerically. 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), i.e. by making the substitution z → 1 − z and P → −P in Eq. (49) simultaneously. This yields where and Summing the expressions in Eqs. (49) and (53), we obtain for the full contribution of the one-loop quark self-energy to the γ * where the NLO form factor V L (a)+(b) can be written as

D. Vertex and instantaneous contributions
We now proceed to calculate the one-loop vertex correction diagram (c) shown in Fig. 3 and the instantaneous diagram (e) shown in Fig. 4. For diagram (c), the momentum space expression of the LCWF can be written as where the light cone energy denominators ED v , ED a and ED LO are given in Eqs. (34), (33) and (27), respectively. 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, where the detailed calculation of the numerator in Eq. (59) is found in Appendix C, and gives with f (D s ) = 1 + (D s − 4)/2. At tree level, the longitudinal photon to quark-antiquark splitting Eq. (31) is proportional to the light cone helicity conserving Dirac structureū(0)γ + v(1) ∼ δ h 0 ,−h 1 . At one loop level, the helicity structure is more complicated, since in addition to a correction to the light cone helicity conserving structure, the result also contains a light cone helicity flip term ∼ δ h 0 ,h 1 . However, to calculate the NLO cross section the one-loop wavefunction in the amplitude will be convoluted with the tree level one in the conjugate amplitude. Since also the eikonal interactions with the target conserve light cone helicity, the two different helicity structures do not interfere. Thus, in fact, the light cone helicity flip term does not contribute to the cross section at this order in perturbation theory. At NNLO the square of this term would contribute to the cross section. It will therefore be convenient to separate out the two helicity structures at this point in the calculation. To do this, we split the result Eq. (60) into two parts where the NLO form factor V L (c) , which factorizes from the LO contribution, is given by with The second term in Eq. (62) is the light cone helicity flip term that appears only in the massive quark case, and it contributes to the form factor S L in Eq. (32). This term can be simplified to with It turns out that the most efficient way to compute the integrals I L (c) and is to rewrite them as a linear combination of tensor and scalar integrals 7 . This procedure, in both regularization (FDH and CDR) schemes, gives and where the one-loop vector integral B i (∆ 1 , ∆ 2 , L) = L i B 1 (∆ 1 , ∆ 2 , L) and the scalar integrals B 0 and B 1 are given in Ref. [16]. The integrals B 0 and B i in Eqs. (67) and (68) are both UV finite, and therefore all the UV divergences are carried by the scalar integral A 0 (∆ 2 ), where the coefficient ∆ 2 (defined in Eq. (37)) is independent of P. This is particularly convenient, since we have to later Fourier transform the P-dependence of the LCWF in Eq. (62). In addition, since only the integral A 0 is UV divergent, the scheme dependent part comes from the term proportional to I L (c) ∼ (D s − 4)A 0 (∆ 2 ) in Eq. (67), and in the other terms we can set D s = 4 and f (D s ) = 1. Since the vector integral B i (∆ 1 , ∆ 2 , L) is proportional to P i , the helicity flip contribution can be written as with S L (c) = 2z The (UV and IR finite) light cone helicity flip contribution (69) will not contribute to the cross section at NLO, since it cannot interfere with the LO LFWF which is helicity non-flip. Hence, we will from here on concentrate only on the helicity conserving part that does contribute to the cross section at NLO. The helicity flip term can be used, as discussed in Appendix H, to rederive the known result for the one-loop Pauli form factor of a massive quark, serving as an additional cross check of our result.
To proceed with the helicity conserving part, we first write the expression in Eq. (63) as a sum of two terms where the first and the second term in r.h.s. of Eq. (71) are the UV divergent and the UV finite pieces of V (c) , respectively. Following the notation of Eq. (67), these two terms can be written as and At this stage, we could carry on and evaluate explicitly the UV divergent term above, but it turns out to be convenient to first add a contribution coming from the instantaneous diagram (e). Using the notation shown in Fig. 4, the contribution of the instantaneous diagram (e) to the LCWF is given by where the numerator can be simplified to As discussed in [16,19] (but now in the massive quark case), the LCWF can be split up into two UV divergent contributions according to the direction of the momentum flow along the instantaneous gluon line. It is straightforward to show that these two contributions take the following form [16]: where, in the first term (e) 1 , the light cone momentum of the gluon line is flowing upwards into the quark line and in the second term (e) 2 in the opposite direction. The UV divergent coefficient corresponding to the first term simplifies to The coefficient V L (e) 2 , where the light cone momentum of the gluon line is flowing downwards into the antiquark line, is obtained by making the substitution z → 1 − z in Eq. (77).
We can now calculate the sum of diagrams (c) and (e) in momentum space. First, summing the UV divergent pieces of the vertex correction diagram in Eqs. (72) and the (e) 1 part of the instantaneous diagram in (77) together we obtain Here, the unphysical term dξ/ξ 2 cancels out between the A 0 (∆ 2 ) terms in the diagrams (c) and (e), and the remaining integrals in Eq. (78) can be performed analytically. Using the expression in Eq. (45) for A 0 (∆ 2 ), we can rewrite the expression above as For the remaining ξ integral it is convenient to first factorize the P independent coefficient ∆ 2 with respect to ξ as where the zeroes in ξ are given by The square root in Eq. (81) is associated with the threshold for massive quark pair creation if we were interested in the timelike photon case. In the spacelike case of interest note that ξ (+) > 1 and ξ (−) < 0 (which can become equalities for massless quarks) so that both zeros of ∆ 2 sit outside of the integration range in ξ. Utilizing these observations, we find the result where closed analytical expressions for the integrals I ξ;1 , I ξ;2 and I ξ;3 are given in Appendix F (see Eqs. (F1), (F2) and (F3)). The computation of the UV finite part of the coefficient in Eq. (73) is a bit more 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 over the required Fourier transform. Instead, we Feynman parametrize the denominator appearing in the B 0 and B i as where the denominator can be rewritten as Now, from Eqs. (83) and (84), we find 8 in which the integral I + is defined as Substituting the expression in Eq. (85) into Eq. (73) yields where we have defined the function I V (c) as with the coefficient The first term in r.h.s. of Eq. (87) is the same integral appearing in the massless case [16]. This integral can be done analytically and it contains both single and double logs in α, but no dependence on P (thus this contribution factors out of the Fourier transform). The second term in r.h.s. of Eq. (87) is an additional UV and ξ → 0 finite contribution coming from the massive quarks. In Eq. (88) we could in principle perform the integration over the ξ and the Feynman parameter x analytically, but it is more convenient to Fourier transform first and then perform the remaining integrals numerically.
Collecting now the contributions from Eqs. (82) and (87) together, we obtain the result where the NLO form factor V L (c)+(e) 1 can be simplified to Here, the function Li 2 (z) is the dilogarithm function, defined as The final result for the full NLO vertex and instantaneous contribution in the momentum space can be obtained by first computing the contribution from the diagrams (d) + (e) 2 , and then adding the obtained result together with the contribution in Eq. (90). The first part of the given task is most easily obtained by using the symmetry between diagrams (c) + (e) 1 and (d) + (e) 2 , i.e. by making the substitution z → 1 − z and P → −P simultaneously in Eq. (90). This yields where the NLO form factor V L (d)+(e) 2 can be written as Here, the function I V (d) is given by where the coefficients ∆ 1 , ∆ 2 and C L m reduce to and All in all, the final result for the full NLO vertex and instantaneous contribution in the momentum space is obtained by summing Eqs. (90) and (93) together. After some amount of algebra, we obtain where the NLO form factor V L (d)+(d)+(e) can be written as Here, we have again used a compact notation by introducing the functions Ω V (γ; z) and L(γ; z), which are given by and L(γ; z) = Li 2 1 1 − 1 2z (1 − γ) with In the massless limit (i.e. γ → 1), the coefficient Ω V , integrals I V (c) , I V (d) and the light cone helicity flip term vanish, and the function L satisfies where the sum of two dilogarithm function can be simplified by using the following identity Using these observations one sees that in the massless limit Eq. (98) simplifies to the result obtained in [16] and [19].
where the full form factor V L at NLO simplifies to with vanish. Using the massless limit of L(γ; z) in Eq. (103) it is then easy to see that the LCWF reduces to the known result of Ref. [16,19]. For clarity, in here we do not show explicitly the light cone helicity flip term Ψ γ * L →qq h.f. , since this contribution vanishes at the cross section level.

B. Fourier transformation to mixed space
We now Fourier transform the full NLO result of γ * L → qq LCWF in Eq. (105) into mixed space by first using the explicit expression Eq. (22) for the leading order LCWF in the factorized form of the NLO result (90). We then factor out the exponential dependence on the center-of-mass coodinate of the dipole and the momentum of the photon as in Eq. (18). This yields where the reduced NLO LCWF in mixed space simplifies to and the Fourier transformed form factor V L is given by Putting these points together, we find the following result where γ E is the Euler's constant and I V is the Fourier transformed version of the integral I V in Eq. (107) given by the following expression with and Here the coefficients C L m and C L m are defined in Eqs. (89) and (97), and κ and κ ′ are defined as: Interestingly, the expression for I V (c) (z, x 01 ) can be simplified by performing the following chain of changes of variables: x → y ≡ ξ + (1 − ξ)x, ξ → η ≡ ξ/y, and finally η → χ ≡ z(1 − η). One then arrives at The corresponding expression is obtained by the replacement z ↔ (1−z), accompanied for convenience by the change of variable χ ↔ (1−χ). This leads to This concludes our calculation of the one loop longitudinal photon to quark-antiquark LCWF for massive quarks, needed for the virtual correction to the DIS cross section. We will now proceed to the quark-antiquark-gluon final state needed for the radiative corrections.

VII. TREE LEVEL GLUON EMISSION WAVE FUNCTION
We then move to the tree level wave functions for gluon emission from a longitudinal photon state, which are needed for the full cross section at NLO. As in the massless case [17,19], we need to calculate two gluon emission diagrams shown in Fig. 5.
Applying the diagrammatic LCPT rules, we obtain for the diagram with gluon emission from the quark, diagram (f), the following expression in momentum space where the energy denominators appearing in Eq. (118) can be written as and the coefficients introduced in the denominators are defined as Similarly, for the diagram with emission from the antiquark (g), we find where the energy denominators in Eq. (122) are given by and the coefficients We can extract the transverse momentum dependence from the spinor and polarization vector structures in Eqs. (118) and (122) by using the spinor matrix element decompositions given in Eq. (A2). This procedure gives and Inserting the expressions above into Eqs. (118) and (122), we find for the sum of the diagrams (f) and (g) the result where the Dirac structures M j (f ) and M j (g) are defined as: (129) A. Fourier transform to coordinate space The Fourier transform to mixed space of the γ * L → qqg LCWF was defined by the expression (17). To simplify the Fourier transformation, we first make the following change of variables (k 1 , k 2 ) → (P, K) for the contribution coming from the diagram (f) in Eq. (128) and correspondingly the following change of variables (k 0 , k 2 ) → (P, K) for the contribution coming from the diagram (g) in Eq. (128): Next, performing the integration in Eq. (17) over the delta function of transverse momenta, we obtain the expression where the LCWF's in Eq. (128) are written in terms of the new variables P and K and the following compact notation is introduced Making these simplifications, we can write the full Fourier transformed quark-antiquark-gluon LCWF as where the reduced wavefunction reads In the above expression, we have introduced the D-dimensional Fourier integrals of the type I i and I, which are defined in Eqs. (E1) and (E2), respectively, in Appendix E. For these integrals, the following compact notation has also been introduced , , ω (g) , λ (g) ), with the coefficients We now have the full expressions for the gluon emission wavefunctions with the massive quarks. In addition, it is straightforward to check that in the massless quark limit the expression in Eq. (136) reduces to the massless quark result obtained in Refs. [17,19].

VIII. THE DIS CROSS SECTION AT NLO
We now use the wave functions to compute the DIS cross section at NLO in the dipole factorization framework.

A. Quark-antiquark contribution
Let us first write down the qq contribution to the DIS cross section at NLO. Applying the formula for the cross section in Eq. (19), we find the following expression At the accuracy we are working in here, i.e. up to terms O(α em α 2 s ), the wave function ψ γ * L →qq should be taken as ψ γ * L →qq NLO , neglecting the α 2 s contribution from the square of the loop corrections. Using the expressions for the LCWFs in Eqs. (31) and (109), it is straightforward to obtain where V L is given by Eq. (111). Adding everything together, the qq contribution to the total cross section can be written as Here, and in the following sections, we use the following the notation to denote (D − 2) and 2-dimensional transverse coordinate integrations B. Quark-antiquark-gluon contribution The qqg contribution to the DIS cross-section at NLO is given by the second term in Eq. (19). This simplifies to Using the expression in Eq. (136) for the gluon emission LCWF in mixed space, we obtain for the LCWF squared the result where we have defined the function K L qqg as The computation of individual terms in Eq. (145) follows closely the detailed derivation presented in the case of massless quarks. Therefore, for a detailed discussion, we refer the reader to our previous work [17,19]. Finally, inserting the result in Eq. (144) into Eq. (143), we obtain (146)

C. UV subtraction
Since the UV renormalization of the coupling g is not relevant at the accuracy of the present calculation, the remaining UV divergences have to cancel between the virtual qq and the real qqg contributions on the cross section level. Due to the complicated analytical structure of the gluon emission contribution in Eq. (146) 9 , 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 Eq. (145), the first and the second term are UV divergent when x 2 → x 0 and x 2 → x 1 , respectively. All the other terms are UV finite and we can immediately take the limit D = 4 at the integrand level. In order to subtract the UV divergences, we will follow the same steps as presented in [17,19]. The general idea used in these works to subtract the UV divergences rely in the following property of Wilson lines at coincident transverse coordinate points which implies that lim x 2 →x 0 S 012 = lim Thus, the UV divergences in Eq. (146) are subtracted by replacing the first and the second term with where the subtraction terms are given in terms of a single function I i UV : Now the function I i UV has to be a good UV approximation of the full integral, i.e. it must satisfy lim from which it follows that: It is important to note that there is no unique choice for the UV divergent subtraction in Eq. (152). The only requirement for the subtraction is that the UV divergence between virtual and real parts needs to cancel. Thus, it is sufficient for the subtraction to approximate the original integrals by any function that has the same value in the UV coordinate limits (for any D). Because of this cancellation, the integrals of the expressions inside the curly brackets in Eqs. (149) and (150) are finite, and one can safely take the limit D s = D = 4 under the x 2 integral.
In an arbitrary dimension D, the integral I i (see Eq. (E5)) is 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 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 [19], 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 → ∞. Another option would be to follow the polynomial subtraction scheme used in [17] (see also discussion in Appendix E of [19]). Here the subtraction function is polynomial in r. This however, introduces a new IR divergence, which must be compensated with another subtraction. For the massive quark case, we present the derivation in the polynomial subtraction scheme in Appendix G.
Proceeding with the exponential subtraction scheme we substitute Eq. (156) into Eq. (155). This gives an explicit expression for the UV approximation of the full integral In Eqs. (149) and (150) we will need the square of the UV approximation I i UV , which must be integrated over r in D − 2 dimensions. This integral can be performed using the following result

D. UV subtracted results
Following the calculations in Sec. VIII C, we then obtain for the UV subtraction terms and Next, introducing the same parametrization as in the qq contribution k + 0 = zq + , k + 1 = (1 − z)q + and k + 2,min = αq + , and changing the variables from (k + 1 , k + 0 ) → z, the sum of Eqs. (159) and (160) yields an expression for the UV divergent qqg subtraction contribution This expression precisely cancels the scheme dependent UV-divergent part in the qq contribution in Eq. (109). The remaining terms in Eq. (141) are UV finite and regularization scheme independent.
In the case of the exponential subtraction scheme, the sum of two UV finite terms in Eqs. (149) and (150) (inside the curly brackets) can be simplified to Here we have introduced the notation: and G (n:m) for the integrals that appear. These integrals could be seen as generalizations of the integral representation of Bessel functions that appear in the massless case. They could doubtlessly be transformed in many ways, but since these integrals are very rapidly converging at both small and large values of the integration variable, they should be well suited for numerical evaluation as is. The corresponding result in the case of the polynomial subtraction scheme is written down in Appendix G.

IX. LONGITUDINAL PHOTON CROSS SECTION
We can now gather here the main result of our paper, which is the the longitudinal virtual photon total cross section at NLO with massive quarks. This cross section can be written as a sum of two UV finite terms where the first term in Eq. (165) is the mass renormalized and the UV subtracted qq contribution, which is obtained by adding the UV subtraction term in Eq. (161) into Eq. (141). This gives where the functions Ω V (γ; z), L(γ; z) and I V (z, x 01 ) are explicitly writen down in Eqs. (100), (101) and (112), respectively. The first term, not proportional to α s , is explicitly the known leading order cross section for massive quarks (see e.g. [34]). The second term in Eq. (165) is the UV finite qqg contribution, which is obtained by replacing the first two terms in Eq. (146) with the subtraction term derived in Eq. (162). This can be simplified to involving the generalized Bessel function integrals from Eqs. (163) and (164). As discussed in more detail in Secs. VI A and VII A, in the limit of zero quark mass these expressions reduce to the known results in Refs. [16,17,19].

X. CONCLUSIONS
In conclusion, we have here calculated, we believe for the first time in the literature, the one-loop light cone wave function for a longitudinal photon splitting into a quark-antiquark pair including quark masses. Such a wavefunction is a central ingredient in any NLO calculation in the small-x dipole factorization formulation for processes involving heavy quarks. In particular, while we concentrated here on the total cross section, this also includes many possible diffractive or exclusive cross sections that will be important parts of the physics program at future DIS facilities.
Our result includes a renormalization of the quark mass in a pole mass or on-shell scheme. The peculiarity of the longitudinal photon polarization state is that vertex correction type diagrams do not contribute to mass renormalization, since at tree level the longitudinal photon vertex does not have a term proportional to the quark mass. This issue will be different for the transverse polarization which we intend to return to in future work. There, one will have to address the issue of the consistency in the mass renormalization between the propagator and vertex correction diagrams. We plan to revisit the issue of quark mass renormalization in LCPT much more thorougly in a separate paper.
After obtaining the one loop LCWF we also Fourier-transformed our result to mixed transverse coordinatelongitudinal momentum space. Combined with the tree level wavefunction for a quark-antiquark-gluon state this enabled us to obtain an explicit expression for the total longitudinal photon cross section in the dipole picture. The obtained cross section still has, just like the one for massless quarks, a high energy divergence when the longitudinal momentum of the gluon becomes small. This divergence (or large logarithm) needs to be further absorbed into a BK-or JIMWLK-evolution of the Wilson lines describing the target. This factorization procedure, together with NLO renormalization group evolution, will need to be developed in a consistent way to confront our calculations with experimental data. We have not elaborated on this issue here, since this will proceed similarly (and have similar problematic issues) as for the case of massless quarks, discussed in previous works [17,18,24,33].
The case of the transverse photon (virtual or real) is even more important for phenomenology, but much more complicated algebraically. It will be addressed in a forthcoming separate publication, where we will calculate the light cone wave function and the DIS cross section for transverse virtual photons. The combination of the two results will enable the calculation of a variety of DIS process cross sections in the dipole picture at NLO with massive quarks. In particular, this includes the heavy quark structure function F c 2 , which could be expected to be crucial observable for the physics of gluon saturation at the future Electron-Ion Collider EIC [69,70]. In this Appendix, we describe how to express QED/QCD emission vertices in the helicity basis by decomposing the spinor structure of given vertex to the light cone helicity non-flip and flip components.
Following the discussion presented in [19], the decomposition for the light cone vertex (without coupling and color structure) in the LC gauge can be expressed as 10 where χ and ω can be either positive or negative massive spinors, i.e. u or v. Note that we are now dealing with on-shell momenta and polarization vectors in the light cone gauge. We also use three-momentum conservation, with the appropriate signs depending on whether χ and ω are negative or positive energy spinors. The photon here is an incoming one with polarization vector ε λ (q), the corresponding expressions for an outgoing photon can be obtained by complex conjugation. To be explicit, momentum conservation means thatq =k 1 +k 2 for pair production χ(1) = v(1), ω(0) = u(0) or χ(1) = u(1), ω(0) = v(0). But for gauge boson absorption by a quark, χ(1) = u(1), ω(0) = u(0) and momentum conservation meansq +k 0 =k 1 . Vice versa, for gauge boson absorption by an antiquark, χ(1) = v(1) and ω(0) = v(0) momentum conservation meansq +k 1 =k 0 . Using the Dirac equation and Clifford algebra, it is straightforward to show that Eq. (A1) can be decomposed into three independent spinor structures where terms appearing in the first line are the light cone helicity non-flip components and terms appearing in the second line are the helicity flip components. The (∓) sign of the mass term is determined from: Here, the obtained result in Eq. (A2) is valid in arbitrary spacetime dimensions, but we emphasize that it relies on plus and transverse momentum conservation.
Appendix B: Numerator for the quark self-energy diagram (a) In this Appendix, we present some details for the calculation of the quark self-energy diagram (a). The numerator for the vertex correction diagram (c) is discussed in the following section.
The numerator for the diagram (a), where sum over the internal helicities, gluon polarization and color is implicit, can be written as Using the decompostion in Eq. (A2) with kinematical variables as in Fig. 2(a), we can express the spinor structures inside the square brakets in the helicity basis as and where the variable K is defined in Eq. (35). Using the completeness relation for the spinors and noting that γ + (k / ′ 0 + m)γ + = 2(k + 0 − k + )γ + and γ + (k / ′′ 0 + m)γ + = 2k + 0 γ + , we obtain In Eq. (B5) the terms linear in the transverse integration variable K vanish due to rotational symmetry. The remaining tensor contractions are evaluated as follows: First, the tensor product of transverse momentum variable K is replaced , which is true under the integration over the (D −2)-dimensional transverse momentum integral The gluon polarization vectors and the transverse gamma matrices appearing in Eq. (B5) are kept in (D s −2) dimension with D s > D. Thus, summing over the helicity states of the gluon and performing the remaining tensor contractions Finally, using the parametrization k + /k + 0 = ξ and definition of the leading order QED photon splitting vertex in Eq. (23) gives Appendix C: Numerator for the vertex correction diagram (c) The numerator for the diagram (c), where again sum over the internal helicities, gluon polarization and color is implicit, can be written as Using the decomposition in Eq. (A2) with the kinematical variables as in Fig. 3(c), we find where we have introduced the variable Substituting Eqs. (C2) and (C3) into (C1) and following the same steps as in B, we obtain the following result where we have defined the coefficients a (c) , b (c) and c (c) as Let us then use these results by considering first the general integral in Eq. (D7) with ∆ = Q 2 + m 2 . In this case, the result reads For the NLO contribution, we also need (in addition to the result in Eqs. (D8)) the following set of transverse Fourier integrals in D = 4: Here, the coefficient Ψ 0 (1) = −γ E (Euler's constant) and the variables (ξ, z, x) ∈ [0, 1], ∆ 1 , ∆ 2 > 0 and L 2 = (1 − ξ) 2 P 2 > 0. The coefficient κ is defined as Appendix E: Fourier transform integrals for the quark-antiquark-gluon Fock state For the gluon emission diagrams from a longitudinal photon state, we need to calculate the following two Fourier integrals and First, using the Schwinger parametrisation, Eq. (D2), for the denominators appearing in Eqs. (E1) and (E2), and then performing the remaining transverse momentum integrals by using the (D − 2)-dimensional Gaussian integral Eq. (D4) yields and By making the change of variables u = s + tω and then changing the order of integration leads to the result In an arbitrary dimension D, the integral over t would then give a dependent of incomplete Gamma function, preventing us from expressing the final result in terms of familiar special functions, e.g. the modified Bessel functions. These forms, however, serve as a sufficient starting point for deducing the appropriate UV subtractions in Sec. VIII C.
This leads to the UV approximation of Eq. (155). This approximation has a Coulomb tail at large distances in r, leading to an IR divergence which is absent when using the original expression from Eq. (155). Hence, in [17], an extra IR subtraction has been included, in order to turn the Coulomb tail into a dipole tail, thus avoiding the appearance of the unphysical IR divergence. All in all, in this scheme, the UV subtraction procedure can be written as while all of the other contributions to the cross section are the same as in the exponential subtraction scheme.
Appendix H: S L and the Pauli form factor In this appendix, a cross-check of our results is provided, by comparison with the literature. The usual parametrization of the γe − e + or γqq vertex function in QED and/or QCD (based on Lorentz and gauge invariance, and discrete symmetries such as parity) can be written as Γ µ (q) =F D (q 2 /m 2 )γ µ + F P (q 2 /m 2 ) q ν 2m iσ µν (H1) with the Dirac and Pauli form factors, and σ µν ≡ (i/2) [γ µ , γ ν ]. The relation (H1) relies on energy and momentum conservation at the vertex, and requires the two external fermion lines to be on mass shell. The photon virtuality q 2 is thus the only scale, apart from the fermion mass m and QCD non-perturbative scales.
Using momentum conservation k µ 0 + k µ 1 = q µ , one can rewrite Eq. (H2) after some algebra as u(0)Γ + (q)v(1) = F D (q 2 /m 2 ) + (q + ) 2 4k + 0 k + 1 F P (q 2 /m 2 ) u(0)γ + v(1) + F P (q 2 /m 2 ) (k + 0 −k + 1 ) 4m using the same notations P and z as in the rest of the present article (see Eq. (28)). Moreover, using momentum conservation k µ 0 + k µ 1 = q µ and the on-shell conditions k 2 0 = k 2 1 = m 2 , the photon virtuality q 2 can be expressed in terms of P and z as On the other hand, in section V, the initial-state LCWF for γ * L → qq as been calculated at NLO accuracy in QCD. It has the general form with V L and S L collecting the helicity non-flip and helicity flip contributions respectively. One recognizes the same two Dirac structures in Eqs. (H3) and (H5). Apart from the normalization, there is however a major difference between the vertex function (H3) and the LFWF (H5) : only the + and transverse components of the momentum is conserved in the splitting in the LFWF. Due to the absence of the conservation of the − component of the momentum, the relation (H4) is not valid for the LFWF, so that V L and S L a priori depend on q 2 , P 2 and z independently. Indeed, the relation (H4) is equivalent to P 2 + m 2 + Q 2 = 0, meaning (ED LO ) = 0.
For that reason, the two coefficients V L and S L contain more information than the Dirac and Pauli form factors, and can be related to them only when imposing the relation (H4). In such a way, one obtains the constraints which can be used to cross-check our results for the γ * L → qq LCWF at NLO with massive quarks. Note that V L and S L depend on q 2 , P 2 and z independently (and m 2 ), and we are then imposing one single relation between them, whereas the form factors depend only q 2 /m 2 , so that the z dependence has to drop. Due to this observation, the relations Eqs. (H6) and (H7) impose very strong constraints on V L and S L .
The form factor S L receives contributions only from the two non-instantaneous vertex correction diagrams Fig. 3 Note that the expression Eq. (H9) for S L (c) is fully finite, both in the UV and at ξ = 0, which is expected since S L is absent at tree level.
Using the Feynman parametrization, one can write B 0 and B j as (H10) Furthermore, using the change of variable x → y = ξ + (1−ξ)x, these expressions can be simplified into where in the second line we have again performed the change of variable ξ → η = ξ y . Up to this point, we have considered S L (c) in the general kinematics relevant for the LCWF. Let us now impose P 2 = −Q 2 − m 2 in order to study the correspondence with the Pauli form factor. In that case, the integral over y In order to have an expression valid also in the time-like case, one can restore the i0 by looking at the relative sign of q 2 and i0 in the energy denominators at the beginning of the calculation. Then, This is indeed the known result for the Pauli form factor at one loop in QCD, which is identical to the QED result [71,72] up to the replacement α s C F ↔ α em e 2 f .