QED corrections to parton distributions and Altarelli–Parisi splitting functions in the polarized case

We discuss the effect of QED corrections in the evolution of polarized parton distributions. We solve the corresponding evolution equations exactly to O (α) and O ( a 2S ) in Mellin N -space, extending the available techniques for pure QCD evolution. To accomplish this, we introduce, for the ﬁrst time, the Altarelli–Parisi polarized kernels at LO in QED. Furthermore, weperform aphenomenological analysis of the QED effects on polarized parton distributions (pPDFs), proposing different scenarios for the polarized photon density. Finally, we quantify the impact of the corresponding QED contributions to the polarized structure function g 1 . We show that the relative corrections to both the pPDFs and the g 1 structure function are approximately at the few percent level, which is the order of magnitude expected considering the value of α .


I. INTRODUCTION
During the last few decades, there has been a great deal of activity in the area of high-energy physics with polarized processes, both from the experimental and theoretical points of view.The interest in this area arises from the need to understand the way in which the proton spin emerges from the share of its constituents, an open question and a key area of research in particle physics.
The spin content of the proton can be encoded in terms of the polarized distributions of quarks and gluons (pPDFs) [1], which can be probed experimentally in high-energy collision processes with polarized nucleons.Despite the significant progress, due to the lack of measurements of observables probing a large kinematical range, associated with the challenges involved in the production of polarized beams, the situation is still unclear, particularly concerning the polarized gluon distribution.
Since early measurements at the polarized high-energy colliders began, the fixed-target polarized lepton-nucleon deep inelastic scattering (DIS) experiments performed have shown that a relatively small amount of the proton spin is carried by the quark and antiquark spins [2][3][4][5][6][7] .
However, these types of measurements are poorly sensitive to gluons.Instead, the best ∆g probes have so far been provided by the polarized proton-proton collisions available at the Relativistic Heavy Ion Collider at Brookhaven National Laboratory [8], where various processes such as jet or hadron production at high p T transverse momentum receive substantial contributions from gluon-induced hard scattering.These measurements improved the description of the gluon spin distribution and allowed a better understanding of the proton spin structure.In particular, global analyses performed by the DSSV [9,10] and NNPDF groups [11], conclude that the contribution of ∆g to the proton spin is not negligible, although providing constraints only for a reduced range of proton momentum fractions.
In this context, the future Electron-Ion Collider (EIC), which will enable a much wider kinematic range and achieve unprecedented precision for polarized measurements [12,13], is expected to provide new insights into the spin share of the proton in terms of its fundamental components.
Considering the forthcoming measurements, it is essential to increase the precision of the theoretical calculations.Since QCD computations are reaching high levels of accuracy, even N3LO in some cases, other effects that were not previously taken into account, e.g.QED effects, start to play a more significant phenomenological role in the theoretical predictions.Given that α ≈ α 2 s , NLO QED corrections compete with those at NNLO in QCD.These corrections typically manifest at the few percent levels for many observables, ultimately becoming quantitatively significant for achieving an accurate description.
In addition to the need to advance in the perturbative computations, it is also important to reach the same level of accuracy on the non-perturbative side, i.e., on the unpolarized and polarized parton distribution functions (PDFs and pPDFs).As mentioned above, the study of PDFs is essential to understand the proton structure.Although PDFs cannot be derived from first principles and are instead fitted using experimental data, their evolution with respect to energy can be computed.This calculation requires the knowledge of the Altarelli-Parisi splitting functions (AP kernels) [14].The computation of the NNLO corrections to the unpolarized splitting functions performed in [15][16][17][18] and the development of modern parton distribution analyses [19][20][21][22][23], allows achieving the required accuracy in QCD.Also, the unpolarized kernels at LO in QED were calculated in [24][25][26], and at NLO in QED and mixed order QED⊗QCD in [27,28].Solutions to the evolution equations including corrections at LO in QED were presented in [24,29,30] and different fits of the parton distributions including these corrections were performed in [31][32][33].
A key density that appears for the first time when QED corrections are taken into account is the photon distribution, i.e, the probability of finding a photon in the nucleon.Previous analyses of this density were performed in [32][33][34][35], either obtaining large uncertainties or relying on phenomenological models for the contribution to the photon PDF from the low-energy regions [36].In this context, a new approach towards obtaining an accurate description of this distribution was presented in [36,37] , where it is shown that the photon density can be computed from the present knowledge of the unpolarized structure functions.
As stated before, the situation is quite different in the polarized case.While the QCD splitting functions are known at LO [14], NLO [38][39][40] and NNLO [41][42][43], the global analysis is still restricted to NLO QCD precision.To improve the accuracy in the polarized case, it is necessary to reach NNLO accuracy in QCD and include the QED corrections, as well.
In this article, we present for the first time, the expressions for polarized splitting functions and the polarized structure function g 1 to order α in QED and study the phenomenological impact of the corresponding corrections after introducing different scenarios for the polarized photon density.
Working on the Mellin space, with a procedure similar to the one developed in [44,45] we solved the evolution equations analytically.This paper is organized as follows.In section II, we present the evolution equations and establish the main notation for this article.In section III, we introduce, for the first time, the QED polarized kernels, and in section IV, we perform the analytical solution of the evolution equations in Mellin space.The phenomenological impact, including the introduction of different scenarios for the polarized photon distributions, is studied in section V. Finally, in section VI, we present the conclusions of this work.

