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

Daniel de Florian<EMAIL_ADDRESS>Lucas Palma Conte
<EMAIL_ADDRESS>
(a) International Center for Advanced Studies (ICAS), ICIFI & ECyT-UNSAM, 25
de Mayo y Francia, (1650) San Martín, Pcia. Buenos Aires, Argentina

Abstract

> We discuss the effect of QED corrections in the evolution of polarized
> parton distributions. We solve the corresponding evolution equations exactly
> to ${\cal O}(\alpha)$ and ${\cal O}(a_{\mathrm{S}}^{2})$ in Mellin
> $N$-space, extending the available techniques for pure QCD evolution. To
> accomplish this, we introduce, for the first time, the Altarelli-Parisi
> polarized kernels at LO in QED. Furthermore, we perform a phenomenological
> 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 $\alpha$.

##  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
originates form the necessity to understand the way in which the proton’s 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) Leader
(2013), 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 Nocera  _et al._ (2012);
Blumlein and Bottcher (2010); Leader  _et al._ (2010); Hirai and Kumano
(2009); de Florian  _et al._ (2008, 2009) . However, these types of
measurements are poorly sensitive to gluons. Instead, the best $\Delta g$
probes have so far been provided by the polarized proton-proton collisions
available at the Relativistic Heavy Ion Collider at Brookhaven National
Laboratory Aschenauer  _et al._ (2015), 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 analysis
performed by the DSSV de Florian  _et al._ (2014); De Florian  _et al._ (2019)
and NNPDF groups Nocera  _et al._ (2014), concludes that the contribution of
$\Delta 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 Accardi  _et al._ (2016), 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. 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
$\alpha\approx\alpha_{s}^{2}$, NLO QED corrections compete with those at NNLO
in QCD. These corrections typically manifest at the few percent level 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)Altarelli
and Parisi (1977). The computation of the NNLO corrections to the unpolarized
splitting functions performed in Moch  _et al._ (2002, 2004); Vogt  _et al._
(2004, 2006) and the development of modern parton distribution analysis Ball
_et al._ (2022a); Hou  _et al._ (2021); Bailey  _et al._ (2021); Abt  _et al._
(2022); Ball  _et al._ (2022b), allows to achieve the required accuracy in
QCD. Also, the unpolarized kernels at LO in QED were calculated in Roth and
Weinzierl (2004); Kripfganz and Perlt (1988); Spiesberger (1995), and at NLO
in QED and mixed order QED$\otimes$QCD in de Florian  _et al._ (2016a, b).
Solutions to the evolution equations including corrections at LO in QED were
presented in Roth and Weinzierl (2004); Bertone  _et al._ (2014); Sadykov
(2014) and different fits of the parton distributions including these
corrections were performed in Carrazza (2015); Martin  _et al._ (2005); Ball
_et al._ (2013).

A key density appearing 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 analysis of this density were performed in Harland-
Lang  _et al._ (2016); Schmidt  _et al._ (2016); Martin  _et al._ (2005); Ball
_et al._ (2013), either obtaining large uncertainties or relying on
phenomenological models for the contribution to the photon PDF from the low-
energy regions Manohar  _et al._ (2016). In this context, a new approach
towards obtaining an accurate description of this distribution was presented
in Manohar  _et al._ (2016, 2017) , 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 Altarelli and Parisi (1977),
NLOMertig and van Neerven (1996); Vogelsang (1996a, b) and NNLO Moch  _et al._
(2014); Vogt  _et al._ (2014); Blümlein  _et al._ (2021), 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
$\alpha$ 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 Vogt (2005); Mottaghizadeh  _et al._ (2017) we solve 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. For the sake of simplicity in 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 $\Delta$ in all
the expressions below.

$\displaystyle\frac{dg}{dt}$ $\displaystyle=$
$\displaystyle\sum_{j=1}^{n_{f}}P_{gq_{j}}\otimes
q_{j}+\sum_{j=1}^{n_{f}}P_{g\bar{q}_{j}}\otimes\bar{q}_{j}+P_{gg}\otimes
g+P_{g\gamma}\otimes\gamma\,,$ $\displaystyle\frac{d\gamma}{dt}$
$\displaystyle=$ $\displaystyle\sum_{j=1}^{n_{f}}P_{\gamma q_{j}}\otimes
q_{j}+\sum_{j=1}^{n_{f}}P_{\gamma\bar{q}_{j}}\otimes\bar{q}_{j}+P_{\gamma
g}\otimes g+P_{\gamma\gamma}\otimes\gamma\,,$ (1)
$\displaystyle\frac{dq_{i}}{dt}$ $\displaystyle=$
$\displaystyle\sum_{j=1}^{n_{f}}P_{q_{i}q_{j}}\otimes
q_{j}+\sum_{j=1}^{n_{f}}P_{q_{i}\bar{q}_{j}}\otimes\bar{q}_{j}+P_{q_{i}g}\otimes
g+P_{q_{i}\gamma}\otimes\gamma\,,$