II. SPLITTING KERNELS AND PARTON DISTRIBUTION BASIS
We start by writing down the general expression for the evolution of gluon, photon and quark distributions.In order to simplify the notation, we write the distributions and kernels in the unpolarized form.The polarized case is recovered by simply considering the corresponding polarized proton distribution and splitting functions, i.e. basically adding ∆ in all the expressions below.
with t = ln (Q 2 ) (Q being the factorization scale), n f is the number of active fermions and P ij the Altarelli-Parisi splitting functions in the space-like region.The analog evolution equations for antiquarks can be obtained by applying charge conjugation.Here, we use the conventional notation to indicate convolutions.We do not include the lepton distributions, since up to the order we reach here they basically factorize from the rest of the distributions * .Along this work we will use the expressions for the splitting functions including QCD and QED perturbative corrections.In this sense, each kernel can be expressed as, where the indices (a, b) indicate the (QCD,QED) perturbative order of the calculation, with a S ≡ α S 4π and a ≡ α 4π .Due to the QED corrections, the kernels can depend on the electric charge of the initiating quarks (up or down type).In general, we have P (n,1) i,j ∼ e 2 q with at least some i, j = q.
The quark splitting functions are decomposed as P q i qk = δ ik P V q q + P S q q , (5) which act as a definition for P V qq and P V q q .In order to minimize the mixing between the different parton distributions in the evolution, it is convenient to introduce the following basis, differentiating the following singlet and the non-singlet pPDF combinations [24]: where ∆ U D could also include the top quark distribution in case of a 6 flavour analysis (adding ∆ ct and t v to complete the basis).
Taking into account that beyond NLO in QCD the singlet non-diagonal terms (P S q q and P S qq ) differ [46], it is useful to define ∆P S ≡ P S qq − P S q q, P S ≡ P S qq where we explicitly use that these contributions do not depend on the quark charge up to the order we reach, since they do not receive QED corrections to O(α).In the expansion in powers of the couplings, the flavor-diagonal ('valence') quantity P V qq in Eq.( 6) begins at first order.On the other hand, P V q q and the flavor-independent ('sea') contributions P S qq and P S q q, and therefore the 'pure-singlet' term P S in Eq.( 12), are of order α 2 s .The kernel ∆P S in Eq.( 12) arises for the first time at the third order.In this work, we reach NLO accuracy in QCD and LO in QED, which implies that ∆P S = 0 for all our calculations.
The evolution equations Eq.( 1) for the parton distributions in the basis of Eq.( 7) read for the non-singlet sector, where in the second and third lines we use the notation P + u(d) , differentiating the kernels in two groups, those corresponding to the up quark flavours (u, c, t) and those corresponding to the down quark flavours (d, s, b).As noted above, the term ∆P S is zero up to the order we work, therefore, in the non-singlet sector we have uncoupled equations for the distributions q v i (Eq.( 13)).For the singlet sector we have, Notice that in the limit of an equal number of u and d quarks (n u = n d ) and the same electric charges (P ug = P dg , P uγ = P dγ , P + u = P + d ), ∆ U D decouples from the other distributions in the evolution, while the singlet evolution recovers the usual pure-QCD expression.

III. QCD-QED POLARIZED SPLITTING KERNELS
Since the calculation at LO in QED is analogous to that in QCD, except for the colour factors and electric charges, the LO Altarelli-Parisi kernels in QED can be derived from those in QCD by adjusting the colour factors as follows [27,28] And, if the Feynman diagram involves a loop of fermions, as in the case of ∆P (1,0) gg , one needs to apply the replacement: where the sum runs over all the fermions involved in the process.To set the correct normalization, we start by reminding the lowest order polarized splitting functions in QCD ∆P with and the usual plus distribution is defined as for any regular test function f .Therefore, setting the factors as Eq.( 20), the lowest order splitting functions in QED P (0,1) ij result in where there is an explicit dependence on the quark electromagnetic charge.Furthermore, the sum over fermion charges in the ∆P (0,1) γγ kernel corresponds to the definition with n f and n L the number of quark and lepton flavours, respectively.

IV. SOLVING THE EVOLUTION EQUATIONS
To solve the evolution equations, we work in Mellin N -space.We define the Mellin transformation, from Bjorken x-space to complex N -moment space, as and its inverse reads here C N denotes a suitable contour in the complex N plane that has an imaginary part ranging from −∞ to ∞ and that intersects the real axis to the right of the rightmost pole of a(N ).In practice, it is beneficial to choose the contour to be bent at an angle < π/2 towards the negative real-N axis [44,47].The integration in Eq.( 27) can then be very efficiently performed numerically by choosing the values of N as the supports for a Gaussian integration [44].
The advantage of working in this space is that the convolutions appearing in the evolution equations can be written simply as products, which makes it easier to manipulate the expressions and solve the corresponding equations analytically at a given order.From now on, all kernels and pPDFs will be expressed in this space, but we will keep the nomenclature a bit loose to avoid complicating the notation.

A. Non-Singlet case
We start with the solution for the non-singlet distributions (Eqs.(13),( 14),( 15)).If we express the kernels with the form of Eq.( 3), we have a generic form for the evolution where f is any distribution of the non-singlet sector, and P is the kernel corresponding to the evolution of that distribution.The solution is customarily obtained by posing an evolution operator that can be decomposed in the product of an LO QED operator and a corresponding LO and NLO one in QCD as, We will make use of the well-known RGE for the strong and electromagnetic coupling constants, to express the Eq.( 29) in terms of a and a S .Inserting the expression in Eq.( 30) into Eq.(29) we can find separate solutions for QED and QCD.Keeping the accuracy only up to NLO in QCD and LO in QED, we get Let's start with the solution for the QCD sector, if we remain at LO accuracy we have which, with the constraint E LO (Q 0 , Q 0 ) = 1, we obtain the well known exponential solution where we define a S ≡ a S (Q) and a S,0 ≡ a S (Q 0 ).For the NLO component one gets with the solution valid up to NLO accuracy as By combining the LO and NLO operators, we arrive at The same procedure can be applied to the QED sector as well.Keeping up to LO, the resulting QED evolution operator is given by where we define a ≡ a(Q) and a 0 ≡ a(Q 0 ).As long as all flavours are active, all combinations evolve as non-singlets.But if one quark flavour is not active, the corresponding ∆ ij combination containing that quark term can be expressed as a combination of singlet and non-singlet terms.