with $t=\mathrm{ln}\left(Q^{2}\right)$ ($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

$(f\otimes g)(x)=\int_{x}^{1}\,\frac{dy}{y}\,f\left(\frac{x}{y}\right)g(y)\,,$
(2)

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***To $\cal{O}(\alpha)$ lepton distributions only couple, in a
trivial way, to the photon density.. 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,

$\displaystyle P_{ij}$ $\displaystyle=$
$\displaystyle\sum_{\\{a,b\\}}a_{\mathrm{S}}^{a}\,a^{b}\,P_{ij}^{(a,b)},$ (3)

where the indices $(a,b)$ indicate the (QCD,QED) perturbative order of the
calculation, with $a_{\mathrm{S}}\equiv\frac{\alpha_{\mathrm{S}}}{4\pi}$ and
$a\equiv\frac{\alpha}{4\pi}$. 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_{i,j}^{(n,1)}\sim e_{q}^{2}$ with at least some $i,j=q$.

The quark splitting functions are decomposed as

$\displaystyle P_{q_{i}\,q_{k}}$ $\displaystyle=$
$\displaystyle\delta_{ik}\,P^{V}_{qq}+P^{S}_{qq}\,,$ (4) $\displaystyle
P_{q_{i}\,\bar{q}_{k}}$ $\displaystyle=$
$\displaystyle\delta_{ik}\,P^{V}_{q\bar{q}}+P^{S}_{q\bar{q}}\,,$ (5)
$\displaystyle P_{q}^{\pm}$ $\displaystyle=$ $\displaystyle P^{V}_{qq}\pm
P^{V}_{q\bar{q}}\,,$ (6)

which acts as a definition for $P^{V}_{qq}$ and $P^{V}_{q\bar{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 Roth and Weinzierl
(2004):

$\displaystyle f^{NS}$ $\displaystyle=$
$\displaystyle\\{u_{v},d_{v},s_{v},c_{v},b_{v},\Delta_{uc},\Delta_{ds},\Delta_{sb}\\},$
$\displaystyle f^{S}$ $\displaystyle=$
$\displaystyle\\{\Delta_{UD},\Sigma,g,\gamma\ \\},$ (7)

where

$\displaystyle q_{v_{i}}$ $\displaystyle=$ $\displaystyle
q_{i}-\bar{q_{i}}\,,$ (8) $\displaystyle\Delta_{uc}$ $\displaystyle=$
$\displaystyle u+\bar{u}-c-\bar{c}\,,$ $\displaystyle\Delta_{ds}$
$\displaystyle=$ $\displaystyle d+\bar{d}-s-\bar{s}\,,$
$\displaystyle\Delta_{sb}$ $\displaystyle=$ $\displaystyle
s+\bar{s}-b-\bar{b}\,,$ (9) $\displaystyle\Delta_{UD}$ $\displaystyle=$
$\displaystyle u+\bar{u}+c+\bar{c}-d-\bar{d}-s-\bar{s}-b-\bar{b}\,,$ (10)
$\displaystyle\Sigma$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{n_{f}}(q_{i}+\bar{q}_{i})\,.$ (11)

$\Delta_{UD}$ could also include the top quark distribution in case of a 6
flavour analysis (adding $\Delta_{ct}$ and $t_{v}$ to complete the basis).

Taking into account that beyond NLO in QCD the singlet non-diagonal terms
($P^{S}_{q\bar{q}}$ and $P^{S}_{qq}$) differ Catani  _et al._ (2004), it is
useful to define

$\displaystyle\Delta P^{S}$ $\displaystyle\equiv$ $\displaystyle
P^{S}_{qq}-P^{S}_{q\bar{q}},$ $\displaystyle P^{S}$ $\displaystyle\equiv$
$\displaystyle P^{S}_{qq}+P^{S}_{q\bar{q}},$ (12)

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
$\cal{O}(\alpha)$. In the expansion in powers of the couplings, the flavor-
diagonal (‘valence’) quantity $P_{qq}^{V}$ in Eq.(6) begins at first order. On
the other hand, $P_{q\bar{q}}^{V}$ and the flavor-independent (‘sea’)
contributions $P^{S}_{qq}$ and $P^{S}_{q\bar{q}}$, and therefore the ‘pure-
singlet’ term $P^{S}$ in Eq.(12), are of order $\alpha_{s}^{2}$. The kernel
$\Delta 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
$\Delta P^{S}=0$ for all our calculations.

The evolution equations Eq.(1) for the parton distributions in the basis of
Eq.(II) read for the non-singlet sector,

$\displaystyle\frac{dq_{v_{i}}}{dt}$ $\displaystyle=$ $\displaystyle
P_{q_{i}}^{-}\otimes q_{v_{i}}+\sum_{j=1}^{n_{f}}\Delta P^{S}\otimes
q_{v_{j}}\,,$ (13) $\displaystyle\frac{d\\{\Delta_{uc},\Delta_{ct}\\}}{dt}$
$\displaystyle=$ $\displaystyle
P_{u}^{+}\otimes\\{\Delta_{uc},\Delta_{ct}\\}\,,$ (14)
$\displaystyle\frac{d\\{\Delta_{ds},\Delta_{sb}\\}}{dt}$ $\displaystyle=$
$\displaystyle P_{d}^{+}\otimes\\{\Delta_{ds},\Delta_{sb}\\}\,,$ (15)

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 $\Delta 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,

$\displaystyle\frac{d\Delta_{UD}}{dt}$ $\displaystyle=$
$\displaystyle\frac{P_{u}^{+}+P_{d}^{+}}{2}\otimes\Delta_{UD}+\frac{P_{u}^{+}-P_{d}^{+}}{2}\otimes\Sigma+(n_{u}-n_{d})P^{S}\,\otimes\Sigma$
(16) $\displaystyle+$ $\displaystyle 2(n_{u}P_{ug}-n_{d}P_{dg})\otimes
g+2(n_{u}P_{u\gamma}-n_{d}P_{d\gamma})\otimes\gamma\,,$
$\displaystyle\frac{d\Sigma}{dt}$ $\displaystyle=$
$\displaystyle\frac{P_{u}^{+}+P_{d}^{+}}{2}\otimes\Sigma+\frac{P_{u}^{+}-P_{d}^{+}}{2}\otimes\Delta_{UD}+n_{f}\,P^{S}\otimes\Sigma$
(17) $\displaystyle+$ $\displaystyle 2(n_{u}P_{ug}+n_{d}P_{dg})\otimes
g+2(n_{u}P_{u\gamma}+n_{d}P_{d\gamma})\otimes\gamma\,.$
$\displaystyle\frac{dg}{dt}$ $\displaystyle=$ $\displaystyle
P_{gq}\otimes\Sigma+P_{gg}\otimes g\,.$ (18) $\displaystyle\frac{d\gamma}{dt}$
$\displaystyle=$ $\displaystyle\frac{P_{\gamma u}-P_{\gamma
d}}{2}\otimes\Delta_{UD}+\frac{P_{\gamma u}+P_{\gamma
d}}{2}\otimes\Sigma+P_{\gamma\gamma}\otimes\gamma\,.$ (19)

Notice that in the limit of equal number of $u$ and $d$ quarks ($n_{u}=n_{d}$)
and same electric charges ($P_{ug}=P_{dg},P_{u\gamma}=P_{d\gamma}$,
$P_{u}^{+}=P_{d}^{+}$), $\Delta_{UD}$ 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
color factors and electric charges, the LO Altarelli-Parisi kernels in QED can
be derived from those in QCD by adjusting the color factors as follows de
Florian  _et al._ (2016a, b)

$C_{A}\rightarrow 0,\;\;\;\;\;\;C_{F}\rightarrow
e_{q}^{2},\;\;\;\;\;\;T_{R}\rightarrow N_{C}\,e_{q}^{2}.\;\;\;\;\;\;$ (20)

And, if the Feynman diagram involves a loop of fermions, as in the case of
$\Delta P_{gg}^{(1,0)}$, one needs to apply the replacement:

$n_{f}\,T_{R}\rightarrow\sum_{f}e_{f}^{2},$ (21)

where the sums run 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 $\Delta P_{ij}^{(1,0)}$ Altarelli and Parisi (1977)

$\displaystyle\Delta P_{qq}^{(1,0)}(x)$ $\displaystyle=$ $\displaystyle
2\,C_{F}\left[\frac{1+x^{2}}{(1-x)_{+}}+\frac{3}{2}\delta(1-x)\right]\equiv
2\,C_{F}\,\left[\delta p_{qq}(x)+\frac{3}{2}\delta(1-x)\right]\,,$
$\displaystyle\Delta P_{qg}^{(1,0)}(x)$ $\displaystyle=$ $\displaystyle
T_{R}\left[4x-2\right]\equiv T_{R}\,\delta p_{qg}(x)\,,$ $\displaystyle\Delta
P_{gq}^{(1,0)}(x)$ $\displaystyle=$ $\displaystyle
C_{F}\left[4-2x\right]\equiv C_{F}\,\delta p_{gq}(x)\,,$ $\displaystyle\Delta
P_{gg}^{(1,0)}(x)$ $\displaystyle=$ $\displaystyle
2\,C_{A}\left[\frac{1}{(1-x)_{+}}-2x+1\right]+\beta_{0}\,\delta(1-x)\equiv
2\,C_{A}\,\delta p_{gg}(x)+\beta_{0}\,\delta(1-x)\,,$ (22)

with $\beta_{0}=\frac{11N_{C}-4n_{f}T_{R}}{3}$ and the usual plus distribution
is defined as

$\int_{0}^{1}dx\,\frac{f(x)}{(1-x)_{+}}=\int_{0}^{1}dx\,\frac{f(x)-f(1)}{1-x}\,,$
(23)

for any regular test function $f$. Therefore, setting the factors as Eq.(20),
the lowest order splitting functions in QED $P_{ij}^{(0,1)}$ result in

$\displaystyle\Delta P_{qq}^{(0,1)}(x)$ $\displaystyle=$ $\displaystyle
2\,e_{q}^{2}\left[\delta p_{qq}(x)+\frac{3}{2}\delta(1-x)\right]\,,$
$\displaystyle\Delta P_{q\gamma}^{(0,1)}(x)$ $\displaystyle=$ $\displaystyle
N_{C}\,e_{q}^{2}\,\delta p_{qg}(x)\,,$ $\displaystyle\Delta P_{\gamma
q}^{(0,1)}(x)$ $\displaystyle=$ $\displaystyle e_{q}^{2}\,\delta p_{gq}(x)\,,$
(24) $\displaystyle\Delta P_{\gamma\gamma}^{(0,1)}(x)$ $\displaystyle=$
$\displaystyle-\frac{4}{3}\sum_{f}e_{f}^{2}\,\delta(1-x)\,,$

where there is an explicit dependence on the quark electromagnetic charge.
Furthermore, the sum over fermion charges in the $\Delta
P^{(0,1)}_{\gamma\gamma}$ kernel corresponds to the definition

$\sum_{f}e_{f}^{2}=N_{C}\sum_{q}^{n_{f}}e_{q}^{2}\,+\,\sum_{l}^{n_{L}}e_{l}^{2}\,,$
(25)

with $n_{f}$ and $n_{L}$ the number of quark and lepton flavours,
respectively.

The expressions for NLO QCD corrections to the splitting functions $\Delta
P_{ij}^{(2,0)}$ can be found in Mertig and van Neerven (1996); Vogelsang
(1996a, b).

##  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

$\displaystyle a(N)=\int_{0}^{1}dx\,x^{N-1}\,a(x),$ (26)

and its inverse reads

$\displaystyle a(x)=\frac{1}{2\pi i}\int_{\textit{C}_{N}}dN\,x^{-N}\,a(N),$
(27)

here $\textit{C}_{N}$ denotes a suitable contour in the complex $N$ plane that
has an imaginary part ranging from $-\infty$ to $\infty$ 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 $<\pi/2$ towards the
negative real-$N$ axis Vogt (2005); Gluck  _et al._ (1990). 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 Vogt (2005).

The advantage of working in this space is that the convolutions appearing in
the evolution equations can be written simply as products,

$(a\otimes b)(N)=a(N)\,b(N),$ (28)

which makes it easier to manipulate the expressions and solve the
corresponding equations analytically at a given order. For 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.

###  IV.1 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

$\displaystyle\frac{df(Q)}{dt}=(a_{\mathrm{S}}P^{(1,0)}+a_{\mathrm{S}}^{2}P^{(2,0)}+aP^{(0,1)}+{\mathcal{O}}(a_{\mathrm{S}}^{3},a^{2},a_{\mathrm{S}}a))f(Q_{0}),$
(29)

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,

$f(Q)=E^{QCD}\left(Q,Q_{0}\right)E^{QED}\left(Q,Q_{0}\right)f\left(Q_{0}\right)=E^{(2,0)}\left(Q,Q_{0}\right)E^{(1,0)}\left(Q,Q_{0}\right)E^{(0,1)}\left(Q,Q_{0}\right)f\left(Q_{0}\right).$
(30)

We will make use of the well known RGE for the strong and electromagnetic
coupling constants,

$\displaystyle\frac{da_{\mathrm{S}}}{dt}$ $\displaystyle=$
$\displaystyle-(a_{\mathrm{S}}^{2}\beta_{0}+a_{\mathrm{S}}^{3}\beta_{1}+...),$
$\displaystyle\frac{da}{dt}$ $\displaystyle=$
$\displaystyle-(a^{2}\beta^{\prime}_{0}+...),$ (31)

to express the Eq.(29) in terms of $a$ and $a_{\mathrm{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

$\displaystyle\frac{dE^{QCD}\left(Q,Q_{0}\right)}{da_{S}}$ $\displaystyle=$
$\displaystyle-\frac{a_{\mathrm{S}}P^{(1,0)}+a_{\mathrm{S}}^{2}P^{(2,0)}}{a_{\mathrm{S}}^{2}(\beta_{0}+a_{\mathrm{S}}\beta_{1})}E^{QCD}\left(Q,Q_{0}\right)$
(32) $\displaystyle\simeq$
$\displaystyle-\frac{1}{a_{\mathrm{S}}\beta_{0}}\left[P^{(1,0)}+a_{\mathrm{S}}\left(P^{(2,0)}-\frac{\beta_{1}}{\beta_{0}}P^{(1,0)}\right)\right]E^{QCD}\left(Q,Q_{0}\right),$
$\displaystyle\frac{dE^{QED}\left(Q,Q_{0}\right)}{da}$ $\displaystyle=$
$\displaystyle-\frac{1}{a\beta^{\prime}_{0}}P^{(0,1)}E^{QED}\left(Q,Q_{0}\right).$

Let’s start with the solution for the QCD sector, if we remain at LO accuracy
we have

$\displaystyle\frac{dE^{(1,0)}(Q,Q_{0})}{da_{\mathrm{S}}}=-\frac{1}{a_{\mathrm{S}}\beta_{0}}P^{(1,0)}\,E^{(1,0)}(Q,Q_{0}),$
(33)

which, with the constraint $E^{LO}(Q_{0},Q_{0})=1$, we obtain the well known
exponential solution

$\displaystyle
E^{(1,0)}(Q,Q_{0})=\left(\frac{a_{\mathrm{S}}}{a_{\mathrm{S},0}}\right)^{\frac{-P^{(1,0)}}{\beta_{0}}},$
(34)

where we define $a_{\mathrm{S}}\equiv a_{\mathrm{S}}(Q)$ and
$a_{\mathrm{S},0}\equiv a_{\mathrm{S}}(Q_{0})$. For the NLO component one gets

$\displaystyle\frac{dE^{(2,0)}(Q,Q_{0})}{da_{\mathrm{S}}}=-\frac{1}{\beta_{0}}\left(P^{(2,0)}-\frac{\beta_{1}}{\beta_{0}}P^{(1,0)}\right)\,E^{(2,0)}(Q,Q_{0}),$
(35)

with the solution valid up to NLO accuracy as

$\displaystyle E^{(2,0)}(Q,Q_{0})$ $\displaystyle=$
$\displaystyle\exp\left[-\frac{a_{\mathrm{S}}-a_{\mathrm{S},0}}{\beta_{0}}\left(P^{(2,0)}-\frac{\beta_{1}}{\beta_{0}}P^{(1,0)}\right)\right]$
(36) $\displaystyle=$ $\displaystyle
1-\frac{a_{\mathrm{S}}-a_{\mathrm{S},0}}{\beta_{0}}\left(P^{(2,0)}-\frac{\beta_{1}}{\beta_{0}}P^{(1,0)}\right)+{\cal
O}(a_{\mathrm{S}}^{2}).$

By combining the LO and NLO operators, we arrive at

$\displaystyle
E^{QCD}\left(Q,Q_{0}\right)=\left(\frac{a_{\mathrm{S}}}{a_{\mathrm{S},0}}\right)^{\frac{-P^{(1,0)}}{\beta_{0}}}\left(1-\frac{a_{\mathrm{S}}-a_{\mathrm{S},0}}{\beta_{0}}\left(P^{(2,0)}-\frac{\beta_{1}}{\beta_{0}}P^{(1,0)}\right)\right)+{\cal
O}(a_{\mathrm{S}}^{2}).$ (37)

The same procedure can be applied to the QED sector as well. The resulting QED
evolution operator is given by

$\displaystyle
E^{QED}\left(Q,Q_{0}\right)=E^{(0,1)}\left(Q,Q_{0}\right)=\left(\frac{a}{a_{0}}\right)^{-\frac{P^{(0,1)}}{\beta^{\prime}}},$
(38)

where we define $a\equiv a(Q)$ and $a_{0}\equiv a(Q_{0})$. As long as all
flavors are active, all combinations evolve as non-singlets. But if one quark
flavor is not active, the corresponding $\Delta_{ij}$ combination containing
that quark term can expressed as a combination of singlet and non-singlet
terms. (e.g. at $Q=1$ GeV there is no charm, so we express
$\Delta_{uc}=\frac{1}{2}(\Sigma+\Delta_{UD})$).

###  IV.2 Singlet case

For the singlet sector, we can express the evolution equations (Eqs.(16),
(17), (18) and (19)) in a matrix form,

$\displaystyle\frac{d\bar{f^{S}}}{dt}=\mathbf{P}.\bar{f^{S}},$ (39)

where $\bar{f^{S}}$ is the vector of singlet combinations of flavors Eq.(II)
and $\mathbf{P}$ is the matrix of kernels,

$\mathbf{P}=\left(\begin{array}[]{cccc}\frac{P_{u}^{+}+P_{d}^{+}}{2}&\frac{P_{u}^{+}-P_{d}^{+}}{2}+(n_{u}-n_{d})P^{S}&2(n_{u}P_{ug}-n_{d}P_{dg})&2(n_{u}P_{u\gamma}-n_{d}P_{d\gamma})\\\
\frac{P_{u}^{+}-P_{d}^{+}}{2}&\frac{P_{u}^{+}+P_{d}^{+}}{2}+n_{f}P^{S}&2(n_{u}P_{ug}+n_{d}P_{dg})&2(n_{u}P_{u\gamma}+n_{d}P_{d\gamma})\\\
0&P_{gq}&P_{gg}&0\\\ \frac{P_{\gamma u}-P_{\gamma d}}{2}&\frac{P_{\gamma
u}+P_{\gamma d}}{2}&0&P_{\gamma\gamma}\end{array}\right),$ (40)

the kernel matrix can be expressed in a perturbative form,

$\displaystyle\mathbf{P}=a_{\mathrm{S}}\mathbf{P}^{(1,0)}+a_{\mathrm{S}}^{2}\mathbf{P}^{(2,0)}+a\mathbf{P}^{(0,1)}+\cdots.$
(41)

For QCD matrices, with no photon content as noticed in the last column/row, we
have,

$\displaystyle\mathbf{P}^{(1,0)}=\left(\begin{array}[]{cccc}P_{qq}^{(1,0)}&0&2\,n_{ud}P_{qg}^{(1,0)}&0\\\
0&P_{qq}^{(1,0)}&2\,n_{f}P_{qg}^{(1,0)}&0\\\
0&P_{gq}^{(1,0)}&P_{gg}^{(1,0)}&0\\\
0&0&0&0\end{array}\right),\,\mathbf{P}^{(2,0)}=\left(\begin{array}[]{cccc}P^{+^{(2,0)}}&n_{ud}P^{S^{(2,0)}}&2\,n_{ud}P_{qg}^{(2,0)}&0\\\
0&P_{qq}^{(2,0)}&2\,n_{f}P_{qg}^{(2,0)}&0\\\
0&P_{gq}^{(2,0)}&P_{gg}^{(2,0)}&0\\\ 0&0&0&0\end{array}\right),$ (50)

where $n_{ud}\equiv n_{u}-n_{d}$ and we use that to this order, we can
express:

$\displaystyle P_{ug}^{(i,0)}$ $\displaystyle=$ $\displaystyle
P_{dg}^{(i,0)}=P_{qg}^{(i,0)},$ $\displaystyle P_{gu}^{(i,0)}$
$\displaystyle=$ $\displaystyle P_{gd}^{(i,0)}=P_{gq}^{(i,0)},$ $\displaystyle
P^{+^{(1,0)}}_{d}$ $\displaystyle=$ $\displaystyle
P^{+^{(1,0)}}_{u}=P^{(1,0)}_{qq},$ $\displaystyle P^{+^{(2,0)}}_{u}$
$\displaystyle=$ $\displaystyle P^{+^{(2,0)}}_{d}\equiv P^{+^{(2,0)}},$
$\displaystyle P_{qq}^{(2,0)}$ $\displaystyle\equiv$ $\displaystyle
P^{+(2,0)}+n_{f}P^{S(2,0)},$

where in the first two lines $i=1,2$. The expression of the kernels
$P^{(1,0)}_{qq}$, $P^{(1,0)}_{qg}$, $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_{q\bar{q}}^{V(1,0)}=0$, and to
construct the LO matrix, we employ that $P^{S(1,0)}=0$. The NLO polarized
kernels $P_{qg}^{(2,0)}$, $P_{gq}^{(2,0)}$, $P^{(2,0)}_{gg}$, $P^{S(2,0)}$ and
$P^{+^{(2,0)}}$ can be found in Mertig and van Neerven (1996); Vogelsang
(1996a, b). For the QED sector, with no gluon content in the third column/row,
we have,

$\displaystyle\mathbf{P}^{(0,1)}=\left(\begin{array}[]{llll}\eta^{+}P_{ff}^{(0,1)}&\eta^{-}P_{ff}^{(0,1)}&0&\delta^{2}P_{f\gamma}^{(0,1)}\\\
\eta^{-}P_{ff}^{(0,1)}&\eta^{+}P_{ff}^{(0,1)}&0&e^{2}_{\Sigma}P_{f\gamma}^{(0,1)}\\\
0&0&0&0\\\ \eta^{-}P_{\gamma f}^{(0,1)}&\eta^{+}P_{\gamma
f}^{(0,1)}&0&P_{\gamma\gamma}^{(0,1)}\end{array}\right),$ (55)

where we use that to this perturbative order, we can express

$\displaystyle P_{u(d)}^{\pm{(0,1)}}$ $\displaystyle=$ $\displaystyle
e_{u(d)}^{2}P_{ff}^{(0,1)},$ (56) $\displaystyle P_{u(d)\gamma}^{(0,1)}$
$\displaystyle=$ $\displaystyle e_{u(d)}^{2}P_{f\gamma}^{(0,1)},$ (57)
$\displaystyle P_{\gamma u(d)}^{(0,1)}$ $\displaystyle=$ $\displaystyle
e_{u(d)}^{2}P_{\gamma f}^{(0,1)},$ (58)

with the definitions of $P_{ff}^{(0,1)}$,$P_{f\gamma}^{(0,1)}$,$P_{\gamma
f}^{(0,1)}$ deduced from Eq.(24) by extracting the electric charge $e_{q}^{2}$
of the corresponding kernels and taking the corresponding transform to Mellin
$N$-space. In the first line, we use that $P_{u(d)}^{\pm(0,1)}=\Delta
P_{qq}^{(0,1)}$ , which arises from the fact that $P^{V(0,1)}_{q\bar{q}}=0$ in
Eq.(6). In addition, we define

$\displaystyle\eta^{\pm}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\,(e_{u}^{2}\pm e_{d}^{2}),$ (59) $\displaystyle
e^{2}_{\Sigma}$ $\displaystyle=$ $\displaystyle
2\,(n_{u}\,e_{u}^{2}+n_{d}\,e_{d}^{2}),$ (60) $\displaystyle\delta^{2}$
$\displaystyle=$ $\displaystyle 2\,(n_{u}\,e_{u}^{2}-n_{d}\,e_{d}^{2}).$ (61)

One non-trivial problem arises from the fact that the evolution matrix in
Eq.(41) has non-commutative properties at different orders. Even the LO QCD
matrix does not conmute with the QED one. This makes more difficult to isolate
the QED and QCD solutions than in the non-singlet case. However, if mixed-
order terms ($\mathcal{O}(a\,a_{\mathrm{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.3. In
summary, we obtain matrix equations analogous to Eq.(32) that, in the case of
QCD, correspond to

$\displaystyle\frac{d\mathbf{E}^{QCD}\left(Q,Q_{0}\right)}{da_{S}}\simeq-\frac{1}{a_{\mathrm{S}}\beta_{0}}\left[\mathbf{P}^{(1,0)}+a_{\mathrm{S}}\left(\mathbf{P}^{(2,0)}-\frac{\beta_{1}}{\beta_{0}}\mathbf{P}^{(1,0)}\right)\right]\mathbf{E}^{QCD}\left(Q,Q_{0}\right).$
(62)

where the matrices $\mathbf{P}^{(1,0)}$ and $\mathbf{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 Vogt (2005); Mottaghizadeh  _et
al._ (2017). The LO solution can be expressed as a closed exponential form,

$\displaystyle\mathbf{E}_{LO}^{QCD}\left(Q,Q_{0}\right)=\left(\frac{a_{\mathrm{S}}}{a_{\mathrm{S},0}}\right)^{-\frac{\mathbf{P^{(1,0)}}}{\beta_{0}}}.$
(63)

Note that in the exponent we have the matrix $\mathbf{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,

$\displaystyle\mathbf{P}^{(1,0)}=\lambda_{1}\,\mathbf{e}_{1}+\lambda_{2}\,\mathbf{e}_{2}+\lambda_{3}\,\mathbf{e}_{3}+\lambda_{4}\,\mathbf{e}_{4},$
(64)

where $\lambda_{i}$ are given by

$\displaystyle\lambda_{1,2}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(P_{qq}^{(1,0)}+P_{gg}^{(1,0)}\pm\sqrt{(P_{qq}^{(1,0)}-P_{gg}^{(1,0)})^{2}+8\,n_{f}\,P_{qg}^{(1,0)}P_{gq}^{(1,0)}}\,\right),$
(65) $\displaystyle\lambda_{3}$ $\displaystyle=$ $\displaystyle
P_{qq}^{(1,0)},\hskip 28.45274pt\lambda_{4}=0,$

and the matrices $\mathbf{e}_{i}$ are the projectors to the subspace
corresponding to each eigenvector,

$\displaystyle\mathbf{e}_{1}$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{cccc}0&\frac{n_{ud}\,(\lambda_{2}-P_{qq}^{(1,0)})}{n_{f}\,(\lambda_{1}-\lambda_{2})}&\frac{n_{ud}\,(\lambda_{2}-P_{qq}^{(1,0)})(P_{qq}^{(1,0)}-\lambda_{1})}{n_{f}\,P_{gq}^{(1,0)}(\lambda_{1}-\lambda_{2})}&0\\\
0&\frac{\lambda_{2}-P_{qq}^{(1,0)}}{\lambda_{2}-\lambda_{1}}&\frac{(\lambda_{2}-P_{qq}^{(1,0)})(\lambda_{1}-P_{qq}^{(1,0)})}{P_{gq}^{(1,0)}(\lambda_{2}-\lambda_{1})}&0\\\
0&\frac{P_{gq}^{(1,0)}}{\lambda_{1}-\lambda_{2}}&\frac{P_{qq}^{(1,0)}-\lambda_{1}}{\lambda_{2}-\lambda_{1}}&0\\\
0&0&0&0\end{array}\right),$ (70) $\displaystyle\mathbf{e}_{2}$
$\displaystyle=$
$\displaystyle\left(\begin{array}[]{cccc}0&\frac{n_{ud}\,(\lambda_{1}-P_{qq}^{(1,0)})}{n_{f}\,(\lambda_{1}-\lambda_{2})}&\frac{n_{ud}\,(\lambda_{2}-P_{qq}^{(1,0)})(P_{qq}^{(1,0)}-\lambda_{1})}{n_{f}\,P_{gq}^{(1,0)}(\lambda_{2}-\lambda_{1})}&0\\\
0&\frac{\lambda_{1}-P_{qq}^{(1,0)}}{\lambda_{1}-\lambda_{2}}&\frac{(\lambda_{2}-P_{qq}^{(1,0)})(\lambda_{1}-P_{qq}^{(1,0)})}{P_{gq}^{(1,0)}(\lambda_{1}-\lambda_{2})}&0\\\
0&\frac{P_{gq}^{(1,0)}}{\lambda_{2}-\lambda_{1}}&\frac{P_{qq}^{(1,0)}-\lambda_{2}}{\lambda_{1}-\lambda_{2}}&0\\\
0&0&0&0\end{array}\right),$ (75) $\displaystyle\mathbf{e}_{3}$
$\displaystyle=$
$\displaystyle\left(\begin{array}[]{cccc}1&-\frac{n_{ud}}{n_{f}}&0&0\\\
0&0&0&0\\\ 0&0&0&0\\\ 0&0&0&0\end{array}\right),\hskip
28.45274pt\mathbf{e}_{4}=\left(\begin{array}[]{cccc}0&0&0&0\\\ 0&0&0&0\\\
0&0&0&0\\\ 0&0&0&1\end{array}\right),$ (84)

with a bit of algebra we arrive at the following expression for the evolution
operator at LO,

$\displaystyle\mathbf{E}_{LO}^{QCD}\left(Q,Q_{0}\right)=\mathbf{e}_{1}\left(\frac{a_{\mathrm{S}}}{a_{\mathrm{S},0}}\right)^{-\frac{\lambda_{1}}{\beta_{0}}}+\mathbf{e}_{2}\left(\frac{a_{\mathrm{S}}}{a_{\mathrm{S},0}}\right)^{-\frac{\lambda_{2}}{\beta_{0}}}+\mathbf{e}_{3}\left(\frac{a_{\mathrm{S}}}{a_{\mathrm{S},0}}\right)^{-\frac{\lambda_{3}}{\beta_{0}}}+\mathbf{e}_{4}\left(\frac{a_{\mathrm{S}}}{a_{\mathrm{S},0}}\right)^{-\frac{\lambda_{4}}{\beta_{0}}}.$
(85)

As mentioned, we can propose a general solution of Eq.(62) employing a series
expansion around the LO result, reading Vogt (2005)

$\displaystyle\mathbf{E}^{QCD}\left(Q,Q_{0}\right)=\mathbf{U}(a_{\mathrm{S}})\,\mathbf{E}_{LO}^{QCD}\left(Q,Q_{0}\right)\,\mathbf{U}^{-1}(a_{\mathrm{S},0}),\hskip
4.2679pt\textup{ with}\hskip
2.84526pt\mathbf{U}(a_{\mathrm{S}})=1+\sum_{k=1}^{\infty}a_{\mathrm{S}}^{k}\,\mathbf{U}_{k}.$
(86)

Keeping terms up to NLO order we have the truncated calculation

$\displaystyle\mathbf{E}^{QCD}\simeq\mathbf{L}+a_{\mathrm{S}}\mathbf{U}_{1}\mathbf{L}-a_{\mathrm{S},0}\mathbf{L}\mathbf{U}_{1}.$
(87)

If we replace this solution in Eq.(62) we can obtain the following relation
for the matrix $\mathbf{U}_{1}$,

$\displaystyle[\mathbf{U}_{1},\mathbf{P}^{(1,0)}]=\mathbf{R}+\mathbf{U}_{1}\textup{,
with}\hskip
8.53581pt\mathbf{R}=\left(\mathbf{P}^{(2,0)}-\frac{\beta_{1}}{\beta_{0}}\mathbf{P}^{(1,0)}\right).$
(88)

By decomposing the matrix $\mathbf{U}_{1}$ into the eigenvalue subspaces of
$\mathbf{P}^{(1,0)}$ and using the relation Eq.(88) we can find an expression
for $\mathbf{U}_{1}$ in terms of the matrix $\mathbf{R}$ and the projectors
$\mathbf{e}_{i}$. Then from Eq.(87) we obtain the evolution operator,

$\displaystyle\mathbf{E}_{QCD}\left(Q,Q_{0}\right)=\left\\{(\frac{a_{\mathrm{S}}}{a_{\mathrm{S},0}})^{-\frac{\lambda_{1}}{\beta_{0}}}\left(\mathbf{e}_{1}-\frac{\left(a_{\mathrm{S},0}-a_{\mathrm{S}}\right)}{\beta_{0}}\mathbf{e}_{1}\mathbf{R}\mathbf{e}_{1}-\frac{1}{\beta_{0}+\lambda_{2}-\lambda_{1}}\left(a_{\mathrm{S}}\left(\frac{a_{\mathrm{S}}}{a_{\mathrm{S},0}}\right)^{\left(\lambda_{2}-\lambda_{1}\right)/\beta_{0}}-a_{\mathrm{S},0}\right)\mathbf{e}_{1}\mathbf{R}\mathbf{e}_{2}\right.\right.$
$\displaystyle\left.-\frac{1}{\beta_{0}+\lambda_{3}-\lambda_{1}}\left(a_{\mathrm{S}}\left(\frac{a_{\mathrm{S}}}{a_{\mathrm{S},0}}\right)^{\left(\lambda_{3}-\lambda_{1}\right)/\beta_{0}}-a_{\mathrm{S},0}\right)\mathbf{e}_{1}\mathbf{R}\mathbf{e}_{3}-\frac{1}{\beta_{0}+\lambda_{4}-\lambda_{1}}\left(a_{\mathrm{S}}\left(\frac{a_{\mathrm{S}}}{a_{\mathrm{S},0}}\right)^{\left(\lambda_{4}-\lambda_{1}\right)/\beta_{0}}-a_{\mathrm{S},0}\right)\mathbf{e}_{1}\mathbf{R}\mathbf{e}_{4}\right)$
$\displaystyle+\left(\frac{a_{\mathrm{S}}}{a_{\mathrm{S},0}}\right)^{-\frac{\lambda_{2}}{\beta_{0}}}\left(\mathbf{e}_{2}-\frac{\left(a_{\mathrm{S},0}-a_{\mathrm{S}}\right)}{\beta_{0}}\mathbf{e}_{2}\mathbf{R}\mathbf{e}_{2}-\frac{1}{\beta_{0}+\lambda_{1}-\lambda_{2}}\left(a_{\mathrm{S}}\left(\frac{a_{\mathrm{S}}}{a_{\mathrm{S},0}}\right)^{\left(\lambda_{1}-\lambda_{2}\right)/\beta_{0}}-a_{\mathrm{S},0}\right)\mathbf{e}_{2}\mathbf{R}\mathbf{e}_{1}\right.$
$\displaystyle\left.-\frac{1}{\beta_{0}+\lambda_{3}-\lambda_{2}}\left(a_{\mathrm{S}}\left(\frac{a_{\mathrm{S}}}{a_{\mathrm{S},0}}\right)^{\left(\lambda_{3}-\lambda_{2}\right)/\beta_{0}}-a_{\mathrm{S},0}\right)\mathbf{e}_{2}\mathbf{R}\mathbf{e}_{3}-\frac{1}{\beta_{0}+\lambda_{4}-\lambda_{2}}\left(a_{\mathrm{S}}\left(\frac{a_{\mathrm{S}}}{a_{\mathrm{S},0}}\right)^{\left(\lambda_{4}-\lambda_{2}\right)/\beta_{0}}-a_{\mathrm{S},0}\right)\mathbf{e}_{2}\mathbf{R}\mathbf{e}_{4}\right)$
$\displaystyle+\left(\frac{a_{\mathrm{S}}}{a_{\mathrm{S},0}}\right)^{-\frac{\lambda_{3}}{\beta_{0}}}\left(\mathbf{e}_{3}-\frac{\left(a_{\mathrm{S},0}-a_{\mathrm{S}}\right)}{\beta_{0}}\mathbf{e}_{3}\mathbf{R}\mathbf{e}_{3}-\frac{1}{\beta_{0}+\lambda_{1}-\lambda_{3}}\left(a_{\mathrm{S}}\left(\frac{a_{\mathrm{S}}}{a_{\mathrm{S},0}}\right)^{\left(\lambda_{1}-\lambda_{3}\right)/\beta_{0}}-a_{\mathrm{S},0}\right)\mathbf{e}_{3}\mathbf{R}\mathbf{e}_{1}\right.$
$\displaystyle\left.-\frac{1}{\beta_{0}+\lambda_{2}-\lambda_{3}}\left(a_{\mathrm{S}}\left(\frac{a_{\mathrm{S}}}{a_{\mathrm{S},0}}\right)^{\left(\lambda_{2}-\lambda_{3}\right)/\beta_{0}}-a_{\mathrm{S},0}\right)\mathbf{e}_{3}\mathbf{R}\mathbf{e}_{2}-\frac{1}{\beta_{0}+\lambda_{4}-\lambda_{3}}\left(a_{\mathrm{S}}\left(\frac{a_{\mathrm{S}}}{a_{\mathrm{S},0}}\right)^{\left(\lambda_{4}-\lambda_{3}\right)/\beta_{0}}-a_{\mathrm{S},0}\right)\mathbf{e}_{3}\mathbf{R}\mathbf{e}_{4}\right)$
$\displaystyle+\left(\frac{a_{\mathrm{S}}}{a_{\mathrm{S},0}}\right)^{-\frac{\lambda_{4}}{\beta_{0}}}\left(\mathbf{e}_{4}-\frac{\left(a_{\mathrm{S},0}-a_{\mathrm{S}}\right)}{\beta_{0}}\mathbf{e}_{4}\mathbf{R}\mathbf{e}_{4}-\frac{1}{\beta_{0}+\lambda_{1}-\lambda_{4}}\left(a_{\mathrm{S}}\left(\frac{a_{\mathrm{S}}}{a_{\mathrm{S},0}}\right)^{\left(\lambda_{1}-\lambda_{4}\right)/\beta_{0}}-a_{\mathrm{S},0}\right)\mathbf{e}_{4}\mathbf{R}\mathbf{e}_{1}\right.$
$\displaystyle\left.\left.-\frac{1}{\beta_{0}+\lambda_{2}-\lambda_{4}}\left(a_{\mathrm{S}}\left(\frac{a_{\mathrm{S}}}{a_{\mathrm{S},0}}\right)^{\left(\lambda_{2}-\lambda_{4}\right)/\beta_{0}}-a_{\mathrm{S},0}\right)\mathbf{e}_{4}\mathbf{R}\mathbf{e}_{2}-\frac{1}{\beta_{0}+\lambda_{3}-\lambda_{4}}\left(a_{\mathrm{S}}\left(\frac{a_{\mathrm{S}}}{a_{\mathrm{S},0}}\right)^{\left(\lambda_{3}-\lambda_{4}\right)/\beta_{0}}-a_{\mathrm{S},0}\right)\mathbf{e}_{4}\mathbf{R}\mathbf{e}_{3}\right)\right\\}.$

Since in the QED sector we restrict our analysis to LO, the evolution equation
becomes simpler and can be expressed as follows

$\displaystyle\frac{d\mathbf{E}^{QED}\left(Q,Q_{0}\right)}{da}$
$\displaystyle=$
$\displaystyle-\frac{1}{a\beta^{\prime}_{0}}\mathbf{P}^{(0,1)}\mathbf{E}^{QED}\left(Q,Q_{0}\right).$
(89)

Again, the solution has an exponential form and we can express it in terms of
the eigenvalues of the matrix $\mathbf{P}^{(0,1)}$ ($\lambda^{\prime}_{i}$)
and the projectors to the corresponding eigenvector subspace
$\mathbf{e}^{\prime}_{1}$, $\mathbf{e}^{\prime}_{2}$,
$\mathbf{e}^{\prime}_{3}$ and $\mathbf{e}^{\prime}_{4}$,

$\displaystyle\mathbf{E}_{LO}^{QED}\left(Q,Q_{0}\right)=\left(\frac{a}{a_{0}}\right)^{-\frac{\mathbf{P}^{(0,1)}}{\beta^{\prime}_{0}}}=\mathbf{e}^{\prime}_{1}.\left(\frac{a}{a_{0}}\right)^{-\frac{\lambda_{1}^{\prime}}{\beta^{\prime}_{0}}}+\mathbf{e}^{\prime}_{2}\left(\frac{a}{a_{0}}\right)^{-\frac{\lambda_{2}^{\prime}}{\beta^{\prime}_{0}}}+\mathbf{e}^{\prime}_{3}\left(\frac{a}{a_{0}}\right)^{-\frac{\lambda_{3}^{\prime}}{\beta^{\prime}_{0}}}+\mathbf{e}^{\prime}_{4}\left(\frac{a_{\mathrm{S}}}{a_{0}}\right)^{-\frac{\lambda_{4}^{\prime}}{\beta^{\prime}_{0}}},$
(90)

with,

$\displaystyle\mathbf{e}^{\prime}_{1}$ $\displaystyle=$
$\displaystyle\frac{-(\mathbf{P}^{(0,1)})^{2}+\mathbf{P}^{(0,1)}(\lambda^{\prime}_{2}+\lambda^{\prime}_{3})-\lambda^{\prime}_{2}\lambda^{\prime}_{3}\,\mathbf{Id}_{4x4}+\mathbf{e}^{\prime}_{4}\,\lambda^{\prime}_{2}\lambda^{\prime}_{3}}{(\lambda^{\prime}_{1}-\lambda^{\prime}_{3})(\lambda^{\prime}_{1}-\lambda^{\prime}_{2})},$
(91) $\displaystyle\mathbf{e}^{\prime}_{2}$ $\displaystyle=$
$\displaystyle\frac{-(\mathbf{P}^{(0,1)})^{2}+\mathbf{P}^{(0,1)}(\lambda^{\prime}_{1}+\lambda^{\prime}_{3})-\lambda^{\prime}_{1}\lambda^{\prime}_{3}\,\mathbf{Id}_{4x4}+\mathbf{e}^{\prime}_{4}\,\lambda^{\prime}_{1}\lambda^{\prime}_{3}}{(\lambda^{\prime}_{2}-\lambda^{\prime}_{1})(\lambda^{\prime}_{2}-\lambda^{\prime}_{3})},$
(92) $\displaystyle\mathbf{e}^{\prime}_{3}$ $\displaystyle=$
$\displaystyle\frac{-(\mathbf{P}^{(0,1)})^{2}+\mathbf{P}^{(0,1)}(\lambda^{\prime}_{1}+\lambda^{\prime}_{2})-\lambda^{\prime}_{1}\lambda^{\prime}_{2}\,\mathbf{Id}_{4x4}+\mathbf{e}^{\prime}_{4}\,\lambda^{\prime}_{1}\lambda^{\prime}_{2}}{(\lambda^{\prime}_{2}-\lambda^{\prime}_{3})(\lambda^{\prime}_{3}-\lambda^{\prime}_{1})},$
(93) $\displaystyle\mathbf{e}^{\prime}_{4}$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{cccc}0&0&0&0\\\ 0&0&0&0\\\ 0&0&1&0\\\
0&0&0&0\end{array}\right),$ (98)

where the eigenvalues are $\lambda_{4}=0$ and $\lambda_{1}$, $\lambda_{2}$
,$\lambda_{3}$ can be found in the mathematica notebook EVOLUTIONQED.nb
accompanying this manuscript. In this notebook, the eigenvalues are calculated
as the three nonzero roots of the characteristic polynomial of the kernel
matrix $\mathbf{P}^{(0,1)}$.

###  IV.3 Combining QED and QCD evolution

Equipped with an operator for the evolution of QED and another for the
evolution of QCD that satisfy two separate differential equations, we can
attempt for a combination of both effects. For simplicity, we will keep the
following discussion only up to the LO order in QCD †††Therefore, an extra
$+{\cal{O}}(a_{\mathrm{S}}^{2})$ is implicit in all the expressions below.,
but the result can be expanded to NLO without difficulties. We start by
proposing the following solution for the Eq.(39),

$\displaystyle\mathbf{E}\left(Q,Q_{0}\right)=\mathbf{E}^{QCD}\left(Q,Q_{0}\right)\mathbf{E}^{QED}\left(Q,Q_{0}\right),$
(99)

where the operators $\mathbf{E}^{QCD}$ and $\mathbf{E}^{QED}$ are solutions of
the equations Eq.(62)(only the LO order) and Eq.(89) respectively. Taking the
derivative with respect to $t$, and using the evolution equations of the
coupling constants Eq.(31), we get

$\displaystyle\frac{d\mathbf{E}\left(Q,Q_{0}\right)}{dt}$ $\displaystyle=$
$\displaystyle\\{(a\mathbf{P}^{(0,1)}+a_{\mathrm{S}}\mathbf{P}^{(1,0)})$ (100)
$\displaystyle-$
$\displaystyle\frac{1}{\beta_{0}\beta^{\prime}_{0}}\left[\mathrm{ln}\left(\frac{a_{\mathrm{S}}}{a_{\mathrm{S},0}}\right)\beta^{\prime}_{0}a\right][\mathbf{P}^{(1,0)},\mathbf{P}^{(0,1)}]\\}\mathbf{E}\left(Q,Q_{0}\right)+\mathcal{O}(a_{\mathrm{S}}^{n}\,a^{m})\,,$

with $n+m=3,\,(n,m\neq 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
$\mathcal{O}(a\,a_{\mathrm{S}})$. As we commented in the section IV, the
kernel matrices of QED ($\mathbf{P}^{(0,1)}$) and QCD ($\mathbf{P}^{(1,0)}$)
do not commute, therefore, by posing in Eq.(99) the operators
$\mathbf{E}^{QCD}\,\mathbf{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 conmutator term of order
$\mathcal{O}(a\,a_{\mathrm{S}})$ that we drop. In fact, by developing the
operators in power series of the couplings (using that,
$\mathrm{ln}\left(\frac{a_{\mathrm{S}}}{a_{\mathrm{S},0}}\right)\simeq-
a_{\mathrm{S}}\,\beta_{0}\,\mathrm{ln}\left(\frac{Q^{2}}{Q^{2}_{0}}\right)$
and
$\mathrm{ln}\left(\frac{a}{a_{0}}\right)\simeq-a\,\beta^{\prime}_{0}\,\mathrm{ln}\left(\frac{Q^{2}}{Q^{2}_{0}}\right)$),
it can be shown that the solutions have the form Bertone  _et al._ (2014),

$\displaystyle\mathbf{E}^{QCD}\mathbf{E}^{QED}=\sum(\mathbf{A}\,a_{\mathrm{S}}+\mathbf{B}\,a)^{n}+\mathbf{C}\,a\,a_{\mathrm{S}},$
(101)
$\displaystyle\mathbf{E}^{QED}\mathbf{E}^{QCD}=\sum(\mathbf{A}\,a_{\mathrm{S}}+\mathbf{B}\,a)^{n}-\mathbf{C}\,a\,a_{\mathrm{S}},$

where we keep only the lowest mixed porder term and

$\displaystyle\mathbf{A}=\frac{1}{n!}\mathbf{P}^{(1,0)}\,\mathrm{ln}\left(\frac{Q^{2}}{Q^{2}_{0}}\right),\hskip
14.22636pt\mathbf{B}=\frac{1}{n!}\mathbf{P}^{(0,1)}\,\mathrm{ln}\left(\frac{Q^{2}}{Q^{2}_{0}}\right),\hskip
14.22636pt\mathbf{C}=\frac{1}{2}[\mathbf{P}^{(1,0)},\mathbf{P}^{(0,1)}]\,\mathrm{ln}^{2}\left(\frac{Q^{2}}{Q^{2}_{0}}\right).$
(102)

Then, if we propose a symmetric solution of the form Bertone  _et al._ (2014)

$\displaystyle\mathbf{\tilde{E}}=\frac{\mathbf{E}^{QED}\mathbf{E}^{QCD}+\mathbf{E}^{QCD}\mathbf{E}^{QED}}{2},$
(103)

we can eliminate the mixed order term, leaving us at a higher accuracy. If we
calculate the derivative with respect to $t$ we get,

$\displaystyle\frac{d\mathbf{\tilde{E}}}{dt}$ $\displaystyle=$
$\displaystyle\\{a\mathbf{P}^{(0,1)}+a_{\mathrm{S}}\mathbf{P}^{(1,0)}$
$\displaystyle-$
$\displaystyle\frac{1}{2\beta_{0}\beta^{\prime}_{0}}\left[\mathrm{ln}\left(\frac{a_{\mathrm{S}}}{a_{\mathrm{S},0}}\right)\beta^{\prime}_{0}a-\mathrm{ln}\left(\frac{a}{a_{0}}\right)\beta_{0}a_{\mathrm{S}}\right][\mathbf{P}^{(1,0)},\mathbf{P}^{(0,1)}]\\}\,\mathbf{\tilde{E}}+\mathcal{O}(a_{\mathrm{S}}^{n}\,a^{m}),$

with $n+m=3,\,(n,m\neq 0)$. Developing the logarithms in power series of the
couplings, we can see that the mixed order term cancels out, and obtain the
required solution,

$\displaystyle\frac{d\mathbf{\tilde{E}}}{dt}$ $\displaystyle=$
$\displaystyle(a\mathbf{P}^{(0,1)}+a_{\mathrm{S}}\mathbf{P}^{(1,0)})\mathbf{\tilde{E}}+\mathcal{O}(a_{\mathrm{S}}^{n}\,a^{m})\,.$
(105)

##  V Phenomenological results

In this section, we will apply the evolution operators to parton distributions
and study their phenomenological implication. 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 De Florian  _et al._ (2019). On the other hand, due
to the lack of experimental polarized data, we do not count with 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 Bertone  _et al._
(2018). 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 Manohar  _et al._ (2016, 2017). However,
the NNPDF3.1luxQED fit is only valid for energies higher than
$Q^{2}=2.72\,\textup{GeV}^{2}$. To obtain the distribution at the initial
scale we want ($Q_{0}^{2}=1\,\textup{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 $\Delta\gamma(x)=\gamma(x)$

  * •

B A moderate scenario, where we assume a linear polarization,
$\Delta\gamma(x)=x\,\gamma(x)$

  * •

C The null scenario, where we take the photon pPDF to vanish at the initial
scale.

Considering that the positivity constraint $|\Delta f(x)|<f(x)$ imposes a
limit on the size of the polarized density distributions Soffer (1995), 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 is roughly reproduced by the
quark densities. At last, in the scheme C, we study the case of setting
$\Delta\gamma(x)$ to zero. As stated before, we combine these three schemes
with the NLO set of polarized distributions from De Florian  _et al._ (2019)
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,

$\displaystyle\delta f=\frac{f_{withQED}-f_{noQED}}{f_{noQED}}.$ (106)

We evolve the DSSV pPDFs supplemented with the photon pPDF from
$Q_{0}^{2}=1\,\textup{GeV}^{2}$ to $Q^{2}=800\textup{GeV}^{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}_{m_{c}}=2\,\textup{GeV}^{2}$ and
$Q^{2}_{m_{b}}=20.25\,\textup{GeV}^{2}$ respectively (we do not take into
account the top quark).

Figure 1: The photon pPDF (top) at $Q^{2}=800\,\textup{GeV}^{2}$. Ratio
between polarized and unpolarized PDF (down). The colors indicate the
different schemes chosen for the photon pPDF as explained in the main text.

In the figure 1 (top) we show the photon pPDF evolved to
$Q^{2}=800\,\textup{GeV}^{2}$ for the different choices we made for the pPDF
at $Q_{0}^{2}=1\,\textup{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.

Figure 2: Relative QED corrections Eq.(106) to the pPDFs at
$Q^{2}=800\,\textup{GeV}^{2}$.

In the figure 2 we show the relative corrections for the parton distributions
of $u$-quark, $d$-quark, $\Sigma$, $\Delta_{UD}$ 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.
In the cases of $\Sigma$ and $\bar{u}$ 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.(106). 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: gray, 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 $\Delta C_{i}^{(i,j)}$ at different perturbative orders, in
Bjorken $x$-space is expressed as,

$\displaystyle g_{1}=$ $\displaystyle\frac{1}{2}\sum_{q}e_{q}^{2}\\{(\Delta
q+\Delta\bar{q})+2\,a_{\mathrm{S}}\,[\Delta C_{q}^{(1,0)}\otimes(\Delta
q+\Delta\bar{q})+\Delta C_{g}^{(1,0)}\otimes\Delta g]$
$\displaystyle+2\,a\,\left[\Delta C_{q}^{(0,1)}\otimes(\Delta
q+\Delta\bar{q})+\Delta C_{\gamma}^{(0,1)}\otimes\Delta\gamma\right]\\}.$

Let us start by recalling the well known coefficients to NLO in QCD
(Vogelsang, 1996b; Lampe and Reya, 2000, and references therein):

$\displaystyle\Delta C_{q}^{(1,0)}(x)$ $\displaystyle=$ $\displaystyle
C_{F}\left[(1+x^{2})\left[\frac{\mathrm{ln}\left(1-x\right)}{1-x}\right]_{+}-\frac{3}{2}\frac{1}{(1-x)_{+}}-\frac{1+x^{2}}{1-x}\mathrm{ln}\left(x\right)+2+x-\left(\frac{9}{2}+\frac{\pi^{2}}{3}\right)\delta(1-x)\right],$
$\displaystyle\Delta C_{g}^{(1,0)}(x)$ $\displaystyle=$ $\displaystyle
2\,T_{R}\,[(2x-1)\left(\mathrm{ln}\left(\frac{1-x}{x}\right)-1\right)+2(1-x)].$

To calculate the QED Wilson coefficients $\Delta C_{i}^{(0,1)}$, it is
necessary to adjust the color factors of the corresponding QCD coefficients,
as described in Eq.(20). Thus, we obtain

$\displaystyle\Delta C_{q}^{(0,1)}(x)$ $\displaystyle=$ $\displaystyle
e_{q}^{2}\left[(1+x^{2})\left[\frac{\mathrm{ln}\left(1-x\right)}{1-x}\right]_{+}-\frac{3}{2}\frac{1}{(1-x)_{+}}-\frac{1+x^{2}}{1-x}\mathrm{ln}\left(x\right)+2+x-\left(\frac{9}{2}+\frac{\pi^{2}}{3}\right)\delta(1-x)\right],$
$\displaystyle\Delta C_{\gamma}^{(0,1)}(x)$ $\displaystyle=$ $\displaystyle
2\,N_{C}\,e_{q}^{2}\,\left[(2x-1)\left(\mathrm{ln}\left(\frac{1-x}{x}\right)-1\right)+2(1-x)\right].$

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.(106)
(bottom). Again, the colors 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 $\Delta
C_{q}^{(0,1)}$ in Eq.(V). In order to compare the size of the contributions of
the Wilson QED coefficients, we defined the ratios,

$\displaystyle R_{q}=\frac{\Delta C_{q}^{(0,1)}\otimes(\Delta
q+\Delta\bar{q})}{\Delta C_{q}^{(0,1)}\otimes(\Delta q+\Delta\bar{q})+\Delta
C_{\gamma}^{(0,1)}\otimes\Delta\gamma}\,\hskip 14.22636pt\textup{ and
}\,\hskip 14.22636ptR_{\gamma}=1-R_{q}.\hskip 28.45274pt$ (108)

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 $\Delta C_{q}^{(0,1)}$ is
significantly higher than that of $\Delta C_{\gamma}^{(0,1)}$ .

Figure 3: Structure function $xg_{1}$ (top) and its relative QED corrections
(down) at $Q^{2}=800\,\textup{GeV}^{2}$. On the inset we show the contribution
of each Wilson coefficient, $\Delta C_{q}^{(0,1)}$ and $\Delta
C_{\gamma}^{(0,1)}$, represented by the ratios given by Eq.(108).

##  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 ${\cal O}(\alpha)$ and ${\cal O}(a_{\mathrm{S}}^{2})$. 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 to QED contributions and from the incorporation of the $\Delta
C_{q}^{(0,1)}$ coefficient in Eq.(V). We show that the QED effects affect both
parton distributions and the structure function at the percent level, the same
order of magnitude expected at NNLO in QCD.

###### Acknowledgements.

We thank Iván Pedron for discussions and useful communications. This work is
partially supported by CONICET and ANPCyT.

## References

  * Leader (2013) E. Leader, Phys. Part. Nucl. 44, 926 (2013), and references therein.
  * Nocera  _et al._ (2012) E. R. Nocera, S. Forte, G. Ridolfi, and J. Rojo, in  _20th International Workshop on Deep-Inelastic Scattering and Related Subjects_ (2012) pp. 937–942, arXiv:1206.0201 [hep-ph] .
  * Blumlein and Bottcher (2010) J. Blumlein and H. Bottcher, Nucl. Phys. B 841, 205 (2010), arXiv:1005.3113 [hep-ph] .
  * Leader  _et al._ (2010) E. Leader, A. V. Sidorov, and D. B. Stamenov, Phys. Rev. D 82, 114018 (2010), arXiv:1010.0574 [hep-ph] .
  * Hirai and Kumano (2009) M. Hirai and S. Kumano (Asymmetry Analysis), Nucl. Phys. B 813, 106 (2009), arXiv:0808.0413 [hep-ph] .
  * de Florian  _et al._ (2008) D. de Florian, R. Sassot, M. Stratmann, and W. Vogelsang, Phys. Rev. Lett. 101, 072001 (2008), arXiv:0804.0422 [hep-ph] .
  * de Florian  _et al._ (2009) D. de Florian, R. Sassot, M. Stratmann, and W. Vogelsang, Phys. Rev. D 80, 034030 (2009), arXiv:0904.3821 [hep-ph] .
  * Aschenauer  _et al._ (2015) E.-C. Aschenauer _et al._ , (2015), arXiv:1501.01220 [nucl-ex] .
  * de Florian  _et al._ (2014) D. de Florian, R. Sassot, M. Stratmann, and W. Vogelsang, Phys. Rev. Lett. 113, 012001 (2014), arXiv:1404.4293 [hep-ph] .
  * De Florian  _et al._ (2019) D. De Florian, G. A. Lucero, R. Sassot, M. Stratmann, and W. Vogelsang, Phys. Rev. D 100, 114027 (2019), arXiv:1902.10548 [hep-ph] .
  * Nocera  _et al._ (2014) E. R. Nocera, R. D. Ball, S. Forte, G. Ridolfi, and J. Rojo (NNPDF), Nucl. Phys. B 887, 276 (2014), arXiv:1406.5539 [hep-ph] .
  * Accardi  _et al._ (2016) A. Accardi _et al._ , Eur. Phys. J. A 52, 268 (2016), arXiv:1212.1701 [nucl-ex] .
  * Altarelli and Parisi (1977) G. Altarelli and G. Parisi, Nucl. Phys. B 126, 298 (1977).
  * Moch  _et al._ (2002) S. Moch, J. A. M. Vermaseren, and A. Vogt, Nucl. Phys. B 621, 413 (2002), arXiv:hep-ph/0110331 .
  * Moch  _et al._ (2004) S. Moch, J. A. M. Vermaseren, and A. Vogt, Nucl. Phys. B 688, 101 (2004), arXiv:hep-ph/0403192 .
  * Vogt  _et al._ (2004) A. Vogt, S. Moch, and J. A. M. Vermaseren, Nucl. Phys. B 691, 129 (2004), arXiv:hep-ph/0404111 .
  * Vogt  _et al._ (2006) A. Vogt, S. Moch, and J. Vermaseren, Acta Phys. Polon. B 37, 683 (2006), arXiv:hep-ph/0511112 .
  * Ball  _et al._ (2022a) R. D. Ball _et al._ (NNPDF), Eur. Phys. J. C 82, 428 (2022a), arXiv:2109.02653 [hep-ph] .
  * Hou  _et al._ (2021) T.-J. Hou _et al._ , Phys. Rev. D 103, 014013 (2021), arXiv:1912.10053 [hep-ph] .
  * Bailey  _et al._ (2021) S. Bailey, T. Cridge, L. A. Harland-Lang, A. D. Martin, and R. S. Thorne, Eur. Phys. J. C 81, 341 (2021), arXiv:2012.04684 [hep-ph] .
  * Abt  _et al._ (2022) I. Abt _et al._ (H1, ZEUS), Eur. Phys. J. C 82, 243 (2022), arXiv:2112.01120 [hep-ex] .
  * Ball  _et al._ (2022b) R. D. Ball _et al._ (PDF4LHC Working Group), J. Phys. G 49, 080501 (2022b), arXiv:2203.05506 [hep-ph] .
  * Roth and Weinzierl (2004) M. Roth and S. Weinzierl, Phys. Lett. B 590, 190 (2004), arXiv:hep-ph/0403200 .
  * Kripfganz and Perlt (1988) J. Kripfganz and H. Perlt, Z. Phys. C 41, 319 (1988).
  * Spiesberger (1995) H. Spiesberger, Phys. Rev. D 52, 4936 (1995), arXiv:hep-ph/9412286 .
  * de Florian  _et al._ (2016a) D. de Florian, G. F. R. Sborlini, and G. Rodrigo, JHEP 10, 056 (2016a), arXiv:1606.02887 [hep-ph] .
  * de Florian  _et al._ (2016b) D. de Florian, G. F. R. Sborlini, and G. Rodrigo, Eur. Phys. J. C 76, 282 (2016b), arXiv:1512.00612 [hep-ph] .
  * Bertone  _et al._ (2014) V. Bertone, S. Carrazza, and J. Rojo (APFEL), Comput. Phys. Commun. 185, 1647 (2014), arXiv:1310.1394 [hep-ph] .
  * Sadykov (2014) R. Sadykov, (2014), arXiv:1401.1133 [hep-ph] .
  * Carrazza (2015) S. Carrazza,  _Parton distribution functions with QED corrections_ , Ph.D. thesis, Milan U. (2015), arXiv:1509.00209 [hep-ph] .
  * Martin  _et al._ (2005) A. D. Martin, R. G. Roberts, W. J. Stirling, and R. S. Thorne, Eur. Phys. J. C 39, 155 (2005), arXiv:hep-ph/0411040 .
  * Ball  _et al._ (2013) R. D. Ball, V. Bertone, S. Carrazza, L. Del Debbio, S. Forte, A. Guffanti, N. P. Hartland, and J. Rojo (NNPDF), Nucl. Phys. B 877, 290 (2013), arXiv:1308.0598 [hep-ph] .
  * Harland-Lang  _et al._ (2016) L. A. Harland-Lang, V. A. Khoze, and M. G. Ryskin, Phys. Rev. D 94, 074008 (2016), arXiv:1607.04635 [hep-ph] .
  * Schmidt  _et al._ (2016) C. Schmidt, J. Pumplin, D. Stump, and C. P. Yuan, Phys. Rev. D 93, 114015 (2016), arXiv:1509.02905 [hep-ph] .
  * Manohar  _et al._ (2016) A. Manohar, P. Nason, G. P. Salam, and G. Zanderighi, Phys. Rev. Lett. 117, 242002 (2016), arXiv:1607.04266 [hep-ph] .
  * Manohar  _et al._ (2017) A. V. Manohar, P. Nason, G. P. Salam, and G. Zanderighi, JHEP 12, 046 (2017), arXiv:1708.01256 [hep-ph] .
  * Mertig and van Neerven (1996) R. Mertig and W. L. van Neerven, Z. Phys. C 70, 637 (1996), arXiv:hep-ph/9506451 .
  * Vogelsang (1996a) W. Vogelsang, Phys. Rev. D 54, 2023 (1996a), arXiv:hep-ph/9512218 .
  * Vogelsang (1996b) W. Vogelsang, Nucl. Phys. B 475, 47 (1996b), arXiv:hep-ph/9603366 .
  * Moch  _et al._ (2014) S. Moch, J. A. M. Vermaseren, and A. Vogt, Nucl. Phys. B 889, 351 (2014), arXiv:1409.5131 [hep-ph] .
  * Vogt  _et al._ (2014) A. Vogt, S. Moch, and J. A. M. Vermaseren, PoS LL2014, 040 (2014), arXiv:1405.3407 [hep-ph] .
  * Blümlein  _et al._ (2021) J. Blümlein, P. Marquard, C. Schneider, and K. Schönwald, Nucl. Phys. B 971, 115542 (2021), arXiv:2107.06267 [hep-ph] .
  * Vogt (2005) A. Vogt, Comput. Phys. Commun. 170, 65 (2005), arXiv:hep-ph/0408244 .
  * Mottaghizadeh  _et al._ (2017) M. Mottaghizadeh, F. Taghavi Shahri, and P. Eslami, Phys. Lett. B 773, 375 (2017), arXiv:1707.00108 [hep-ph] .
  * Catani  _et al._ (2004) S. Catani, D. de Florian, G. Rodrigo, and W. Vogelsang, Phys. Rev. Lett. 93, 152003 (2004), arXiv:hep-ph/0404240 .
  * Gluck  _et al._ (1990) M. Gluck, E. Reya, and A. Vogt, Z. Phys. C 48, 471 (1990).
  * Bertone  _et al._ (2018) V. Bertone, S. Carrazza, N. P. Hartland, and J. Rojo (NNPDF), SciPost Phys. 5, 008 (2018), arXiv:1712.07053 [hep-ph] .
  * Soffer (1995) J. Soffer, Phys. Rev. Lett. 74, 1292 (1995), arXiv:hep-ph/9409254 .
  * Lampe and Reya (2000) B. Lampe and E. Reya, Phys. Rept. 332, 1 (2000), arXiv:hep-ph/9810270 .