B. Singlet case
For the singlet sector, we can express the evolution equations (Eqs.(16), ( 17), ( 18) and ( 19)) in a matrix form, where f S is the vector of singlet combinations of flavors Eq.( 7) and P is the matrix of kernels, the kernel matrix can be expressed in a perturbative form, For QCD matrices, with no photon content as noticed in the last column/row, we have, where n ud ≡ n u − n d and we use that to this order, we can express: where in the first two lines i = 1, 2. The expression of the kernels P (1,0) qq , P (1,0) gq , and P (1,0) gg are obtained by performing the corresponding transform to Mellin N -space on Eq.( 22).In the third line, we use that P V (1,0) q q = 0, and to construct the LO matrix, we employ that P S(1,0) = 0.
The NLO polarized kernels P (2,0) qg , P (2,0) gq , P (2,0) gg , P S(2,0) and P + (2,0) can be found in [38][39][40].For the QED sector, with no gluon content in the third column/row, we have, where we use that to this perturbative order, we can express with the definitions of P (0,1) f γ ,P (0,1) γf deduced from Eq.( 24) by extracting the electric charge e 2 q of the corresponding kernels and taking the corresponding transform to Mellin N -space.In the first line, we use that P ±(0,1) u(d) = ∆P (0,1) qq from Eq.( 24) with e 2 q = e 2 u (e 2 d ) , which arises from the fact that P V (0,1) q q = 0 in Eq.( 6).In addition, we define One non-trivial problem arises from the fact that the evolution matrix in Eq.( 41) has noncommutative properties at different orders.Even the LO QCD matrix does not commute with the QED one.This makes more difficult to isolate the QED and QCD solutions than in the nonsinglet case.However, if mixed-order terms (O(a a S )) are thrown away, which is correct at the order we are working, we can separate the QED and QCD equations.A more detailed discussion is presented in section IV C. In summary, we obtain matrix equations analogous to Eq.( 32) that, in the case of QCD, correspond to where the matrices P (1,0) and P (2,0) do not commute, so the equation cannot be solved in a closed exponential form.The usual way to solve the equation is to combine the exact solution at LO with a power expansion in the coupling for the rest [44,45].The LO solution can be expressed as a closed exponential form, Note that in the exponent we have the matrix P (1,0) , to be able to operate in the following steps it is convenient to rewrite this expression by applying the eigenvalue decomposition of the LO splitting function matrix, where λ i are given by and the matrices e i are the projectors to the subspace corresponding to each eigenvector, , with a bit of algebra, we arrive at the following expression for the evolution operator at LO, As mentioned, we can propose a general solution of Eq.( 50) employing a series expansion around the LO result, reading [44] Keeping terms up to NLO we have the truncated calculation If we replace this solution in Eq.( 50) we can obtain the following relation for the matrix U 1 , By decomposing the matrix U 1 into the eigenvalue subspaces of P (1,0) and using the relation Eq.( 58) we can find an expression for U 1 in terms of the matrix R and the projectors e i .Then from Eq.( 57) we obtain the evolution operator, (aS,0−aS) (aS,0−aS) (aS,0−aS) Since in the QED sector we restrict our analysis to LO, the evolution equation becomes simpler and can be expressed as follows Again, the solution has an exponential form and we can express it in terms of the eigenvalues of the matrix P (0,1) (λ ′ i ) and the projectors to the corresponding eigenvector subspace e ′ 1 , e ′ 2 , e ′ 3 and −(P (0,1) ) 2 + P (0,1) (λ , , where the eigenvalues are λ 4 = 0 and λ 1 , λ 2 ,λ 3 can be found in the Jupyter notebook EVOLU-TIONQED.ipynb accompanying this manuscript.In this notebook, the eigenvalues are calculated as the three nonzero roots of the characteristic polynomial of the kernel matrix P (0,1) .

C. Combining QED and QCD evolution
Equipped with an operator for the evolution of QED and another for the evolution of QCD that satisfies two separate differential equations, we can attempt a combination of both effects.
For simplicity, we will keep the following discussion only up to the LO order in QCD † , but the result can be expanded to NLO without difficulties.We start by proposing the following solution for the Eq.( 39), where the operators E QCD and E QED are solutions of the equations Eq.( 50)(only the LO) and Eq.( 59) respectively.Taking the derivative with respect to t, and using the evolution equations of the coupling constants Eq.( 31), we get Therefore, an extra +O(a 2 S ) is implicit in all the expressions below.
with n + m = 3, (n, m ̸ = 0).We note that by neglecting the last term on the right-hand side we can get the solution we were looking for, since QED and QCD effects decouple.This term can be neglected since it is of order O(a a S ).As we mentioned in the section IV, the kernel matrices of QED (P (0,1) ) and QCD (P (1,0) ) do not commute, therefore, by posing in Eq.( 62) the operators E QCD E QED in different order we arrive at a different possible solution.However, by inverting the operators, the only change is a sign in the commutator term of order O(a a S ) that we drop.In fact, by developing the operators in power series of the couplings (using that, ln a S a S,0 ), it can be shown that the solutions have the form [29], where we keep only the lowest mixed order term, and we define Then, if we propose a symmetric solution of the form [29] Ẽ we can eliminate the mixed order term, leaving us at a higher accuracy.If we calculate the derivative with respect to t we get, with n + m = 3, (n, m ̸ = 0).Developing the logarithms in the power series of the couplings, we can see that the mixed order term cancels out, and we obtain the required solution, V. PHENOMENOLOGICAL RESULTS In this section, we will apply the evolution operators to parton distributions and study their phenomenological implications.Although we are not attempting to perform a global analysis, these results will be useful as a first visualization of the phenomenological effects of the QED corrections on the polarized parton distributions and the structure function.To carry out the analysis, we use as a probe for quarks and gluon pPDFs the DSSV18 set of helicity parton densities [10].On the other hand, due to the lack of experimental polarized data, we do not count on a distribution for the polarized photon.Therefore, in order to analyze the impact of the QED corrections, we present different scenarios for the photon pPDF, based on the unpolarized distribution NNPDF3.1luxQED[48].In that work, they performed a global analysis complemented by the LUXqed constraint that relates the photon density distribution to the lepton-proton scattering structure functions [36,37].
However, the NNPDF3.1luxQEDfit is only valid for energies higher than Q 2 = 2.72 GeV 2 .To obtain the distribution at the initial scale we want (Q 2 0 = 1 GeV 2 ), we simply evolve it to that scale.Having the photon PDF, we propose three possible schemes to model the polarized pPDF at an initial scale Q 2 0 : • A An extreme scenario, where we take the polarized pPDF identical to the unpolarized one ∆γ(x) = γ(x) • B A moderate scenario, where we assume a linear polarization, ∆γ(x) = x γ(x) • C The null scenario, where we take the photon pPDF to vanish at the initial scale.
Considering that the positivity constraint |∆f (x)| < f (x) imposes a limit on the size of the polarized density distributions [49], scheme A proposes to take the photon pPDF (at the initial scale) as the upper limit of this constraint.The scheme B assumes a simple linear relation between the polarized and unpolarized distributions, that is roughly reproduced by the quark densities.At last, in the scheme C, we study the case of setting ∆γ(x) to zero.As stated before, we combine these three schemes with the NLO set of polarized distributions from [10] and observe how the distributions are affected in the evolution due to the QED corrections and the polarized photon density.
In order to analyze the results, we define the relative QED corrections as, We evolve the DSSV pPDFs supplemented with the photon pPDF from Q 2 0 = 1 GeV 2 to Q 2 = 800GeV 2 .In the evolution, we use a variable flavor number scheme with a minimum of three flavors and set the charm and bottom threshold at Q 2 mc = 2 GeV 2 and Q 2 m b = 20.25 GeV 2 respectively (we do not take into account the top quark).
In figure 1 (top) we show the photon pPDF evolved to Q 2 = 800 GeV 2 for the different choices we made for the pPDF at Q 2 0 = 1 GeV 2 .As expected, the extreme scheme A differs significantly from the other scenarios, which provide similar results, since for them the photon distribution at high scales is mostly dominated by the emission from quarks more than from the assumption at the initial scale.In the bottom part of the figure 1 we show the ratio between the polarized and unpolarized photon distributions.
In the figure 2 we show the relative corrections for the parton distributions of u-quark, d-quark, Σ, ∆ U D and gluon as a function of x.The QED corrections are large in the regions of x where the partons distributions are small.The effect is typically of the order of the few percent, the same order of magnitude expected from the NNLO QCD corrections [50].In the cases of Σ and ū pPDFs, there is a change of the sign in the distribution, forcing the relative effect to be very large due to the smallness of the corresponding QCD distribution, Eq.(69).We also show the results for the different scenarios chosen to model the photon pPDFs.In particular, the scheme A (red) produces significantly different results from the other schemes (scheme B: grey, scheme C: blue) in the small x regions.
In order to estimate the effect of QED corrections directly on physical observables, we first calculate the structure function g 1 to first order in QED.The function g 1 can be expressed in terms of the Wilson coefficients ∆C (i,j) i at different perturbative orders, in Bjorken x-space is expressed as, Let us start by recalling the well-known coefficients to NLO in QCD [40, 51, and references therein]: To calculate the QED Wilson coefficients ∆C (0,1) i , it is necessary to adjust the colour factors of the corresponding QCD coefficients, as described in Eq.( 20).Thus, we obtain ∆C (0,1) q (x) = e 2 q (1 + In figure 3 we present the structure function xg 1 (with QED corrections) as a function of x (top), and the relative QED corrections using Eq.(69) (bottom).Again, the colours indicate the different schemes we chose to model the photon pPDF at the initial scale.Relative corrections are expected to be larger for values of x where the function g 1 approaches zero.Nevertheless, we also observe corrections at the percent level in the region of x where xg 1 reaches its peak value.On the other hand, in this particular region, there do not seem to be significant differences between the different schemes A, B and C.This suggests that the main correction in g 1 arises from the modification of the quark densities due to QED effects and, primarily, due to the incorporation of the coefficient ∆C (0,1) q in Eq.(70).In order to compare the size of the contributions of the Wilson QED coefficients, we defined the ratios, R q = ∆C (0,1) q ⊗ (∆q + ∆q) ∆C (0,1) q ⊗ (∆q + ∆q) + ∆C (0,1) γ ⊗ ∆γ and R γ = 1 − R q . (71) In the bottom right-hand side of the figure 3 we show the R i ratios as a function of x.Note that the contribution of ∆C (0,1) q is significantly higher than that of ∆C (0,1) γ . The study presented in [52] estimates the potential impact of the EIC data on both the helicity proton distributions and the , represented by the ratios given by Eq.(71).

VI. CONCLUSIONS
In this paper, we present the first-order polarized kernels in QED.Working in Mellin space, we solved the partonic evolution equations, including the photon pPDF, at O(α) and O(a 2 S ).After introducing three possible scenarios for the unknown polarized photon distribution, we first analyze the phenomenological impact of the polarized parton distributions, using the set of pPDFs DSSV18 as a test probe.Finally, we calculate the QED corrections to the structure function g 1 .
We found that the main correction in g 1 arises from the modification of the quark densities due

FIG. 1 .
FIG. 1.The photon pPDF (top) at Q 2 = 800 GeV 2 .The ratio between polarized and unpolarized PDF (down).The colours indicate the different schemes chosen for the photon pPDF as explained in the main text.