The set $\mathcal{P}_{\text{FKS}}(f_r)$ thus can be written in a generalised form with respect to that of Ref. <cit.> which is restricted to QCD interactions,

\begin{align}
\begin{split}
\mathcal{P}_{\text{FKS}}(f_r)=\left\{(i,j)\text{ : }1\le i\le n_I+n_e+n_L+n_H+n_W, \right. \\
\left. n_I+1\le j\le n_I+n_{e(L)}+n_L,\quad i\ne j,\right.\\
\left.\mathcal{R}(f_r)\rightarrow\infty\quad\text{  if  }\quad k^0_j\rightarrow0\quad\text{  or  }\quad \vec{k}_i\parallel\vec{k}_j,\right.\\
\left.(i,j)\leftrightarrow(j,i)\text{ non-redundancy conditions},\right.\\
\left. \{q_s,q_e\}\text{ correction type conditions}\right\}
\end{split}
\label{PFKSset}
\end{align}

The $(i,j)\leftrightarrow(j,i)$ non-redundancy conditions are constructed by the following considerations.
In order to not double-count contributions where by interchanging $(i,j)\leftrightarrow(j,i)$ for $i,j\ge n_I+1$ real squared matrix elements $\mathcal{R}_{\tilde{\alpha}}(f_r)$ are the same,  all pairs $(i,j)$ are hence removed from the set of FKS pairs which are redundant. This is done according to

\begin{align}
&P_j=g, P_i\ne g,\quad (j,i)\in \mathcal{P}_{\text{FKS}}\quad \Rightarrow\quad (i,j) \notin \mathcal{P}_{\text{FKS}}\quad \text{if}\quad n_I+1\leq i\\
&P_j\neq g, P_i= g,\quad (j,i)\in \mathcal{P}_{\text{FKS}}\quad \Rightarrow\quad (i,j) \notin \mathcal{P}_{\text{FKS}}\quad \text{if}\quad n_I+1\leq i<j
\end{align}

if QCD splittings are allowed and

\begin{align}
&P_j=\gamma, P_i\ne \gamma,\quad (j,i)\in \mathcal{P}_{\text{FKS}}\quad \Rightarrow\quad (i,j) \notin \mathcal{P}_{\text{FKS}}\quad \text{if}\quad n_I+1\leq i\\
&P_j\neq \gamma, P_i= \gamma,\quad (j,i)\in \mathcal{P}_{\text{FKS}}\quad \Rightarrow\quad (i,j) \notin \mathcal{P}_{\text{FKS}}\quad \text{if}\quad n_I+1\leq i<j
\end{align}

for the QED case, respectively. The pair $(g,g)$ is a special case for which $\mathcal{R}_{\tilde{\alpha}}$ is completely symmetric with respect to $(i,j)\leftrightarrow(j,i)$ if both gluons are final states. For this case a symmetry factor has to be applied. For final-state pairs $(P_i,P_j)=(f,\bar{f})$ due to $\gamma\rightarrow f\bar{f}$ splittings with $f$ denoting all charged fermions of the SM, where $(q,\bar{q})$ additionally can originate from $g\rightarrow q\bar{q}$ splittings, a similar symmetry factor as for $(g,g)$ has to be applied. Alternatively, and what is used in the following, an ordering of $i$ and $j$ according to the particles' identity, i. e.

\begin{align}
P_j=f, P_i=\bar{f},\quad (j,i)\in \mathcal{P}_{\text{FKS}}\quad \Rightarrow\quad (i,j) \notin \mathcal{P}_{\text{FKS}}\quad \text{if}\quad n_I+1\leq i
\end{align}

must be imposed.

All possible splittings allowed from SM vertices in terms of FKS pairs can be summarised in three groups with sets $G$ containing tuples $(P_i,P_j)$ with $(i,j)\in \mathcal{P}_{\text{FKS}}(f_r)$, each subdivided into two cases regarding QCD and QED splittings.

1. 	The first group comprises the emission of gluons off quarks and photons off charged fermions or massive gauge bosons, i. e.

\begin{align}
&G_1^{\text{QCD}}=\{(q,g),\text{c.-c.}\} & G_1^{\text{QED}}=\{(f,\gamma),(W^{+},\gamma),\text{c.-c.}\}
\label{GroupG1}
\end{align}

with `c.-c.' denoting the corresponding charge-conjugates of the previous tuples. In every tuple of these two sets at least one singularity exists which comes from either soft gluons or photons. For massless quarks for the QCD or massless fermions for the QED case there is an additional collinear singularity with respect to the FKS pairs $(q,g)$ or $(f,\gamma)$ and their corresponding charge-conjugates.

2. The second group contains all splittings inducing both a soft and a collinear divergence of the real squared amplitudes. Since there are no (tree-level) self-interactions of photons in the SM $(\gamma,\gamma)$ pairs do not exist, and we are left with the sets

\begin{align}
&G_2^{\text{QCD}}=\{(g,g)\} &G_2^{\text{QED}}=\emptyset\quad.
\end{align}

3. All splittings inducing $f\bar{f}$ pairs belong to the third group for which collinear singularities in $\mathcal{R}_{\tilde{\alpha}}$ occur if $m_f=0$, but never a soft singularity. For this group the sets

\begin{align}
&G_3^{\text{QCD}}=\{(q,\bar{q})\} &G_3^{\text{QED}}=\{(f,\bar{f})\}
\label{GroupG3}
\end{align}

can be defined. Note that if the third bullet point of the $\{q_s,q_e\}$ correction type conditions as stated above applies, the intersection $G_3^{\text{QCD}}\cap G_3^{\text{QED}}=\{(q,\bar{q})\}$ demands for subtraction of both, QCD and QED singularities.

The partition function $\mathcal{S}_{\tilde{\alpha}}$ introduced above is equivalent to $\mathcal{S}_{\mathcal{I}_{\tilde{\alpha}}}=\mathcal{S}_{(i,j)}$ with $\mathcal{I}_{\tilde{\alpha}}\in\mathcal{P}_{\text{FKS}}(f_r)$ since there is at most one FKS pair associated with one singular phase-space region $\tilde{\alpha}$.
According to Eq. (<ref>) the condition

\begin{align}
\sum_{\tilde{\alpha}}\mathcal{S}_{\tilde{\alpha}}(f_r)=\sum_{(i,j)\in \mathcal{P}_{\text{FKS}}(f_r)}\mathcal{S}_{(i,j)}(f_r)=1
\end{align}

has to be imposed. With this constraint $\mathcal{S}_{(i,j)}$ can be explicitly constructed. To a certain extent, there is an arbitrariness in the concrete choice of the $\mathcal{S}_{\tilde{\alpha}}$ function describing the collinear and soft limits. However, according to Ref. <cit.>, independent of this choice, this function has to fulfil the following further conditions. First of all, it is required that $\mathcal{S}_{(i,j)}$ vanishes for all singular phase-space regions of FKS pairs distinct from $(i,j)$,

\begin{align}
&\lim_{\vec{k}_k\parallel\vec{k}_l}\mathcal{S}_{(i,j)}=0\qquad \forall (k,l)\ne (i,j) \quad\text{ with }\quad(k,l)\in \mathcal{P}_{\text{FKS}}\text{ and } m_k=m_l=0\\
&\lim_{k^0_l\rightarrow 0}\mathcal{S}_{(i,j)}=0\qquad  \forall (k,l)\text{, }l\ne j\quad \text{ with }\quad P_l=g \text{ or }P_l=\gamma\text{ and }(k,l)\in \mathcal{P}_{\text{FKS}}.
\end{align}

In the collinear limit, the condition

\begin{align}
\lim_{\vec{k}_i\parallel\vec{k}_j}\mathcal{S}_{(i,j)}=h_{ij}(z_{ij}), \qquad z_{ij}=\frac{E_j}{E_j+E_i} \quad\text{ if }m_i=m_j=0
\end{align}

has to be fulfilled. $h_{ij}(z)$ is a symmetrisation factor defined in $0\le z\le 1$ with

\begin{align}
\lim_{z\rightarrow 0}h_{ij}(z)=1, \quad \lim_{z\rightarrow 1}h_{ij}(z)=0, \quad h_{ij}(z)+h_{ij}(1-z)=1,
\end{align}

which is thus a measure of the softness of the radiated particle $j$.
In soft regions of the phase-space it is required that

\begin{align}
\lim_{k^0_j\rightarrow 0}\mathcal{S}_{(i,j)}=c_{ij} \qquad \text{ if } P_j=g\text{ or }P_j=\gamma \quad \text{ and } \sum_{i,(i,j)\in \mathcal{P}_{\text{FKS}}}c_{ij}=1
\end{align}

with $c_{ij}$ defined in $0<c_{ij}\le 1$. To meet these requirements the usual definition of $\mathcal{S}_{\tilde{\alpha}}$ used for example in <cit.> is

\begin{align}
&\mathcal{S}_{(i,j)}=\frac{1}{\mathcal{D}}\frac{h_{ij}(z_{ij})}{d_{ij}}, \qquad (i,j)\in \mathcal{P}_{\text{FKS}} \label{SalphaDef}\\
&\mathcal{D}\equiv \sum_{kl}d^{-1}_{kl}h_{kl}(z_{kl}),\quad (k,l)\in \mathcal{P}_{\text{FKS}}
\end{align}

where $d_{ij}$ must vanish if $k^0_i=0$, $k^0_j=0$ or $\vec{k}_i\parallel\vec{k}_j$, such that it can be expressed as

\begin{align}
d_{ij}=2k_i\cdot k_j \frac{k^0_ik^0_j}{(k^0_i+k^0_j)^2}\quad.
\end{align}

In the soft limit of $j$, i. e. $k^0_j\rightarrow 0$, $\mathcal{D}d_{ij}$ which appears in the denominator of Eq. (<ref>) gets ill-defined. Therefore, for these regions the modified weights

\begin{align}
d_{ij}^{\text{soft}}=\frac{2k_i\cdot\hat{k}_j}{k^0_i} \qquad\text{ with }\quad\hat{k}_j=\frac{k_j}{k^0_j}
\end{align}

have to be used.
The function $h_{ij}(z_{ij})$ can be chosen with respect to the asymmetrised set $\mathcal{P}_{\text{FKS}}$ of Eq. (<ref>) as

\begin{align}
\begin{cases}
z_{ij}=\frac{k^0_j}{k^0_j+k^0_i} & \text{ if }(P_i,P_j)=(g,g)\\
1 & \text{ else}\quad.
\end{cases}
\end{align}

§.§ Subtraction of real-emission singularities

In this section the subtraction of singularities to real squared amplitudes in the infrared limits within the real subtracted contribution $d\sigma^{n+1}_{(q_s,q_e)}$ introduced in Eq. (<ref>) is defined.
In general, we note that the singularities in the FKS scheme appear in the limits $\xi\rightarrow0$ for soft radiation and $y\rightarrow1$ for collinear radiation off the final state. Both limits at the same time describes soft-collinear radiation. Due to this, real squared amplitudes  can be expressed in a regularised form, $\hat{\mathcal{R}}_{\tilde{\alpha}}$, by writing

\begin{align}
\label{regularisedRalpha}
\end{align}

If both initial state particles are potential emitters of a radiated particle $j$ these initial states have to be treated equivalently with respect to the subtraction of collinear singularities due to the origin of the emitted parton being non-reconstructable from the experimental point of view. For this reason, the FKS pairs $(i,j)$ with $1\le i\le 2$ and $P_j=g$ or $P_j=\gamma$ induce singularities in the collinear limits $y=\pm1$. Without loss of generality and by exchanging $(1-y)\rightarrow(1-y^2)$ in Eq. (<ref>) and in the following formulae all collinear singularities can be regulated as in the special case of final state radiation.

The integral over the phase-space $\Phi_{\text{rad}}^{(E_{\tilde{\alpha}})}$ of real squared matrix elements with all infrared singularities subtracted can be formulated as

\begin{align}
\int d\Phi_{\text{rad}}^{(E_{\tilde{\alpha}})}(\xi,y,\phi)\left({\mathcal{R}}_{\tilde{\alpha}}(\xi,y)- {\mathcal{R}}_{\tilde{\alpha}}(0,y)-{\mathcal{R}}_{\tilde{\alpha}}(\xi,1)+{\mathcal{R}}_{\tilde{\alpha}}(0,1)\right)
\end{align}

which in regularised form, by means of Eq. (<ref>), gives

\begin{align}
\begin{split}
%		\int \frac{d\Phi_{\text{rad}}^{(E_{\tilde{\alpha}})}}{\xi^2(1-y)}\left(\hat{\mathcal{R}}_{\tilde{\alpha}}(\xi,y)- \hat{\mathcal{R}}_{\tilde{\alpha}}(0,y)-\hat{\mathcal{R}}_{\tilde{\alpha}}(\xi,1)+\hat{\mathcal{R}}_{\tilde{\alpha}}(0,1)\right)\\
\int^{2\pi}_{0}d\phi \int_{-1}^{1}dy\int^{\xi_{\text{max}}}_{0}d\xi  \frac{\mathcal{J}^{(E_{\tilde{\alpha}})}(\xi,y,\phi)}{\xi^2(1-y)}\left(\hat{\mathcal{R}}_{\tilde{\alpha}}(\xi,y)- \hat{\mathcal{R}}_{\tilde{\alpha}}(0,y)-\hat{\mathcal{R}}_{\tilde{\alpha}}(\xi,1)+\hat{\mathcal{R}}_{\tilde{\alpha}}(0,1)\right)
\end{split}
\label{realsubtraction}
\end{align}

The integrand of this equation is thus well-defined in the collinear and soft limit and can be evaluated by numerical means[Note, that for a numerical Monte-Carlo evaluation of the integrals, $\xi$ obtains a rescaling according to $\xi=\xi_{\text{max}}(y)\tilde{\xi}$ with new integration variable $\tilde{\xi}\in [0,1]$.  According to Sec. 4.5 of Ref. <cit.>, due to this variable transformation additional terms arise with an integrand factor $\log\xi_{\text{max}}(y)$ in the soft and $\log \xi_{\text{max}}(1)$ in the soft-collinear limit. Contrary to FSR, for which $\xi_{\text{max}}(y)=1$, for the case of ISR these terms in general do not vanish.]. As explained in Sec. (<ref>) the subtraction term associated with a factorising Born process $f_{b,\tilde{\alpha}}$ in the soft and/or collinear limit is uniquely defined in the correction type per region $\tilde{\alpha}$. With  emitter $E_{\tilde{\alpha}}$ and splitting $s_{\tilde{\alpha}}$ declared in Eq. (<ref>) this can be expressed by the numbers $\mathcal{C}^{\text{QCD}},\mathcal{C}^{\text{QED}}\in \{0,1\}$ as

\begin{align}
\begin{split}
\text{QCD splitting:  }\quad	\mathcal{C}^{\text{QCD}}(s_{\tilde{\alpha}})=1,  \mathcal{C}^{\text{QED}}(s_{\tilde{\alpha}})=0 &\quad\text{if } (P_i,P_j)\in G^{\text{QCD}} \text{ and } P_{E_{\tilde{\alpha}}}\ne \gamma\\
\text{QED splitting:  }\quad \mathcal{C}^{\text{QCD}}(s_{\tilde{\alpha}})=0, \mathcal{C}^{\text{QED}}(s_{\tilde{\alpha}})=1 &\quad\text{if } (P_i,P_j)\in G^{\text{QED}} \text{ and } P_{E_{\tilde{\alpha}}}\ne g\\
&\text{with }\quad P_i,P_j\in f_r \text{ and } P_{E_{\tilde{\alpha}}}\in f_b
\end{split}
\end{align}

From the integrand of Eq. (<ref>) the subtraction terms $S^{(E_{\tilde{\alpha}},{s}_{\tilde{\alpha}})}$ introduced in Eq. (<ref>) can be extracted

\begin{align}
\begin{split}
&=	\frac{\mathcal{J}^{(E_{\tilde{\alpha}})}(\xi,y,\phi)}{\xi^2(1-y)}\left( \hat{\mathcal{R}}_{\tilde{\alpha}}(0,y)+\hat{\mathcal{R}}_{\tilde{\alpha}}(\xi,1)-\hat{\mathcal{R}}_{\tilde{\alpha}}(0,1)\right)\quad.
\end{split}
\end{align}

The soft subtraction terms are proportional to $\hat{\mathcal{R}}_{\tilde{\alpha}}(0,y)$ which can be rewritten as

\begin{align}
\begin{split}
\hat{\mathcal{R}}_{\tilde{\alpha},{(q_s,q_e)}}(0,y)=\lim_{\xi\rightarrow0}\hat{\mathcal{R}}_{\tilde{\alpha},{(q_s,q_e)}}(\xi,y)=\lim_{\xi\rightarrow0}4\pi \mathcal{S}_{\tilde{\alpha}}^{\text{soft}}\xi^2(1-y)\sum_{k,l}^{n}\frac{k_k\cdot k_l}{(k_k\cdot k_j)(k_l\cdot k_j)}\\
\left[\alpha_s\mathcal{B}^{\text{QCD}}_{kl,{(q_s-1,q_e)}}(f_{b,\tilde{\alpha}})\mathcal{C}^{\text{QCD}}(s_{\tilde{\alpha}})+\alpha\mathcal{B}^{\text{QED}}_{kl,{(q_s,q_e-1)}}(f_{b,\tilde{\alpha}})\mathcal{C}^{\text{QED}}(s_{\tilde{\alpha}})\right]
\label{softnonint}
\end{split}
\end{align}

with the partition function $\mathcal{S}_{\tilde{\alpha}}^{\text{soft}}$ depending purely on weights $d^{\text{soft}}_{ij}$ for soft radiated gluons or photons labeled by $j$.

$\mathcal{B}^{\text{QCD}}_{ij}$ represents the colour-correlated squared matrix element and $\mathcal{B}^{\text{QED}}_{ij}$ the charge-correlated, respectively. Formally, these quantities read

\begin{align}
\begin{split}
\mathcal{B}^{\text{QCD}}_{kl}=-\left(\mathcal{M}_{kl}^{(n)}\right)_{\text{QCD}}J_b= - \mathcal{N}_{B}J_b\sum_{
\text{spin, col.}}\text{Re}& \left[
\mathcal{A}^{(n)}_{c_1,\ldots,b_k,\ldots,b_l,\ldots,c_n}\right.\\ &\times \left.T^a_{b_kd_k}T^a_{b_ld_l}\mathcal{A}^{\,*(n)}_{c_1,\ldots,d_k,\ldots,d_l,\ldots,c_n} \right]\\
\mathcal{B}^{\text{QED}}_{kl}=-\left(\mathcal{M}_{kl}^{(n)}\right)_{\text{QED}}J_b= - \mathcal{N}_{B}J_b\sum_{
\text{spin, col.}}\text{Re}&\left[\mathcal{A}^{( n)}Q_kQ_l\mathcal{A}^{\,*(n)}\right]
\label{QEDchargecorrelations}
\end{split}
\end{align}

$T^a_{b_id_i}$ denotes the generator of the $SU(3)$ group in the fundamental representation and $Q_i$ the charge operator, respectively. The latter is defined as

\begin{align}
Q_i=(-1)^{s_i}e(P_i) \quad, \qquad s_i=\left\{\begin{array}{ll}
2, & \quad \text{outgoing (anti-)particle}\\
1, & \quad \text{incoming (anti-)particle}
\end{array}\right. \quad ,
\end{align}

with $e(P_i)$ the electric charge of the particle $P_i$ which can be read off from Table <ref>.

In the collinear limit of $k_i$ and $k_j$ of the FKS pair, first of all considering the case of ISR, we can define the energy fraction $z$ as

\begin{align}
&z\equiv \frac{k_{E_{\tilde{\alpha}}}^0}{{k}_i^0} &1-z=\frac{k_j^0}{{k}_{i}^0}\quad.
\label{defzandonemz}
\end{align}

with $k^{\mu}_{E_{\tilde{\alpha}}}$ denoting the momentum of the initial state emitter after radiation of a parton.
By requiring the on-shell condition of the external states before and after radiation, the momentum $k^{\mu}$ of the parton which enters the underlying Born process can be parametrised as follows.
By fixing $j$ as the radiated parton we choose $k^{\mu}=k^{\mu}_{E_{\tilde{\alpha}}}$ with $k^2<0$ and the momentum of the incoming parton $p^{\mu}=k^{\mu}_i$ and demand $\bar{k}_{E_{\tilde{\alpha}}}^2=k_j^2=k_i^2=0$ such that $k^{\mu}$ can be described as

\begin{align}
&k^{\mu}=k_{E_{\tilde{\alpha}}}^{\mu}=zp^{\mu}+k^{\mu}_{\perp}-\frac{|\vec{k}_{\perp}|^{2}}{(1-z)}\frac{\eta^{\mu}}{2p\cdot \eta}
\end{align}

with $\eta^{\mu}$ an auxiliary light-like momentum and $k^{\mu}_{\perp}$ the transverse momentum of the radiated parton with respect to the emitter particle. In summary, for these momenta the conditions

\begin{align}
k_{\perp}k_j=0\qquad k_{\perp}\eta=0\qquad \eta^2=0
\end{align}

must be satisfied. Furthermore, we define the normalised transverse momentum

\begin{align}
\hat{k}_{\perp}^{\mu}=\frac{k^{\mu}_{\perp}}{|\vec{k}_{\perp}|}\quad.
\end{align}

In the collinear limit the splitting of a massless particle can be described by the unregularised Altarelli-Parisi splitting functions <cit.> and written exclusive in spin states of the splitting parton <cit.>, i. e.

\begin{align}
&\hat{P}^{ss^{\prime}}_{{E_{\tilde{\alpha}}\rightarrow (i,j)}}(z,\hat{k}_{\perp})=\biggl\langle s\biggl\lvert\hat{P}_{{E_{\tilde{\alpha}}\rightarrow (i,j)}}(z,\hat{k}_{\perp})\bigg\rvert s^{\prime}\biggr\rangle \quad \text{with } P_{E_{\tilde{\alpha}}}=f,\bar{f}\\ &\hat{P}^{\mu\nu}_{{E_{\tilde{\alpha}}\rightarrow (i,j)}}(z,\hat{k}_{\perp})=\biggl\langle \mu\biggl\lvert\hat{P}_{{E_{\tilde{\alpha}}\rightarrow (i,j)}}(z,\hat{k}_{\perp})\bigg\rvert \nu \biggr\rangle\quad \text{with } P_{E_{\tilde{\alpha}}}=\gamma,g\quad.
\end{align}

For simplicity, we restrict ourselves to the QED case for the collinear subtraction terms in the following. The QCD terms are analogous, replacing $\alpha$ by $\alpha_s$, $\mathcal{C}^{\text{QED}}$  by $\mathcal{C}^{\text{QCD}}$ and the QED with the corresponding QCD splitting functions.
The polarised splitting functions enter the regularised real squared matrix element in the collinear limit, $\hat{\mathcal{R}}_{\tilde{\alpha}}(\xi,y_{\text{coll}})$ with $y_{\text{coll}}=\pm 1$, which can be written in a factorised form as

\begin{align}
&\hat{\mathcal{R}}_{\tilde{\alpha}}(\xi,y_{\text{coll}})=\lim_{y\rightarrow y_{\text{coll}}}\frac{8\pi \alpha}{-k^2}\xi^2(1-y^2)\hat{P}^{ss^{\prime}}_{{E_{\tilde{\alpha}}\rightarrow (i,j)}}(z,\hat{k}_{\perp})\mathcal{B}_{ss^{\prime}}(zp)\mathcal{C}^{\text{QED}}(s_{\tilde{\alpha}}) \quad \text{with }P_{E_{\tilde{\alpha}}}={f,\bar{f}} \label{regRcollf}\\
&\hat{\mathcal{R}}_{\tilde{\alpha}}(\xi,y_{\text{coll}})=\lim_{y\rightarrow y_{\text{coll}}}\frac{8\pi \alpha}{-k^2}\xi^2(1-y^2)\hat{P}^{\mu\nu}_{{E_{\tilde{\alpha}}\rightarrow (i,j)}}(z,\hat{k}_{\perp})\mathcal{B}_{\mu\nu}(zp)\mathcal{C}^{\text{QED}}(s_{\tilde{\alpha}}) \quad \text{with } P_{E_{\tilde{\alpha}}}=\gamma\label{regRcollgam}
\end{align}

where $\xi=1-z$ and $k^2$ depends on $y$ as

\begin{align}
-2p^0k^0_j(1- y)=-2(p^0)^2(1-z)(1-y) & \text{if } y\rightarrow1\\
-2p^0k^0_j(1+ y)=-2(p^0)^2(1-z)(1+y) & \text{if } y\rightarrow-1
\end{cases}
\label{k2collISR}
\end{align}

Hence, Eq. (<ref>) and (<ref>) are eventually rendered finite in the collinear limit.
The spin-correlated squared matrix elements $\mathcal{B}_{ss^{\prime}}$ and $\mathcal{B}_{\mu\nu}$ entering these equations are defined as

\begin{align}
\end{align}

Moreover, by using the definitions above, the polarised splitting functions for ISR can be expressed as

\begin{align}
\hat{P}^{ss^{\prime}}_{f\rightarrow f\gamma}(z,\hat{k}_{\perp})&=\delta^{ss^{\prime}}e(f)^2\frac{1+z^2}{1-z} \label{ftofgam}\\
\hat{P}^{ss^{\prime}}_{f\rightarrow \gamma f}(z,\hat{k}_{\perp})&=\delta^{ss^{\prime}}n_c(f)e(f)^2\left(z^2+(1-z)^2\right)\label{ftogamf}\\
\hat{P}^{\mu\nu}_{\gamma\rightarrow f f}(z,\hat{k}_{\perp})&=e(f)^2\left(-g^{\mu\nu}z+\frac{4(1-z)}{z}\hat{k}_{\perp}^{\mu}\hat{k}_{\perp}^{\nu}\right)
\label{gamtoff}
\end{align}

with $n_c(f)$ the colour degrees of freedom of $f$. The splitting functions of Eq. (<ref>) and (<ref>) by substituting $f\leftrightarrow q$ and $\gamma\leftrightarrow g$ are identical to the analogous QCD splitting functions up to the factor $e(f)^2$ which is replaced therein by the QCD Casimir operator $C_F$. The function $\hat{P}^{ss^{\prime}}_{q\rightarrow g q}$ agrees with Eq. (<ref>) if replacing $n_c(f)e(f)^2$ with $T_F$, respectively. In addition, in QCD, $g\rightarrow gg$ splittings exists and the corresponding splitting function reads

\begin{align}
\hat{P}^{\mu\nu}_{g\rightarrow g g}(z,\hat{k}_{\perp})=C_A\left(-2\left[\frac{z}{1-z}+z(1-z)\right]g^{\mu\nu}+\frac{4(1-z)}{z}\hat{k}_{\perp}^{\mu}\hat{k}_{\perp}^{\nu}\right)\quad.
\end{align}

Analogously, for the case of FSR we redefine the energy fraction $z$ as

\begin{align}
&z\equiv\frac{k^0_i}{\bar{k}^0_i} &1-z=\frac{k_j^0}{\bar{k}^0_i}\quad.
\label{varFSRz}
\end{align}

Furthermore, $k^{\mu}$ is the momentum of the massless splitting parton. With $p^{\mu}=k^{\mu}_i$ and the on-shell conditions for the external states of the Born and real process, it can be defined as

\begin{align}
k^{\mu}=k_{E_{\tilde{\alpha}}}^{\mu}=\bar{k}^{\mu}_i=\frac{p^{\mu}}{z}+k^{\mu}_{\perp}-\frac{z|\vec{k}_{\perp}|^{2}}{(1-z)}\frac{\eta^{\mu}}{2p\cdot \eta}\quad.
\end{align}

The regularised real squared matrix elements in the collinear limit of FSR, i. e. $\hat{\mathcal{R}}_{\tilde{\alpha}}(\xi,y_{\text{coll}})$ with $y_{\text{coll}}=1$, factorises as

\begin{align}
&\hat{\mathcal{R}}_{\tilde{\alpha}}(\xi,y_{\text{coll}})=\lim_{y\rightarrow y_{\text{coll}}}\frac{8\pi \alpha}{k^2}\xi^2(1-y)\hat{P}^{ss^{\prime}}_{{E_{\tilde{\alpha}}\rightarrow (i,j)}}(z,\hat{k}_{\perp})\mathcal{B}_{ss^{\prime}}\left(\frac{p}{z}\right)\mathcal{C}^{\text{QED}}(s_{\tilde{\alpha}}) \quad \text{with }P_{E_{\tilde{\alpha}}}={f,\bar{f}} \label{regRcollfFSR}\\
&\hat{\mathcal{R}}_{\tilde{\alpha}}(\xi,y_{\text{coll}})=\lim_{y\rightarrow y_{\text{coll}}}\frac{8\pi \alpha}{k^2}\xi^2(1-y)\hat{P}^{\mu\nu}_{{E_{\tilde{\alpha}}\rightarrow (i,j)}}(z,\hat{k}_{\perp})\mathcal{B}_{\mu\nu}\left(\frac{p}{z}\right)\mathcal{C}^{\text{QED}}(s_{\tilde{\alpha}}) \quad \text{with } P_{E_{\tilde{\alpha}}}=\gamma\label{regRcollgamFSR}
\end{align}

\begin{align}
\end{align}

and the final-state polarised splitting functions for QED

\begin{align}
\hat{P}^{ss^{\prime}}_{f\rightarrow f\gamma}(z,\hat{k}_{\perp})&=\delta^{ss^{\prime}}e(f)^2\frac{1+z^2}{1-z} \label{ftofgamFSR}\\
%		\hat{P}^{ss^{\prime}}_{f\rightarrow \gamma f}(z,\hat{k}_{\perp})&=e(f)^2\frac{1+z^2}{1-z}\label{ftogamfFSR}\\
\hat{P}^{\mu\nu}_{\gamma\rightarrow f f}(z,\hat{k}_{\perp})&=n_c(f)e(f)^2\left(-g^{\mu\nu}-{4z(1-z)}\hat{k}_{\perp}^{\mu}\hat{k}_{\perp}^{\nu}\right)\quad.
\label{gamtoffFSR}
\end{align}

In the same way as for ISR explained above, the analogous QCD splitting functions follow from exchanging $e(f)^2$ with $C_F$ in Eq. (<ref>)
and $n_c(f)e(f)^2$ with $T_F$ in Eq. (<ref>), respectively. The additional $g\rightarrow gg$ splitting function in the QCD case for FSR reads

\begin{align}
\hat{P}^{\mu\nu}_{g\rightarrow g g}(z,\hat{k}_{\perp})=C_A\left(-2\left[\frac{z}{1-z}+\frac{1-z}{z}\right]g^{\mu\nu}+4z{(1-z)}\hat{k}_{\perp}^{\mu}\hat{k}_{\perp}^{\nu}\right)\quad.
\end{align}

In the soft-collinear limit the regularised real squared matrix elements follow from taking the soft limit $\xi\rightarrow0$ of Eq. (<ref>) for ISR and Eq. (<ref>) for FSR. In both cases $\hat{\mathcal{R}}_{\tilde{\alpha}}(0,y_{\text{coll}})$ are rendered finite by considering the fact that $1-z$ and $k^0_j$ are proportional to $\xi$.

§.§ Non-singular tree-level $(n+1)$-body contributions

For processes with at least one jet in hadron collisions and three jets in lepton collisions, considering EW corrections real-emission processes can occur to the same coupling order as tree-level LO $(n+1)$-body processes without singular regions.
Both processes lead to same hadronic external states and thus experimentally cannot be distinguished.
Contributions of non-singular tree-level partonic processes can occur if interferences from amplitudes to different powers in $g_s$ and $e$ are not colour-forbidden. For an example we consider the process $pp\rightarrow Zj$ at $\mathcal{O}(\alpha_s\alpha)$ which at NLO EW can get an additional jet with respect to the Born process.
Protons and jets always involve massless quarks. Thus, the sub-processes $q_u\bar{q}_d\rightarrow Zq_u\bar{q}_d$ and its charge conjugates with $m_q=0$ at $\mathcal{O}(\alpha_s\alpha^2)$ contribute in the same way as all real emission processes  which are sub-processes of $pp\rightarrow Zjj$. These additional non-vanishing contributions are due to the interference of a diagram with $W^{\pm}$ exchange leading to $\mathcal{O}(e^3)$ amplitudes with a corresponding one with gluon exchange leading to $\mathcal{O}(g_s^2e)$ amplitudes, respectively. Exemplarily these diagrams are depicted in Fig. <ref>. Due to the absence of collinear or soft splittings for which the squared amplitudes of an underlying Born process can be factored out in the corresponding infrared limits there are effectively no phase-space regions for which the interference of these amplitudes diverges.

The sum over all regions $\tilde{\alpha}$ in Eq. (<ref>) implicitly contains the sum over all possible real process structures $f_r$ introduced in Sec. <ref> for which a singular structure of the squared amplitudes in the phase space is manifest. The last term of this equation however comprises all contributions from the non-singular processes at the same NLO coupling order, $\mathcal{O}(\alpha_s^{q_s}\alpha^{q_e})$, and can be expressed as a sum over the non-singular process particle structures $f_{r,\text{non-s.}}$,

\begin{align}
\mathcal{R}_{\text{non-sing},(q_s,q_e)}\left(\Phi_{n+1}\right)=\sum_{f_{r,\text{non-s.}}}\mathcal{R}_{(q_s,q_e)}\left(f_{r,\text{non-s.}},\Phi_{n+1}\right)\quad.
\end{align}

[above left=1cm of a] (i1)$q_u$;
[above right=1cm of a] (f1)$q_u$;
[below=0.7cm of a] (b);
[below left=1cm of b] (i2)$\bar{q}_d$;
[below right=1cm of b] (f2)$\bar{q}_d$;
[above right=1em of a] (ii1);
[below right=1.5em of ii1] (ii2)$Z$;

(i1) – [fermion,small] (a) – [fermion,small] (ii1) – [fermion,small] (f1),
(a) – [gluon,small] (b),
(i2) – [anti fermion,small] (b) – [anti fermion,small] (f2),
(ii1) – [photon,small] (ii2),

[above left=1cm of a] (i1)$q_u$;
[below left=1cm of a] (i2)$\bar{q}_d$;
[right=1cm of a] (b);
[above right=1cm of b] (f1)$q_u$;
[below right=1cm of b] (f2)$\bar{q}_d$;
[above right=1em of b] (ii1);
[below right=1.5em of ii1] (ii2)$Z$;

(i1) – [fermion,small] (a) – [fermion,small] (i2),
(a) – [photon] (b),
(f2) – [fermion,small] (b) – [fermion,small] (ii1) – [fermion,small] (f1),
(ii1) – [photon,small] (ii2),

Tree-level interfering diagrams with gluon exchange at $\mathcal{O}(g_s^2e)$ and with $W^{\pm}$ exchange at $\mathcal{O}(e^3)$ leading to $\mathcal{O}(\alpha_s\alpha^2)$ contributions.

§.§ Integrated subtraction terms of the $n$-body phase-space

The formalism for the FKS terms presented in this section is based on the approach following Ref. <cit.>.

The subtraction terms in integrated form $I_{(q_s,q_e)}(\bar{\Phi}_{n})$ from Eq. (<ref>) to be integrated over $n$ particle phase-space dimensions and expanded as a series in $1/\varepsilon$-poles can be achieved by making use of dimensional regularisation writing the radiated phase-space measure of Eq. (<ref>) in $d=4-2\varepsilon$ dimensions,

\begin{align}
d\Phi_{\text{rad}}=\frac{d^{3-2\varepsilon}k_{N+1}}{2k^0_{N+1}(2\pi)^{3-2\varepsilon}}=\frac{s^{1-\varepsilon}}{(4\pi)^{3-2\varepsilon}}\xi^{1-2\varepsilon}(1-y^2)^{-\varepsilon}d\xi dyd\Omega^{(2-2\varepsilon)}\quad.
\label{intRealdimreg}
\end{align}

With this phase-space measure and Eq. (<ref>) the diverging integral of the real emission amplitudes for a singular $\tilde{\alpha}$ FKS region over the $n+1$ phase-space can be rewritten as

\begin{align}
\begin{split}
\int d\Phi_{\text{rad}}^{(E_{\tilde{\alpha}})}{\mathcal{R}}_{\tilde{\alpha}}(\xi,y)=\frac{s^{1-\varepsilon}}{(4\pi)^{3-2\varepsilon}}\int d\Omega^{(2-2\varepsilon)}\int_{-1}^{1}dy(1-y)^{-1-\varepsilon}(1+y)^{-\varepsilon}\\
\int_{0}^{\xi_{\text{max}}}d\xi \xi^{-1-2\varepsilon}\hat{\mathcal{R}}_{\tilde{\alpha}}(\xi,y)
\end{split}
\label{realdimreg}
\end{align}

We define the plus-distributions

\begin{align}
\int_{-1}^{1}dy\left(\frac{g(y)}{1\mp y}\right)_{\delta_{o/i}} f(y)&=\int_{-1}^{1}dyg(y)\frac{f(y)-f(\pm 1)\Theta(\delta_{o/i}-1\pm y)}{1\mp y}\label{plusdistrydelta}\\
\int_{0}^{\xi_{\text{max}}}d\xi\left(\frac{g(\xi)}{\xi}\right)_c f(\xi)&=\int_{0}^{\xi_{\text{max}}}d\xi g(\xi)\frac{f(\xi)-f(0)\Theta(\xi_c-\xi)}{\xi}\label{plusdistrxicut}
\end{align}

with arbitrary FKS parameters $\xi_c\in[0,\xi_{\text{max}}]$ and $\delta_{o/i} \in [0,2]$. Setting the values[The parameters $\xi_c$, $\delta_o$ and $\delta_i$ are set to these values as a default in .] $\xi_c=\xi_{\text{max}}$, $\delta_o =2$ and $\delta_i =2$, results in the simplified plus-distributions,

\begin{align}
\int_{-1}^{1}dy\left(\frac{g(y)}{1-y}\right)_+ f(y)&=\int_{-1}^{1}dyg(y)\frac{f(y)-f(1)}{1-y} \label{plusdistributionsy}\\
\int_{0}^{\xi_{\text{max}}}d\xi\left(\frac{g(\xi)}{\xi}\right)_+ f(\xi)&=\int_{0}^{\xi_{\text{max}}}d\xi g(\xi)\frac{f(\xi)-f(0)}{\xi}\quad.
\label{plusdistributions}
\end{align}

For simplicity, we use these plus-distributions and leave out the dependence on the arbitrary FKS parameters $\xi_c$ and $\delta_o$ ($\delta_i$ in the case of ISR) in the following. We come back to them only for the final explicit formulas.

By using the general trick

\begin{align}
\int dx \frac{f(x)}{x^{1+\varepsilon}}=\int dx \frac{f(x)-f(0)}{x^{1+\varepsilon}}+f(0)\int dx~ {x^{-1-\varepsilon}}
\end{align}

and the simplified plus-distributions of Eqs. (<ref>) and (<ref>),
the factors $(1-y)^{-1-\varepsilon}$ and $\xi^{-1-2\varepsilon}$ in the integrand of Eq. (<ref>) expanded in ${\varepsilon}$ can be expressed as

\begin{align}
&\xi^{-1-2\varepsilon}=-\frac{1}{2\varepsilon}\delta(\xi)+\left(\frac{1}{\xi}\right)_+-2\varepsilon\left(\frac{\log\xi}{\xi}\right)_++\mathcal{O}(\varepsilon^2) \label{xiexpansion}\\
\end{align}

Inserting these identities into Eq. (<ref>) the resulting term which is free of $\delta(\xi)$ and $\delta(1-y)$ functions and for which $\mathcal{O}(\varepsilon)$ and higher contributions are dropped is

\begin{align}
\int_{0}^{2\pi}d\phi\int_{-1}^{1}dy\int^{\xi_{\text{max}}}_{0}d\xi  {\mathcal{J}^{(E_{\tilde{\alpha}})}(\xi,y,\phi)}\left(\frac{1}{\xi}\right)_+\left(\frac{1}{1-y}\right)_+\hat{\mathcal{R}}_{\tilde{\alpha}}(\xi,y)\quad.
\end{align}

By the definition of plus distributions, this term is identical to the finite real-subtracted part of the cross section given in Eq. (<ref>), such that we are left with the subtraction terms containing $\delta(\xi)$ for the soft, $\delta(1-y)$ for the collinear and $\delta(\xi)\delta(1-y)$ for the soft-collinear description of the limits.

The soft subtraction terms in integrated form can be obtained from Eq. (<ref>) by replacing $\xi^{-1-2\varepsilon}$ with the first term of Eq. (<ref>) which yields

\begin{align}
\begin{split}
I^{\text{soft}}_{\tilde{\alpha},(q_s,q_e)}(\bar{\Phi}_{n})=&-\frac{1}{2\varepsilon}\frac{s^{1-\varepsilon}}{(4\pi)^{3-2\varepsilon}}\int d\Omega^{(2-2\varepsilon)}\int_{-1}^{1}dy(1-y^2)^{-\varepsilon}\int_{0}^{\xi_{\text{max}}}d\xi \delta(\xi) \xi^2{\mathcal{R}}_{\tilde{\alpha}}(\xi,y)\\
=&-\frac{1}{2\varepsilon}\frac{s^{1-\varepsilon}}{(4\pi)^{3-2\varepsilon}}\int d\Omega^{(2-2\varepsilon)}\int_{-1}^{1}dy(1-y^2)^{-\varepsilon}\left[\lim_{\xi\rightarrow 0}\xi^2{\mathcal{R}}_{\tilde{\alpha}}(\xi,y)\right]
\end{split}
\end{align}

Performing now the sum over all regions $\tilde{\alpha}$ corresponding to one underlying Born structure $f_b$ and using the factorised form of the real squared amplitudes in the soft limit according to Eq. (<ref>) we get

\begin{align}
\begin{split}
\sum_{\tilde{\alpha}\in \mathcal{P}_{\text{FKS}}(f_b)}I^{\text{soft}}_{\tilde{\alpha},(q_s,q_e)}(\bar{\Phi}_{n})=\sum_{k,l}^{n}\left[\frac{\alpha_s}{2\pi}\mathcal{B}^{\text{QCD}}_{kl,{(q_s-1,q_e)}}(f_b,\bar{\Phi}_{n})+\frac{\alpha}{2\pi}\mathcal{B}^{\text{QED}}_{kl,{(q_s,q_e-1)}}(f_b,\bar{\Phi}_{n})\right]\\
\underbrace{\lim_{\xi\rightarrow0}\left(-\frac{2^{2\varepsilon}}{2\varepsilon}\frac{s^{1-\varepsilon}\mu_R^{2\varepsilon}}{(2\pi)^{1-2\varepsilon}}\right)\int d\Omega^{(2-2\varepsilon)}\int_{-1}^{1}dy(1-y^2)^{-\varepsilon}\frac{\xi^2}{4}\frac{k_k\cdot k_l}{(k_k\cdot k_j)(k_l\cdot k_j)}}_{=\sum_{\rho}\mathcal{E}^{(m_k,m_l)}_{kl,\rho}}.
\label{softSubtractiontermint}
\end{split}
\end{align}

With the eikonal integral factors $\mathcal{E}^{(m_k,m_l)}_{kl,\rho}$ which come along with poles $\varepsilon^{\rho}$ for $\rho=-2,-1,0$ we arrive at the form of the soft subtraction terms for which the finite part is represented by Eq. (<ref>). The divergent terms proportional to $\mathcal{E}^{(m_k,m_l)}_{kl,-1}\varepsilon^{-1}$ and $\mathcal{E}^{(m_k,m_l)}_{kl,-2}\varepsilon^{-2}$ exactly cancel with the corresponding divergent terms from virtual soft and soft-collinear emissions of massless partons. These are contained in $\mathcal{V}^{(1)}$ and $\mathcal{V}^{(2)}$ for which the regularisation scheme, here the conventional dimensional regularisation (CDR) scheme, has to be applied consistently to all terms. For the finite remainder $I^{(0)}_{(q_s,q_e)}$ the explicit form of the eikonal integral $\mathcal{E}^{(m_k,m_l)}_{kl,0}$ is given in App. <ref>, where $\mathcal{V}^{(0)}$ are the finite remnants of one-loop integrals which are provided by OLPs as for example . In addition to the choice of the regularisation scheme, the same normalisation factor $\mathcal{N}$ has to be applied to the subtraction terms as well as to the virtual one-loop matrix elements. The normalisation of the subtraction terms in  is chosen according to the BLHA standard[Note, that the convention for the normalisation factor can differ with respect to the OLP which is interfaced. Remarks on the convention of  distinct to the one of the BLHA interface and the consequential adjustments in  are given in <cit.>.] <cit.>, i. e.

\begin{align}
\mathcal{N}(\varepsilon)=\frac{(4\pi)^{\varepsilon}}{\Gamma(1-\varepsilon)}\left(\frac{\mu_R^2}{Q^2}\right)^{\varepsilon}\quad,
\label{normalisationfactor}
\end{align}

where $Q$ denotes the Ellis-Sexton scale <cit.>.

The procedure to obtain the collinear subtraction terms in integrated form follows similarly as stated above. Since collinear singularities from ISR need a special treatment due to extra counterterms coming from the PDF evolution which will be adressed in Sec. (<ref>), this section is restricted to the subtraction terms from collinear FSR. Note, that the soft-collinear singularities are already included by the terms proportional to $1/\varepsilon^2$ in the soft subtraction terms $I^{\text{soft}}_{\tilde{\alpha},(q_s,q_e)}$. Therefore the following definition is used, representing the expression of Eq. (<ref>) for which the $1/\varepsilon$ term is subtracted and $\mathcal{O}(\alpha^2)$ terms are dropped,

\begin{align}
\begin{split}
\mathcal{P}_+(\xi)=\left(\frac{1}{\xi}\right)_+-2\varepsilon\left(\frac{\log\xi}{\xi}\right)_+
\label{Pplusdistr}
\end{split}
\end{align}

In the r. h. s. of this equation the parametrisation $k_j^0=(1-z)\bar{k}^0_i$ according to Eq. (<ref>) with $j$ the radiated parton and the FSR condition $\xi_{\text{max}}=2\bar{k}_i/\sqrt{s}$ is used. Replacing $(1-y)^{-1-\varepsilon}$ with the first term of Eq. (<ref>) and $\xi^{-1-2\varepsilon}$ with $\mathcal{P}_+(\xi)$ in Eq. (<ref>) yields

\begin{align}
\begin{split}
I^{\text{coll}}_{\tilde{\alpha}}(\bar{\Phi}_{n})=\frac{s^{1-\varepsilon}}{(4\pi)^{3-2\varepsilon}}\int d\Omega^{(2-2\varepsilon)}\int_{-1}^{1}dy\int_{0}^{\xi_{\text{max}}}d\xi {\mathcal{P}_+(\xi)}\\
\times\left[-\frac{2^{-\varepsilon}}{\varepsilon}\delta(1-y)\right](1-y)\xi^2{\mathcal{R}}_{\tilde{\alpha}}(\xi,y)
\end{split}
\end{align}

which by integrating over the Born phase-space gives

\begin{align}
\begin{split}
\int d\bar{\Phi}_{n}I^{\text{coll}}_{\tilde{\alpha}}(\bar{\Phi}_{n})=&\int (2\pi)^d\delta^{(d)}\left(k_{\oplus}+k_{\ominus}-\sum_{l=n_I+1}^{N+1}k_l\right)\prod_{l\neq i,j}d\Phi_l\\
\label{collIntegral}
\end{split}
\end{align}

The real squared amplitudes in the collinear limit, i. e. the expression in brackets, can be further written in a factorised form according to Eqs. (<ref>) and (<ref>) and space-averaged due to the angular integration over $d\Omega_j^{(2-2\varepsilon)}$,

\begin{align}
\begin{split}
\biggr\langle\lim_{y\rightarrow1}(1-y){\mathcal{R}}_{\tilde{\alpha}}(\xi,y)\biggr\rangle_{\Omega}=&\frac{8\pi\mu_R^{2\varepsilon}}{2(\bar{k}^0_i)^2z(1-z)}\langle \hat{P}\rangle_{\tilde{\alpha}}(z,\varepsilon)\\
&\times \left(\alpha_s\mathcal{B}^{\text{QCD}}_{(q_s-1,q_e)}(f_{b,\tilde{\alpha}})+\alpha\mathcal{B}^{\text{QED}}_{(q_s,q_e-1)}(f_{b,\tilde{\alpha}})\right)\quad.
\label{spaceavRalpha}
\end{split}
\end{align}

Here, $\langle \hat{P}\rangle_{\tilde{\alpha}}$ denotes the space-averaged unregularised Altarelli-Parisi splitting function for which explicit formulas are given in Eqs. (<ref>) to (<ref>). Inserting Eqs. (<ref>) and (<ref>) into Eq. (<ref>), using the relations $k_i^0=z\bar{k}_i^0$ and $k_j^0=(1-z)\bar{k}^0_i$ and performing the integration over $d\Omega^{(2-2\varepsilon)}_j$ yields

\begin{align}
\begin{split}
\int d\bar{\Phi}_{n}I^{\text{coll}}_{\tilde{\alpha},(q_s,q_e)}(\bar{\Phi}_{n})=&-\frac{2^{-\varepsilon}}{\varepsilon}\frac{(4\pi)^{\varepsilon}}{\Gamma(1-\varepsilon)}\left(\frac{\mu_R}{\bar{k}^0_i}\right)^{2\varepsilon}\int d\bar{\Phi}_n\\
&\times \int_{0}^{1}dz~z^{-2\varepsilon}(1-z)\left[(1-z)^{-1-2\varepsilon}+\frac{\xi_{\text{max}}^{2\varepsilon}}{2\varepsilon}\delta(1-z)\right]\langle\hat{P}\rangle_{\tilde{\alpha}}(z,\varepsilon).
\label{intCollsub}
\end{split}
\end{align}

The prefactor on the r. h. s. of this equation can be further modified to

\begin{align}
\end{align}

such that the normalisation $\mathcal{N}(\varepsilon)$ to be used consistently for $I^{\text{soft}}$, $I^{\text{coll}}$ and $\mathcal{V}$ is factored out.
The resulting integrals corresponding to the two terms in brackets of Eq. (<ref>), which can be performed in an analytical way, are

\begin{align}
\label{collinearintegralb}
\end{align}

By evaluating the integrals in these terms, as e. g. done in <cit.>, and expanding

\begin{align}
\mathcal{T}_E=-\sum_{\tilde{\alpha}\in\mathcal{P}_{\text{FKS}}}\delta_{EE_{\tilde{\alpha}}}\frac{1}{\varepsilon}\left(\frac{2Q^2}{s}\right)^{\varepsilon}\xi_{\text{max}}^{-2\varepsilon}\left(I^{(0)}_{z,\tilde{\alpha}}+I^{(1)}_{z,\tilde{\alpha}}\right)
\end{align}

as a series in $\varepsilon$, the collinear integrated subtraction terms $I^{\text{coll}}_{\tilde{\alpha},(q_s,q_e)}$ summed over all singular regions $\tilde{\alpha}$ corresponding to one Born structure $f_b$ yields

\begin{align}
\begin{split}
\sum_{\tilde{\alpha}\in \mathcal{P}_{\text{FKS}}(f_b)}I^{\text{coll}}_{\tilde{\alpha},(q_s,q_e)}(\bar{\Phi}_{n})=\mathcal{N}(\varepsilon)\sum_{E(f_b)}\left(\frac{\alpha_s}{2\pi}\mathcal{B}^{\text{QCD}}_{(q_s-1,q_e)}(f_{b})\mathcal{T}^{\text{QCD}}_{E(f_b)}\right.\\
\left.+\frac{\alpha}{2\pi}\mathcal{B}^{\text{QED}}_{(q_s,q_e-1)}(f_{b})\mathcal{T}^{\text{QED}}_{E(f_b)}\right)\quad.
\end{split}
\end{align}

$\mathcal{T}_E$ depends on the correction type via the corresponding splitting functions and explicitly reads

\begin{align}
\mathcal{T}_{E(f_b)}^{\text{QCD}}=\frac{1}{\varepsilon}\left(-2\log \xi_{\text{max}}C^{\text{QCD}}_{E(f_b)}+\gamma^{\text{QCD}}_{E(f_b)} \right)+\mathcal{Q}^{\text{QCD}}_{E(f_b)}
\end{align}

\begin{align}
\mathcal{Q}^{\text{QCD}}_{E(f_b)}=2\log\xi_{\text{max}}\left(\log \xi_{\text{max}}-\log \frac{Q^2}{s}\right)C_{E(f_b)}^{\text{QCD}}+\left(\log \frac{Q^2}{s}-2\log\xi_{\text{max}}\right)\gamma_{E(f_b)}^{\text{QCD}}+\gamma^{\prime \text{QCD}}_{E(f_b)}
\label{Qcoll}
\end{align}

where $\mathcal{T}_{E(f_b)}^{\text{QED}}$ and $\mathcal{Q}_{E(f_b)}^{\text{QED}}$ are analogous by substituting the superscripts $\text{`QCD'}\leftrightarrow\text{`QED'}$. The Casimir operators $C_{E(f_b)}^{\text{QCD}}$ and $C_{E(f_b)}^{\text{QED}}$ and definitions $\gamma^{\text{QCD}}_{E(f_b)}$, $\gamma^{\text{QED}}_{E(f_b)}$, $\gamma^{\prime\text{QCD}}_{E(f_b)}$ and $\gamma^{\prime\text{QED}}_{E(f_b)}$ are well-defined by the particle identity $P_{E(f_b)}$ of the corresponding emitter. Their explicit values are given in App. <ref>. The term proportional to $1/\varepsilon$ in $\mathcal{T}_{E(f_b)}$ in the FKS scheme is cancelled by the corresponding terms of $\mathcal{V}^{(1)}$ describing the collinear limit of virtual emitted partons. The finite remnant of the collinear subtraction terms is finally given by the general form already introduced in Eq. (<ref>). If keeping a dependence on FKS parameters $\xi_c$ and $\delta_o$ via the plus distributions of Eq. (<ref>) and (<ref>) explicit, the expression of $\mathcal{Q}_{E(f_b)}$ changes by the substitution

\begin{align}
\frac{Q^2}{s}\rightarrow \frac{2Q^2}{s\delta_o}
\end{align}

as well as by adding the term

\begin{align}
2\left(\log \frac{2Q^2}{s\delta_o} \log \xi_{c}-\log^2\xi_c\right)C_{E(f_b)} \quad.
\end{align}

By setting the values $\delta_o=2$ and $\xi_c=1$, Eq. (<ref>) is recovered again.

§.§ Initial state remnant terms

The integrated subtraction terms for collinear ISR can be constructed in a similar way for collinear FSR, as discussed in the last section. However, different to final-state splittings, for which the collinear singularities of the real contribution are cancelled by counter-terms from the virtual one, the corresponding divergences of the initial state arise from collinearly splitting partons of the beam and not of the hard process. Thus, the cancellation of collinear initial-state singularities requires additional counterterms coming from the DGLAP evolution of the PDFs. The formalism presented here is based on the considerations of <cit.>.

If the beam particles $B_{\oplus}$ and $B_{\ominus}$ for the cross section introduced in Eq. (<ref>) are associated with partons $d$ – which is possible due to the universality of writing cross sections in terms of convolution products with PDFs – the partonic densities up to NLO in $\alpha_s$ can be written as

\begin{align}
\Gamma^{(d)}_a(x)=\delta_{ad}\delta(1-x)-\frac{\alpha_s}{2\pi}\left(\frac{1}{\bar{\varepsilon}}P^{\text{QCD}}_{ad}(x,0)-K^{\text{QCD}}_{ad}(x)\right) +\mathcal{O}(\alpha_s^2)\quad.
\label{PDFpartonic}
\end{align}

Since the expansion of the PDFs up to NLO in $\alpha$ as well as its consequences on the following formulas is completely analogous to this expansion in $\alpha_s$ for simplicity the QED case will be left out and all superscripts `QCD' dropped for the derivation of the ISR remnant.

The function $K_{ad}(x) $ depends on the renormalisation scheme

$K_{ad}(x)=0$ in case the $\overline{\text{MS}}$ scheme is chosen. $P_{ad}(x,0)$ are the regularised unpolarised Altarelli-Parisi splitting functions in four dimensions and depend on the unregularised ones $\hat{P}_{ad}$ as <cit.>

\begin{align}
\label{regSplittingfun}
\end{align}

Related to this, the following identity for the Altarelli-Parisi kernels must be obeyed,

\begin{align}
\label{casimirsplittingfun}
\end{align}

The quantities $\gamma_a$ and $C_a$ as well as their explicit values are defined in <ref>.

For partons treated as beam particles Eq. (<ref>) takes the form

\begin{align}
d\tilde{\sigma}^{(de)}=\sum_{ab}\Gamma^{(d)}_a\star \Gamma^{(e)}_b \star d{\sigma}_{ab}
\end{align}

By denoting $d{\sigma}_{ab}$ as the subtracted and $d\hat{\sigma}_{ab}$ as the unsubtracted cross section these quantities on partonic level can be perturbatively expanded as

\begin{align}
&d\sigma_{ab}=d\sigma_{ab}^{(0)}+d\sigma_{ab}^{(1)} &d\hat{\sigma}_{ab}=d\hat{\sigma}_{ab}^{(0)}+d\hat{\sigma}_{ab}^{(1)}
\label{partonicbeamcrosssec}
\end{align}

with superscript (0) for the LO and (1) for the NLO, respectively. The subtracted cross sections $d\sigma_{ab}$ to each order can be obtained by subtracting the $1/\bar{\varepsilon}$ terms according to

\begin{align}
&d\sigma_{ab}^{(1)}=d\hat{\sigma}_{ab}^{(1)}+\underbrace{\frac{\alpha_s}{2\pi}\sum_d\left(\frac{1}{\bar{\varepsilon}}P_{da}-K_{da}\right)\star d\hat{\sigma}_{db}^{(0)}}_{d\sigma_{ab}^{(cnt,\oplus)}}+\underbrace{\frac{\alpha_s}{2\pi}\sum_e\left(\frac{1}{\bar{\varepsilon}}P_{eb}-K_{eb}\right)\star d\hat{\sigma}_{ae}^{(0)}}_{d\sigma_{ab}^{(cnt,\ominus)}}
\label{shortdistcross1}
\end{align}

where the counterterms $d\sigma_{ab}^{(cnt,\oplus)}$ and $d\sigma_{ab}^{(cnt,\ominus)}$ can be extracted. The quantity $d\hat{\sigma}_{ab}^{(0)}$ is associated with the spin- and colour-averaged Born squared matrix elements summed up for a specific initial state $\{ab\}$, i. e.

\begin{align}
\end{align}

The final aim is to compute the component $d\sigma_{ab}^{(1)}$ which contains the remnant of the collinear initial-state subtraction terms in integrated form and is introduced as $d\bar{\sigma}^{n+1}_{ab}$ with the notion `degenerate $n+1$ phase-space contribution' or `DGLAP remnant' in Sec. (<ref>). The unsubtracted cross section $d\hat{\sigma}_{ab}^{(1)}$ can be constructed in the same way as the final state collinear subtraction terms $\sum_{\tilde{\alpha}}I^{\text{coll}}_{\tilde{\alpha}}$ in the previous section. For this reason, analogously to Eq. (<ref>), we start by expanding

\begin{align}
(1-y^2)^{-1-\varepsilon}=-\frac{4^{-\varepsilon}}{2\varepsilon}\left[\delta(1-y)+\delta(1+y)\right] + \frac{1}{2}\left[\left(\frac{1}{1-y}\right)_++\left(\frac{1}{1+y}\right)_+\right]+\mathcal{O}(\varepsilon).
\label{yISexpansion}
\end{align}

Using $\mathcal{P}_+(\xi)$ defined in Eq. (<ref>) and one of the terms of Eq. (<ref>) containing a $\delta$ function, the integrated subtraction terms in the collinear limit $y=+1$ for emitter ${\oplus}$ and $y=-1$ for emitter ${\ominus}$, respectively, can be written as

\begin{align}
\begin{split}
I^{\circled{\scriptsize$\pm$}\text{ coll}}_{\tilde{\alpha}}(\bar{\Phi}_{n})=\frac{s^{1-\varepsilon}}{(4\pi)^{3-2\varepsilon}}\int d\Omega^{(2-2\varepsilon)}\int_{-1}^{1}dy\int_{0}^{\xi_{\text{max}}}d\xi {\mathcal{P}_+(\xi)}\\
\times\left[-\frac{4^{-\varepsilon}}{2\varepsilon}\delta(1\mp y)\right](1-y^2)\xi^2{\mathcal{R}}_{\tilde{\alpha}}(\xi,y)\quad.
\end{split}
\label{IcollIS}
\end{align}

Due to the angular integration measure the regularised real squared amplitudes in the collinear limit $y=+1$ can be space-averaged. By means of Eqs. (<ref>) to (<ref>) and identifying $p^0=\sqrt{s}/2$ and $k_j^0=\xi\sqrt{s}/2$ it takes the factorised form

\begin{align}
\begin{split}
\biggr\langle\lim_{y\rightarrow1}(1-y^2)\xi^2{\mathcal{R}}_{\tilde{\alpha}}(\xi,y)\biggr\rangle_{\Omega}={8\pi\alpha_s\mu_R^{2\varepsilon}}\frac{4}{s}\xi\langle \hat{P}\rangle_{\tilde{\alpha}}(z,\varepsilon)\mathcal{B}(f_{b,\tilde{\alpha}})\quad.
\label{spaceavRalphaIS}
\end{split}
\end{align}

with the unpolarised unregularised initial state splitting functions $\langle \hat{P}\rangle$ defined in eqs. (<ref>) to (<ref>). By integrating over the angular phase-space measure $d\Omega^{(2-2\varepsilon)}$ we get the factor

\begin{align}
\int d\Omega^{(2-2\varepsilon)}=\frac{2\pi^{1-\varepsilon}}{\Gamma(1-\varepsilon)}\quad.
\label{angularmeasure}
\end{align}

Inserting Eq. (<ref>) and (<ref>) into Eq. (<ref>) with the substitution $z=1-\xi$, expanding the $\varepsilon$-dependent factor as

\begin{align}
\frac{1}{\varepsilon}\frac{(4\pi)^{\varepsilon}}{\Gamma(1-\varepsilon)}\left(\frac{\mu_R^2}{s}\right)^{\varepsilon}=\frac{1}{\varepsilon}-\gamma_E+\log 4\pi+\log\frac{\mu_R^2}{s} +\mathcal{O}(\varepsilon)
\end{align}

and summing over all $\tilde{\alpha}$ regions with emitter $E_{\tilde{\alpha}}=\oplus$ leading to one initial state $\{ab\}$ yields the collinear subtraction term of the form

\begin{align}
\begin{split}
\sum_{\tilde{\alpha}}I^{\oplus,\text{coll}}_{\tilde{\alpha},ab}=I^{\oplus,\text{coll}}_{ab}=-\sum_{d}\int_{0}^{\xi_{\text{max}}}d\xi~ \left[\frac{1}{\bar{\varepsilon}}+\log\frac{\mu_R^2}{s}\right]{\mathcal{P}_+(\xi)}\xi\langle \hat{P}_{da}\rangle(1-\xi,\varepsilon)\\
\times \frac{\alpha_s}{2\pi}\mathcal{B}_{db}(1-\xi,\bar{\Phi}_n)
\end{split}
\end{align}

using the definition of the $\overline{\text{MS}}$-pole

\begin{align}
\frac{1}{\bar{\varepsilon}}=\frac{1}{\varepsilon}-\gamma_E+\log 4\pi\quad.
\end{align}

Inserting Taylor-expanded splitting functions $\langle \hat{P}_{da}\rangle(1-\xi,\varepsilon)$ for $\varepsilon\rightarrow 0$ and dropping all terms of $\mathcal{O}(\varepsilon)$ results in

\begin{align}
\begin{split}
I^{\oplus,\text{coll}}_{ab}=-\sum_{d}\int_{0}^{\xi_{\text{max}}}d\xi~\xi\left[\frac{1}{\bar{\varepsilon}}\left(\frac{1}{\xi}\right)_+\langle \hat{P}_{da}\rangle(1-\xi,0)+\left(\frac{1}{\xi}\right)_+\langle \hat{P}_{da}\rangle(1-\xi,0)\log \frac{\mu_R^2}{s}\right.\\
\left.+\left(\frac{1}{\xi}\right)_+\frac{d\langle \hat{P}_{da}\rangle}{d\varepsilon}(1-\xi,\varepsilon)\biggr\rvert_{\varepsilon=0}-2\left(\frac{\log\xi}{\xi}\right)_+\langle \hat{P}_{da}\rangle(1-\xi,0) \right]\\
\times\frac{\alpha_s}{2\pi}\mathcal{B}_{db}(1-\xi,\bar{\Phi}_n)
\end{split}
\label{Icollab}
\end{align}

Following Eqs. (<ref>) and (<ref>), related to the construction of the phase-space for ISR with massless emitters in Sec. <ref>, we observe that the maximal energy which can be radiated off one initial-state emitter in the strict collinear limit is constrained to $\xi_{\text{max}}=1-x^{\circled{\scriptsize$\pm$}}$. In the same way, this constraint leads to a lower bound for the PDF rescaling variable $z^{\circled{\scriptsize$\pm$}}$, i. e.

\begin{align}
\label{zrescalingrestriction}
x^{\circled{\scriptsize$\pm$}}\le z^{\circled{\scriptsize$\pm$}}\le 1\quad.
\end{align}

Due to this, the integral from the convolution product in the counterterm $d\sigma_{ab}^{(cnt,\oplus)}$, which formally reads

\begin{align}
d\sigma_{ab}^{(cnt,\oplus)}\propto \int^{1}_{0}dz~ \left(\frac{1}{\bar{\varepsilon}}P_{da}(z,0)-K_{da}(z)\right) \mathcal{B}_{db}(z,\bar{\Phi}_n)\quad,
\end{align}

has to be transformed in an integration over  $\xi=1-z^{\oplus}$.
Using Eq. (<ref>) expressed in terms of $\langle \hat{P}_{da}\rangle$ with the transformation of the plus distribution for $z=1-\xi$ in the collinear limit

\begin{align}
\left(\frac{1}{1-z}\right)_+=\left(\frac{1}{\xi}\right)_++\delta(\xi)\log \xi_{\text{max}}
\end{align}

as well as the identity of Eq. (<ref>) we get

\begin{align}
\begin{split}
d\sigma_{ab}^{(cnt,\oplus)}=\sum_{d}\int_{0}^{1-x^{\oplus}}d\xi \left[\frac{1}{\bar{\varepsilon}}\left(\frac{1}{\xi}\right)_+\langle \hat{P}_{da}\rangle(1-\xi,0)+\frac{1}{\bar{\varepsilon}}(\gamma_d+2C_d\log \xi_{\text{max}})\delta_{ad}\delta(\xi)\right.\\
\left.-K_{da}(1-\xi)\right]
\times \frac{\alpha_s}{2\pi}\mathcal{B}_{db}(1-\xi,\bar{\Phi}_n).
\end{split}
\label{counterterm+}
\end{align}

Since the subtracted short distance cross section $d\sigma^{(1)}_{ab}$ according to Eq. (<ref>) depends on the terms $I^{\oplus,\text{coll}}_{ab}$ and $d\sigma_{ab}^{(cnt,\oplus)}$as

\begin{align}
d\sigma^{(1)}_{ab}\propto I^{\oplus,\text{coll}}_{ab}+d\sigma_{ab}^{(cnt,\oplus)}
\end{align}

we extract the initial-state remnant for emitter $B_{\oplus}=a$

\begin{align}
\begin{split}
d\sigma_{ab}^{\oplus \text{ISR} }=&I^{\oplus,\text{coll}}_{ab}+d\sigma_{ab}^{(cnt,\oplus)}\\
=&\frac{\alpha_s}{2\pi}\frac{1}{\bar{\varepsilon}}(\gamma_a+2C_a\log \xi_{\text{max}})\mathcal{B}_{ab}(\bar{\Phi}_n)\\%+\underbrace{\frac{\alpha_s}{2\pi}\sum_{d}\mathcal{K}^{\text{QCD}}_{da}\star \mathcal{B}^{\text{QCD}}_{db}}_{=d\bar{\sigma}^{n+1,\oplus}_{ab}}
&+\frac{\alpha_s}{2\pi}\sum_{d}\int_{0}^{1-x^{\oplus}}d\xi~\xi\left(\langle \hat{P}_{da}\rangle(1-\xi,0)\left[\left(\frac{1}{\xi}\right)_+\log \frac{s}{\mu_R^2}+2\left(\frac{\log\xi}{\xi}\right)_+\right]\right. \\
&\left.- \left(\frac{1}{\xi}\right)_+\frac{d\langle \hat{P}_{da}\rangle}{d\varepsilon}(1-\xi,\varepsilon)\biggr\rvert_{\varepsilon=0}-K_{da}(1-\xi)\right)\mathcal{B}_{db}(1-\xi,\bar{\Phi}_n)
\end{split}
\label{collinearremnantISR}
\end{align}

with the last term, free of $1/\bar{\varepsilon}$-poles, identified as

\begin{align}
\begin{split}
d\bar{\sigma}^{n+1,\oplus}_{ab}=\frac{\alpha_s}{2\pi}\sum_{d}\mathcal{K}_{da}\star \mathcal{B}_{db}.
\end{split}
\label{degenerate+}
\end{align}

Hence, the $1/\bar{\varepsilon}$ term of $I^{\oplus,\text{coll}}_{ab}$ in Eq. (<ref>) is cancelled by the corresponding one of $d\sigma_{ab}^{(cnt,\oplus)}$ in Eq. (<ref>). Eventually, the pole term of Eq. (<ref>) which comes from the flavour-diagonal part of the regularised splitting functions $P_{ab}$ is cancelled exactly by the corresponding term of the virtual part $\mathcal{V}$ in the soft limit.
Using the definition of the $\overline{\text{MS}}$-scheme pole in the exact form and an additional dependence on the arbitrary Ellis-Sexton scale $Q$ the pole $1/\bar{\varepsilon}$ can be replaced by

\begin{align}
\frac{1}{\bar{\varepsilon}}=\frac{(4\pi)^{\varepsilon}}{\Gamma(1-\varepsilon)}\left(\frac{\mu_R^2}{Q^2}\right)^{\varepsilon}\frac{1}{\varepsilon}-\log \frac{\mu_R^2}{Q^2}=\mathcal{N}(\varepsilon)\frac{1}{\varepsilon}-\log \frac{\mu_R^2}{Q^2}\quad.
\end{align}

Due to this, the finite part of the first term of Eq. (<ref>) can be identified as a contribution to the remainders of the soft subtraction terms,

\begin{align}
d\sigma_{ab,S}^{\oplus\text{ISR remn.}}=- \frac{\alpha_s}{2\pi}(\gamma_a+2C_a\log \xi_{\text{max}})\mathcal{B}_{ab}(\bar{\Phi}_n)\log \frac{\mu_R^2}{Q^2}
\end{align}

Varying the scale $Q$ where $\mu_R\ne Q$ should not change the complete NLO result and thus yields a sanity check on the implemented FKS scheme.

The initial-state remnant for the emitter $B_{\ominus}=b$ is almost exactly identical to $d\sigma_{\text{ISR}}^{\oplus }$
and analogously to Eq. (<ref>) reads

\begin{align}
\begin{split}
d\sigma_{ab}^{\ominus \text{ISR} }
=&\frac{\alpha_s}{2\pi}\frac{1}{\bar{\varepsilon}}(\gamma_b+2C_b\xi_{\text{max}})\mathcal{B}_{ab}(\bar{\Phi}_n)\\%+\underbrace{\frac{\alpha_s}{2\pi}\sum_{d}\mathcal{K}^{\text{QCD}}_{da}\star \mathcal{B}^{\text{QCD}}_{db}}_{=d\bar{\sigma}^{n+1,\oplus}_{ab}}
&+\frac{\alpha_s}{2\pi}\sum_{d}\int_{0}^{1-x^{\ominus}}d\xi~\xi\left(\langle \hat{P}_{db}\rangle(1-\xi,0)\left[\left(\frac{1}{\xi}\right)_+\log \frac{s}{\mu_R^2}+2\left(\frac{\log\xi}{\xi}\right)_+\right]\right. \\
&\left.- \left(\frac{1}{\xi}\right)_+\frac{d\langle \hat{P}_{db}\rangle}{d\varepsilon}(1-\xi,\varepsilon)\biggr\rvert_{\varepsilon=0}-K_{db}(1-\xi)\right)\mathcal{B}_{ad}(1-\xi,\bar{\Phi}_n)
\end{split}
\end{align}

Furthermore, accounting for additional QED corrections, if required by the coupling power conditions, the analogous terms have to be added up. For these the definitions $\gamma_a^{\text{QED}}$, the underlying Born $\mathcal{B}^{\text{QED}}_{ab}$ and expressions $\mathcal{K}^{\text{QED}}_{da}$ containing the QED splitting functions must be used.

Summarising the finite part of the initial-state remnant, summing over all flavour combinations $\{ab\}$ which are possible from the initial state, we have

\begin{align}
d\sigma^{\text{ISR remn.} }=d\sigma_{S}^{\text{ISR remn.}}+d\bar{\sigma}^{n+1}=\sum_{a,b}\left(d\sigma_{ab,S}^{\text{ISR remn.}}+d\bar{\sigma}^{n+1}_{ab}\right)
\end{align}

with the soft ISR remnant

\begin{align}
\begin{split}
d\sigma_{S}^{\text{ISR remn.}}=-\sum_{a,b}\left[\frac{\alpha_s}{2\pi}\left(\gamma^{\text{QCD}}_a+\gamma^{\text{QCD}}_b+ 2\xi_{\text{max}}(C_a^{\text{QCD}}+C_b^{\text{QCD}})\right)\mathcal{B}^{\text{QCD}}_{ab}\right.\\
\left.+\frac{\alpha}{2\pi}\left(\gamma^{\text{QED}}_a+\gamma^{\text{QED}}_b+2\xi_{\text{max}}(C_a^{\text{QED}}+C_b^{\text{QED}})\right)\mathcal{B}^{\text{QED}}_{ab}\right]\log \frac{\mu_R^2}{Q^2}
\end{split}
\end{align}

and the $n+1$ degenerate contribution

\begin{align}
\begin{split}
=&\sum_{a,b}\sum_{d}\left[\frac{\alpha_s}{2\pi}\left(\mathcal{K}^{\text{QCD}}_{da}\star \mathcal{B}^{\text{QCD}}_{db}+\mathcal{K}^{\text{QCD}}_{db}\star \mathcal{B}^{\text{QCD}}_{ad}\right)\right.\\
&\left.+\frac{\alpha}{2\pi}\left(\mathcal{K}^{\text{QED}}_{da}\star \mathcal{B}^{\text{QED}}_{db}+\mathcal{K}^{\text{QED}}_{db}\star \mathcal{B}^{\text{QED}}_{ad}\right)\right]\quad.
\end{split}
\end{align}

The explicit form of the ISR remnant depending on the set of arbitrary FKS parameters $\xi_c$ and $\delta_i$, the ISR counterpart to $\delta_o$, is provided by replacing the plus distributions by the modified ones of Eq. (<ref>) and (<ref>), substituting

\begin{align*}
\frac{s}{\mu_R^2} \rightarrow\frac{s\delta_i}{2\mu^2_R}
\end{align*}

and by adding the term

\begin{align*}
\end{align*}

CHAPTER: EW PRECISION CALCULATIONS FOR LEPTON COLLISIONS – COLLINEAR ISR

In order to describe SM precision physics at lepton colliders, i. e. future colliders as for example ILC or CLIC, generally the same methods and schemes for NLO computations are applicable as for hadron collisions. However, due to the colliding leptons being colourless states of QED interactions beam effects are described fundamentally different than those of QCD initial states. These differences are manifest by the distinct perturbative nature of both of these SM theories.
In QED, low-energy external leptons and photons which are identified as on-shell particles are fully describable by perturbative expansions due to the abelian group structure of the theory. In this section the two concepts of higher order computations treating the leptonic initial-states either in the massive or massless approach will be presented. Since these concepts essentially are equivalent for all leptons which differ only in numerical values of their masses we refer to all of these by using notions related to electrons in the following.

In Sec. <ref> it will be discussed how a resummation of large logarithms coming from ISR effects in contrast to pure fixed order corrections in $\alpha$ affects cross sections of lepton collision processes.
The universal treatment of collinear ISR effects for MC integrated results by embedding electron PDFs into cross section computations is challenged by numerical issues of the evaluation and structure of the PDFs. Methods to solve these are suggested in Sec. <ref>.

§ DESCRIPTION OF COLLINEAR RADIATION OFF THE LEPTONIC INITIAL STATE

Adding electroweak corrections to cross sections of $e^+e^-$ processes naively is pictured as considering matrix elements for Feynman diagrams with loops or photon radiation leading to contributions of the desired higher order in $\alpha$, which are then integrated over the phase space. However, due to collinear photon emission off the electron or positron, logarithms of the ratio $Q^2/m^2$, with $Q$ the energy scale and $m$ the electron mass, lead to large terms such that the perturbation series in $\alpha$ would be spoiled if $Q^2\gg m^2$. Analogously to QCD, a way to solve this is through the structure-function approach: By means of the DGLAP evolution equations <cit.> the logarithms can be absorbed into universal functions, the structure functions (or PDFs), which can be factored out due to factorisation theorem. Whereas for QCD the initial conditions for the PDFs are non-perturbative functions and determined only by fits to experimental data, in QED these can be calculated perturbatively.

At LO the initial conditions of the electron PDFs – due to the description of electrons as on-shell elementary particles in the low energy limit – are trivially given by a delta Dirac function $\delta(1-z)$, with $z$ denoting the energy fraction which is radiated away. This description however is accurate only by the  caveat that logarithms $\log Q^2/m^2$ appearing at higher orders in $\alpha$ of the cross sections are sufficiently suppressed by the coupling, i. e. at low scales $Q$.
For large energy scales $Q$ the Dirac delta functional behaviour of the LO PDFs is smoothened to a peak in the vicinity of $z=1$ with a tail reaching down to $z=0$ if the universal scale-dependency of terms of a perturbation series is absorbed into the PDFs which is understood in the following sense.
A class of terms containing logarithms $\log Q^2/m^2$ appearing at each order $k$ of the cross section perturbatively expanded in $\alpha$, e. g. all terms at the leading-logarithmic (LL) level,

\begin{align}
\alpha^k \log^k \frac{Q^2}{m^2}C_{\text{LL}}^{(k)}\quad,
\label{LLres}
\end{align}

at the next-to-leading logarithmic (NLL) level,

\begin{align}
\alpha^k \log^k \frac{Q^2}{m^2}C_{\text{LL}}^{(k)}+\alpha^k \log^{k-1} \frac{Q^2}{m^2}C_{\text{NLL}}^{(k-1)}\quad,
\label{NLLres}
\end{align}

and so forth, can be resummed by factorising them from the matrix elements at an arbitrary scale $\mu$. In this way the matrix elements keep a dependence on logarithms $\log Q^2/\mu^2$, whereby the factorised sum rewritten as an exponential factor contains the electron mass dependency in terms of $\log \mu^2/m^2$ and is thus included into the PDF definition of the electron. This is the basic concept of writing a cross section in the form of Eq. (<ref>) which in the case of $e^+e^-$ collisions reads

\begin{align}
\begin{split}
d\tilde{\sigma}^{e^+e^-}(K_{\oplus},K_{\ominus},m^2)=\sum_{i,j=e^{\pm},\gamma}\int dx_{\oplus}dx_{\ominus}\Gamma^{e^+}_i(x_{\oplus},\mu^2,m^2)\Gamma_j^{e^-}(x_{\ominus},\mu^2,m^2)\\
\times d\sigma_{ij}(x_{\oplus}K_{\oplus},x_{\ominus}K_{\ominus},\mu^2)
\end{split}
\label{e+e-crosssection}
\end{align}

where the universal scale dependence of the PDFs is encoded in the solution of the DGLAP evolution equations

\begin{align}
\frac{\partial\Gamma_i(z,\mu^2)}{\partial\log \mu^2}=\frac{\alpha(\mu)}{2\pi}\left[P_{ij}\star \Gamma_j\right](z,\mu^2)\quad.
\label{DGLAPevolutioneq}
\end{align}

The PDFs here are with respect to the electron, i. e. $\Gamma_i=\Gamma_i^{e^-}$, where for convenience the superscripts are neglected.
If choosing a scale $\mu\sim Q$ we find the partonic cross sections $d\sigma_{ij}$ being completely describable by perturbative means. In QED, this is true for initial conditions of the PDFs as well, which is just a consequence of the nature of the abelian group theory.

The analytical form of the electron PDFs at LL precision has been first derived in Refs. <cit.> by using the following approach.

In general, the customary approach to solve Eq. (<ref>) is by transforming it from the partonic basis, here $i={e^{\pm},\gamma}$, into the one of singlet (S) and non-singlet (NS) combinations such that the corresponding PDFs are of the form

\begin{align}
&\Gamma_{\text{S}}=\Gamma_{e^-}+\Gamma_{e^+} &\Gamma_{\text{NS}}=\Gamma_{e^-}-\Gamma_{e^+}\quad.
\label{singlet-NSPDFs}
\end{align}

In this way we find the simplified equations

\begin{align}
\frac{\partial}{\partial\log \mu^2}\left(\begin{array}{c}
\Gamma_{\text{S}}\\
\Gamma_{\gamma}
\end{array}\right)&=\frac{\alpha(\mu)}{2\pi}\mathbf{P}_{\text{S}}\star \left(\begin{array}{c}
\Gamma_{\text{S}}\\
\Gamma_{\gamma}
\end{array}\right)\\
\frac{\partial\Gamma_{\text{NS}}}{\partial\log \mu^2}&=\frac{\alpha(\mu)}{2\pi}P_{\text{NS}}\star \Gamma_{\text{NS}}
\label{NSevolutionseq}
\end{align}

for which the non-singlet component decouples from the singlet-$\gamma$ system. The kernels in these equations read

\begin{align}
\mathbf{P}_{\text{S}}&=\left(\begin{array}{c c}
P_{SS} &P_{S\gamma}\\
P_{\gamma S} &P_{\gamma\gamma}
\end{array}\right),\qquad P_{SS}=P_{e^{\pm}e^{\pm}}+P_{e^{\pm}e^{\mp}},\quad P_{S\gamma}=2n_FP_{e^{\pm}\gamma},\quad P_{\gamma S}=P_{\gamma e^{\pm}}\\
\end{align}

with $n_F$ the number of lepton flavour families. $n_F=1$ is assumed in the following accounting for a single (massless) lepton flavour family in the calculation, e. g. the electron/positron in an $e^+e^-$ collision, which is a sufficiently accurate description according to Ref. <cit.>.

§.§ LO-LL approximation

Approximately, for the LL solution of the electron PDFs a one-dimensional flavour space can be assumed, i. e. $i=e^-$, such that the problem simplifies to solving Eq. (<ref>) for the non-singlet case, a scalar integro-differential equation, whereby

\begin{align}
\Gamma_{\text{NS}}(x,\mu^2)\simeq \Gamma_{e^-}(x,\mu^2)
\end{align}

As the initial condition of the electron PDF at LO is the Dirac delta function, as mentioned above, i. e.

\begin{align}
\Gamma_{e}(x,\mu_0^2)=\Gamma_{e}(x,m^2)=\delta(1-x),
\end{align}

the non-singlet evolution equation can be written in the form

\begin{align}
\Gamma_{e}(x,\mu^2)=\delta(1-x)+\int_{m^2}^{\mu^2}\frac{dq^2}{q^2}\frac{\alpha(q^2)}{2\pi}\int_x^1dzP_{ee}(z)\Gamma_{e}\left(\frac{x}{z},q^2\right)
\label{gameeEvol}
\end{align}

where for convenience we use `$e$' as subscript denoting the parton identity of the electron $e^-$. Using an auxiliary function $G(x,\mu^2)$ which is related to $\Gamma_e(x,\mu^2)$ as

\begin{align}
&G(x,\mu^2)=\int^1_xdt\Gamma_e(t,\mu^2) &\Gamma_e(x,\mu^2)=-\frac{\partial}{\partial x}G(x,\mu^2)
\label{gammaandG}
\end{align}

and the form of Eq. (<ref>) we find

\begin{align}
\label{Gevolution}
\end{align}

This equation can be solved for $G$ which provides the solution for the PDF evolution equations by the relation on the r. h. s. of Eq. (<ref>). By means of this, for the simplest case, using the one-loop accurate regularised (unpolarised) Altarelli-Parisi kernels

\begin{align}
& P_{ee}(z)=\langle \hat{P}_{ee}\rangle(z)-\delta(1-z)\int_{0}^{1}dt\langle \hat{P}_{ee}\rangle(t), &\langle \hat{P}_{ee}\rangle(z)=\frac{1+z^2}{1-z}\quad,
\label{Peeregularised}
\end{align}

the solution of Gribov-Lipatov resummation to all orders in $\alpha$ in the asymptotic region $x\simeq1$ <cit.>, i. e. for soft-collinear radiation of photons, can be obtained,

\begin{align}
\Gamma_e(x,\mu^2)=\frac{e^{\eta (\frac{3}{4}-\gamma_E)}}{\Gamma(1+\eta)}\eta (1-x)^{\eta-1}
\label{GribovLipatovSol}
\end{align}

with the definitions

\begin{align}
&\eta=\frac{\alpha}{\pi}L, &L=\log\frac{\mu^2}{m^2}
\label{etadefinition}
\end{align}

and $\gamma_E=0.577...$, the Euler-Mascheroni constant. Identically, this analytic form can be achieved by transforming the integro-differential evolution equations of Eq. (<ref>) into Mellin space by

\begin{align}
&M[f]\equiv f_N=\int_{0}^{1}dzz^{N-1}f(z) &M[g\star h]=M[g]M[h] \quad.
\end{align}

This simplifies the problem to solving ordinary differential equations, using the asymptotic limit $N\rightarrow\infty$ analogously to $z\rightarrow 1$ in $z$-space and an analytical Mellin inversion of the resulting solution according to Ref. <cit.>.

The recursive solution delivers an option to find $\Gamma_e(x,\mu^2)$ for arbitrary values $x<1$. Conceptionally, the procedure is as follows.

The function $G(x,\mu^2)$ as well as the PDFs $\Gamma_e(x,\mu^2)$ generally can be written as a perturbative series expressed as

\begin{align}
&G(x,\mu^2)=\sum_{n=0}^{\infty}G^{(n)}(x,\mu^2)=\sum_{n=0}^{\infty}\frac{\eta^n}{2^n n!}I_n(x),\label{Gseries}\\
&\Gamma_e(x,\mu^2)=\sum_{n=0}^{\infty}\Gamma_e^{(n)}(x,\mu^2)=\sum_{n=0}^{\infty}\frac{\eta^n}{2^n n!}\frac{\partial I_n(x)}{\partial x}\quad.
\label{GammaSeries}
\end{align}

Inserting $G(x,\mu^2)$ of this form in Eq. (<ref>) yields the recurrence relation

\begin{align}
\end{align}

which with Eq. (<ref>) and the boundary conditions $G^{(0)}(x,\mu^2)=G(x,m^2)=1$ implicating $I_0(x)=1$ and $I_1(x)=\int_{x}^{1}dz P(z)$ leads to

\begin{align}
I_n(x)=\int_{x}^{1}dz\langle\hat{P}_{ee}\rangle(z) \left[I_{n-1}\left(\frac{x}{z}\right)-I_{n-1}(x)\right] +I_{n-1}(x)I_1(x)\quad.
\end{align}

The divergence of $\langle\hat{P}_{ee}\rangle(z)$ for $z\rightarrow1$ in the first term is regularised in this way by the difference $I_{n-1}\left(x/z\right)-I_{n-1}(x)$. By this relation, $I_n$ in principle can be computed for each order $n$ iteratively which, however, gets increasingly complicated with each new iteration. This leads to the solution for $G(x,\mu^2)$ and by its derivative to that for $\Gamma^e_e(x,\mu^2)$. In order to match this iterative solution for arbitrary $x$ computed to a certain order in $\alpha$ with the asymptotic $x\simeq1$ Gribov-Lipatov solution to all orders in $\alpha$ described above the following has to be considered. Adding the solution of Eq. (<ref>) for $x\simeq1$ to the iterative one for all $x$ leads to a double counting of the soft terms in the $x\rightarrow 1$ regime up to the order of which the iterative solution is computed. Due to this, all terms of Eq. (<ref>) expanded again as a series in $\eta$ up to that order have to be subtracted again from the sum of the asymptotic $x\simeq1$ and the `all $x$' solution.
We thus find the standard LL accurate solution of electron PDFs with Gribov-Lipatov resummation to all orders in $\alpha$ matched with the `all $x$' solution up to $\mathcal{O}(\alpha^3)$ with the explicit form <cit.>

\begin{align}
\begin{split}
\Gamma_e(x,\mu^2)=&\frac{e^{\eta (\frac{3}{4}-\gamma_E)}}{\Gamma(1+\eta)}\eta (1-x)^{\eta-1}-\frac{1}{2}\eta(1+x)\\
&+\frac{1}{2^22!}\eta^2\left[-4(1+x)\log(1-x)+3(1+x)\log x-4\frac{\log x}{1-x}-5-x\right]\\
& + \frac{1}{1-x}\left[-\frac{3}{2}(1+8x+3x^2)\log(x)-6(x+5)(1-x)\log(1-x)\right.\\
&\left. -12(1+x^2)\log(x)\log(1-x)+\frac{1}{2}(1+7x^2)\log^2(x)\right.\\
\end{split}
\label{LLsoftPDF}
\end{align}

Note, that this is the electron PDF (or structure function) at LO-LL in $\alpha$ which matches the accuracy level of matrix elements at LO such that the overall cross section observables are LO-LL accurate.

§.§ NLO and NLO-NLL approximation

In order to increase the precision of observables to NLO in $\alpha$, terms up to (at least) the same order have to be taken into account coming from both the expansion of the PDFs and of the short distance partonic cross section. Schematically, by using perturbatively expanded forms of partonic cross sections and PDFs represented by terms of Eqs. (<ref>) and (<ref>), we symbolically rewrite Eq. (<ref>) for $\mu=Q$ as an expansion in $\alpha$,

\begin{align}
\begin{split}
\Gamma^{e^+}(Q^2)\star \Gamma^{e^-}(Q^2)\star d\sigma(Q^2)= \left(C^{(0)}_{\text{LL}}+\frac{\alpha}{2\pi} \log \frac{Q^2}{m^2} C^{(1)}_{\text{LL}}+\frac{\alpha}{2\pi} C_{\text{NLL}}^{(0)}+\mathcal{O}(\alpha^2)\right)\\
\star\left(d\sigma^{(0)}+\frac{\alpha}{2\pi} d\sigma^{(1)}+\mathcal{O}(\alpha^2)\right)\\
=\underbrace{C^{(0)}_{\text{LL}}\star d\sigma^{(0)}+\frac{\alpha}{2\pi}\left[\left(\log \frac{Q^2}{m^2} C^{(1)}_{\text{LL}}+C_{\text{NLL}}^{(0)}\right)\star d\sigma^{(0)}+C^{(0)}_{\text{LL}}\star d\sigma^{(1)}\right]}_{=d\sigma_{\text{NLO}}}+\mathcal{O}(\alpha^2).
\end{split}
\end{align}

We thus observe that the NLO fixed order cross section $d\sigma_{\text{NLO}}$, written symbolically in this form, implicitly contains the first $\mathcal{O}(\alpha)$ term of the expanded LL accurate PDFs by the term proportional to $\alpha/(2\pi)\log ({Q^2}/{m^2})$ as well as the initial conditions of NLL PDFs by $C_{\text{NLL}}^{(0)}$. By quantifying the factor

\begin{align}
\left[\frac{\alpha}{2\pi} \log \frac{Q^2}{m^2}\right]^2\sim\left[\frac{\alpha}{2\pi} \log \frac{s}{m^2}\right]^2
\label{estimatefNLOmassiveIS}
\end{align}

which appears in the leading logarithmic terms at each order in $\alpha$, which are not taken into account in the NLO cross section we gain a first estimate on the reliability of a fixed-order NLO massive initial state approximation. How the electron mass can affect fixed-order calculations with higher order perturbative corrections in $\alpha$ for $e^+e^-$ annihilation is studied in detail in Refs. <cit.>.

Since this fixed-order approximation however is highly dependent on the physics related to that observable the recommended method for a generic approach is to resum all logarithms at NLL accuracy into universal NLO-NLL accurate PDFs by solving the DGLAP evolution equations with kernels and initial conditions at NLO in $\alpha$. Convoluted with partonic cross sections at NLO for which the initial-state mass dependencies are removed this leads to perturbatively reliable NLO-NLL accurate predictions for observables at lepton colliders. Particularly, this level of precision is desired for high centre-of-mass energies of future collider setups as $Q^2\sim s \gg m^2$ implies an enhancement of the logarithmic terms.

As already mentioned above the initial conditions of NLL PDFs are of NLO accuracy in $\alpha$ and are different to the Dirac $\delta$ function as in the LO case since the leptons at NLO QED inherently mix with photons in the initial state. The solution for the NLO initial conditions of electron and photon PDFs are derived in Ref. <cit.> by means of an explicit short-distance cross section computation for a specific (but arbitrary) process, for which $e^+e^-\rightarrow u\bar{u}(\gamma)$ is chosen. This is achieved by computing the particle-level cross section, i. e. the l. h. s. of Eq. (<ref>), including terms up to $\mathcal{O}(m/\sqrt{s})$, on the one hand, and the parton-level cross sections with massless electrons, on the other hand. Solving Eq. (<ref>) with these analytical expressions for the electron PDFs yields the NLO initial conditions which – by applying the shorthand notation $\Gamma_i\equiv\Gamma^{e^-}_i$ – read

\begin{align}
\Gamma_i=\Gamma^{(0)}_i+\frac{\alpha}{2\pi}\Gamma^{(1)}_i+\mathcal{O}(\alpha^2)
\label{NLOIScond}
\end{align}

\begin{align}
\Gamma^{(0)}_i(x,\mu_0^2)&=\delta_{e^-i}\delta(1-x),\\
\Gamma^{(1)}_{e^-}(x,\mu_0^2)&=\left[\frac{1+x^2}{1-x}\left(\log\frac{\mu_0^2}{m^2}-2\log(1-x)-1\right)\right]_++K_{ee}(x),\\
\Gamma^{(1)}_{\gamma}(x,\mu_0^2)&=\frac{1+(1-x)^2}{x}\left(\log\frac{\mu^2_0}{m^2}-2\log x-1\right)+K_{\gamma e}(x),\label{gammaPDF}\\
\Gamma^{(1)}_{e^+}(x,\mu_0^2)&=0 \label{NLOinitconde+}
\end{align}

Here, the functions $K_{ee}(x)$ and $K_{\gamma e}(x)$ are renormalisation scheme dependent and equal to zero if the $\overline{\text{MS}}$ scheme is chosen. This is analogous to the definition of $K_{ad}$ in Eq. (<ref>)).
More details on their derivation and an alternative ansatz which is beyond the scope of this thesis are given in the reference stated above.

The only difference of the NLL compared to the LL evolution of the PDFs are the QED splitting kernels being of NLO accuracy, for which the two-loop corrections are given in Ref. <cit.>. Moreover,
the same methods as suggested for the LL evolution in Sec. <ref> can be applied: For the asymptotic $x\simeq 1$ solution the NS evolution operator obtained by solving the non-singlet DGLAP equations of Eq. (<ref>) in Mellin $N$-space with Altarelli-Parisi kernels at NLO in $\alpha$ is evaluated in the limit $N\rightarrow\infty$. This evolution operator Mellin inverted into $x$-space and convoluted with the NLO initial conditions then yields the non-singlet solution for the PDFs $\Gamma_{\text{NS}}$ in the limit $x\rightarrow1$.
Since the NLO kernel $\mathbf{P}_{\text{S}}=\mathbf{P}_{\text{S}}^{(0)}+\alpha/(2\pi)\mathbf{P}_{\text{S}}^{(1)}$ transformed into Mellin N-space according to Eq. (5.71) in Ref. <cit.> gets diagonal in the limit $N\rightarrow\infty$, the singlet PDF $\Gamma_{\text{S}}$ decouples from the photon PDF $\Gamma_{\gamma}$. Furthermore, the singlet-singlet component of the $\mathbf{P}_{\text{S}}$ kernel, i. e. $\lim_{N\to \infty}M[P_{SS}]$, yields the same expression as the corresponding one of the non-singlet case. Thus, the evolution operator for the decoupled singlet equations and – due to Eqs. (<ref>) to (<ref>) – the complete solution $\Gamma_{\text{S}}(x,\mu^2)$ are identical to those resulting from the non-singlet equations. Owing to these considerations as well as Eq. (<ref>) the asymptotic NLO-NLL accurate electron PDF describing soft radiation of photons with resummation to all orders in $\alpha$ is constructed as

\begin{align}
\Gamma^{\text{NLL}}_{e^-,\text{asy}}(x,\mu^2)=\left(\Gamma^{\text{NLO}}_{\text{NS}}(\mu_0^2)\star M^{-1}\left[\lim_{N\to \infty}E_{\text{NS},N}(\mu_0^2,\mu^2)\right]\right) (x,\mu^2)
\end{align}

with the NS evolution operator $E_{\text{NS},N}$ in Mellin $N$-space and

\begin{align}
\Gamma^{\text{NLO}}_{\text{NS}}(\mu_0^2)=\Gamma^{\text{NLO}}_{e^-}(\mu_0^2)=\Gamma^{(0)}_{e^-}(\mu_0^2)+\frac{\alpha}{2\pi}\Gamma^{(1)}_{e^-}(\mu_0^2)\quad.
\end{align}

The explicit form including a running of $\alpha$ is given by <cit.>

\begin{align}
\begin{split}
\Gamma^{\text{NLL}}_{e^-,\text{asy}}(x,\mu^2)=&\frac{e^{\hat{\xi}_1-\gamma_E\xi_1}}{\Gamma(1+\xi_1)}\xi_1 (1-x)^{\xi_1-1}\\
&\times \biggr\{1+\frac{\alpha(\mu_0)}{\pi}\left[\left(\log\frac{\mu^2_0}{m^2}-1\right)\left(A(\xi_1)+\frac{3}{4}\right)-2B(\xi_1)+\frac{7}{4}\right.\\
\end{split}
\label{NLLsoftPDF}
\end{align}

with quantities

\begin{align}
\xi_1&=2t-\frac{\alpha(\mu)}{4\pi^2b_0}\left(1-e^{-2\pi b_0t}\right)\left(\frac{20}{9}n_F+\frac{4\pi b_1}{b_0}\right),\\
\hat{\xi}_1&=\frac{3}{2}t-\frac{\alpha(\mu)}{4\pi^2b_0}\left(1-e^{-2\pi b_0t}\right)\left(\lambda_1+\frac{3\pi b_1}{b_0}\right),\\
t&=\frac{\alpha(\mu)}{2\pi}L-\frac{\alpha^2(\mu)}{4\pi}\left(b_0L^2-\frac{2b_1}{b_0}L\right)+\mathcal{O}(\alpha^2),\qquad L =\log\frac{\mu^2}{\mu^2_0} \label{tvariable}\\
b_0&=\frac{n_F}{3\pi},\quad b_1=\frac{n_F}{4\pi^2},\quad \lambda_1=\frac{3}{8}-\frac{\pi^2}{2}+6\zeta(3)-\frac{n_F}{18}(3+4\pi^2)
\end{align}

\begin{align}
A(\kappa)&=-\gamma_E-\psi_0(\kappa), \\
\end{align}

where $\psi_j(\kappa)$ is defined as:

\begin{align}
\psi_j(\kappa)&=\frac{d^{j+1}\log \Gamma(\kappa)}{d\kappa^{j+1}}\quad.
\end{align}

For small $\xi_1$, e. g. for a collider setup with energies at the order of hundred GeV where $\xi_1\sim 0.05$, the expansions

\begin{align}
&A(\xi_1)=\frac{1}{\xi_1}+\mathcal{O}(\xi_1), &B(\xi_1)=-\frac{\pi^2}{6}+2\zeta(3)\xi_1+\mathcal{O}(\xi_1^2)
\end{align}

show that the term proportional to $A(\xi_1)\log(1-x)$ is the dominant of Eq. (<ref>) even in extreme regions of $x$ in the limit $x\rightarrow1$. Due to $\xi_1\simeq \eta$ this is the essential difference of the asymptotic NLL solution compared to the LL counterpart, i. e. the first term of Eq. (<ref>).

For the solution of the electron PDFs at NLL accuracy valid for all $x$ values we adopt the iterative approach as suggested for the LL PDFs in Sec. <ref>. Comparing the resummed LL accurate terms of Eq. (<ref>) with those of NLL-accuracy in Eq. (<ref>), the manifest difference is the second term in the latter. Hence, the representation of the function $G$ and the PDFs $\Gamma_e$ in a perturbative series for NLL accuracy analogously to Eqs. (<ref>) and (<ref>) read

\begin{align}
&G^{\text{NLL}}(x,\mu^2)=\sum_{n=0}^{\infty}G^{\text{NLL},(n)}(x,\mu^2)=\sum_{n=0}^{\infty}\frac{\eta^n}{2^n n!}\left(I^{\text{LL}}_n(x)+\frac{\alpha}{2\pi}I^{\text{NLL}}_n(x)\right),\label{GseriesNLL}\\
&\Gamma^{\text{NLL}}_e(x,\mu^2)=\sum_{n=0}^{\infty}\Gamma_e^{{\text{NLL}},(n)}(x,\mu^2)=\sum_{n=0}^{\infty}\frac{\eta^n}{2^n n!}\left(\frac{\partial I^{\text{LL}}_n(x)}{\partial x}+\frac{\alpha}{2\pi}\frac{\partial I^{\text{NLL}}_n(x)}{\partial x}\right).
\label{GammaSeriesNLL}
\end{align}

This is used for the approach of solving the equation

\begin{align}
\label{GevolutionNLL}
\end{align}

with Altarelli-Parisi kernels at NLO in $\alpha$, $P(z)=P^{(0)}+\alpha/(2\pi)P^{(1)}$, an the relations

\begin{align}
&G^{\text{NLL}}(x,\mu^2)=\int^1_xdt\Gamma^{\text{NLL}}(t,\mu^2) &\Gamma^{\text{NLL}}(x,\mu^2)=-\frac{\partial}{\partial x}G^{\text{NLL}}(x,\mu^2)\quad.
\label{gammaandGNLL}
\end{align}

Considering a NLL evolution operator formally it is not justified to neglect the running of $\alpha$ anymore as in the case of LL PDFs. Moreover, the running of $\alpha$ from the Thompson to the EW scale phenomenologically yields a significant difference, i. e.

\begin{align}
\frac{\Delta \alpha}{\alpha(0)}=\frac{ \alpha(M_Z)-\alpha(0)}{\alpha(0)}\sim 0.07
\end{align}

It is thus convenient to use the variable <cit.>

\begin{align}
t=\frac{1}{2\pi b_0}\log \frac{\alpha(\mu)}{\alpha(\mu_0)}
\end{align}

which in an expanded form by means of

\begin{align}
\frac{\partial \alpha (\mu)}{\partial \log \mu^2}=\beta(\alpha)=b_0\alpha^2+b_1\alpha^3+\ldots
\end{align}

is already introduced in Eq. (<ref>), and we rewrite Eq. (<ref>) as

\begin{align}
G^{\text{NLL}}(x,t)=G^{\text{NLL}}(x,0)+\int_{0}^{t}du\frac{b_0\alpha^2(u)}{\beta(\alpha(u))}[P\star G^{\text{NLL}}]\left(x,u\right)\quad.
\end{align}

By this variable transformation, we find the perturbative series of Eq. (<ref>) modified to

\begin{align}
G^{\text{NLL}}(x,t)=\sum_{n=0}^{\infty}\frac{t^n}{ n!}\left(J^{\text{LL}}_n(x)+\frac{\alpha(t)}{2\pi}J^{\text{NLL}}_n(x)\right)\quad.
\end{align}

Eventually, we find the recursive relations

\begin{align}
J^{\text{LL}}_n &=P^{(0)}\star J^{\text{LL}}_{n-1},\\
\begin{split}
J^{\text{NLL}}_n&=(-1)^n(2\pi b_0)^nG^{\text{NLL},(1)}(\mu_0^2)+\sum_{p=0}^{n-1}(-1)^p(2\pi b_0)^p\biggr\{P^{(0)}\star J^{\text{NLL}}_{n-1-p}\\
&\qquad\qquad\qquad\qquad\qquad\qquad\quad+\left[P^{(1)}-\frac{2\pi b_1}{b_0}P^{(0)}\right]\star J^{\text{LL}}_{n-1-p}\biggr\}
\end{split}
\end{align}

These equations are solved iteratively for the non-singlet and singlet-$\gamma$ evolutions up to $J^{\text{LL}}_3$ and $J^{\text{NLL}}_2$ in Ref. <cit.> yielding NLL accurate non-singlet, singlet and photon PDFs $\Gamma^{\text{NLL}}_{\text{rec}}$ with resummation up to the respective order by means of Eq. (<ref>). Since the analytical expressions are lengthy and cover several pages they are not presented here explicitly.

Analogously to the LL accurate PDFs, a matching of the recursive `all $x$' solution, described by $\Gamma^{\text{NLL}}_{\text{rec}}$, to the asymptotic $x\simeq1$ one, $\Gamma^{\text{NLL}}_{\text{asy}}$ respectively, has to be applied. This can be achieved in an additive manner as proposed in <cit.> by

\begin{align}
\Gamma^{\text{NLL}}_{\text{mtc}}(x)= \Gamma^{\text{NLL}}_{\text{rec}}(x)+\left[\Gamma^{\text{NLL}}_{\text{asy}}(x)-\Gamma^{\text{NLL}}_{\text{subt}}(x)\right]F(x)
\label{PDFmatched}
\end{align}

with an arbitrary function $F(x)$ fulfilling the condition

\begin{align}
\lim_{x\rightarrow 1}F(x)=1
\end{align}

and terms $\Gamma^{\text{NLL}}_{\text{subt}}(x)$ which are present in both contributions, $\Gamma^{\text{NLL}}_{\text{rec}}(x)$ and $\Gamma^{\text{NLL}}_{\text{asy}}(x)$, and need to be subtracted exactly once to avoid the double counting.

Another approach in order to solve the DGLAP equations with kernels at NLO in $\alpha$ is by purely numerical means. For this, the DGLAP equations are solved in Mellin space for small scale intervals into which the scale evolution range is partitioned. Multiplying up all thus obtained evolution operators per scale interval yields the complete evolution operator in Mellin space. By Mellin inverting the product of the initial conditions in $N$-space and this evolution operator finally results in the solution for the PDFs in $x$ space. Even though this procedure is beneficial in a theoretical sense, such that NLL resummation can be achieved to all orders in $\alpha$ for the whole $x$ range, it involves the drawback on the statistics on which the numerical integration over the Mellin contour depends on.

A few remarks on the photon PDFs are left to be said: Formally for cross sections of $e^+e^-$ collisions – ignoring photons from beam spectra for the moment – photon PDFs would play a role only at the NNLO in $\alpha$: The initial conditions of Eq. (<ref>) start to contribute to $\Gamma_{\gamma}$ at $\mathcal{O}(\alpha)$ and partonic cross sections $d\sigma_{\gamma e^-}$,  $d\sigma_{e^+\gamma}$ and $d\sigma_{\gamma \gamma}$ on the r. h. s. of Eq. (<ref>) do not exist at the LO. Naively, one may conclude by this that the contributions from photon-induced channels are suppressed for an NLO calculation.
However, photon PDFs can have a sizeable impact if soft photon radiation effects are enhanced. The photon PDFs are dominant and largely enhanced against the electron PDFs for small $x$ values <cit.>. It is thus phenomenologically improper to neglect them for an NLO calculation.

The derivation of the solution for the NLL QED PDFs in this section is described in a rather schematic way. For more technical details and complete formulas the reader is encouraged to consider Refs. <cit.>.

§ MC INTEGRATION METHODS FOR ELECTRON PDFS

Challenges of the embedding of electron PDFs into numerical cross section computations in a Monte-Carlo generator framework encompass two aspects.
The numerical evaluation of NLL accurate electron PDFs per phase-space point requires an amount of computing time which is inappropriate for the purpose of a Monte-Carlo integration. This is facilitated by proper grid interpolation measures (similar to the case of proton PDFs) which will be covered by Sec. <ref>. Moreover, as leptonic parton distribution functions feature a divergent but integrable structure in the radiated energy an appropriate phase-space mapping is a necessary ingredient for their inclusion in a Monte-Carlo integration. How this can be done systematically for LO-LL and LO-NLL observables and as an approach for NLO-NLL cross sections in the FKS framework which will be covered in Sec. <ref>.

§.§ Electron PDF grid interpolation measures

The application of PDF grids and adequate interpolation measures for the electron PDFs are required for the same reason as it is commonly done for proton PDFs: The evaluation of PDFs depending on beam energy fraction $x$ and factorisation scale $\mu$ at higher orders in the coupling demands a tremendous amount of computing time due to the involved evolution procedure beyond LO. This is outlined in Sec. <ref> with respect to the NLL electron PDFs for which the evaluation of lengthy analytical expressions or the integration over a Mellin contour pose a time-consuming issue if used for a Monte Carlo phase-space integration. This especially concerns Monte-Carlo samplings with dynamic scales for which the evolution scale at which the PDFs are evaluated varies point by point.

A grid for which PDF values for fixed points of the $x$ range and the scale $\mu$ are saved as data sets with appropriate interpolation and optional extrapolation routines offers a possibility to circumvent this problem and is the well-known method of the  library <cit.>.
An adequate measure in order to interpolate points on a two-dimensional regular grid thereby is represented by the application of bicubic interpolation. Regarding the $x$-binning and interpolation for proton PDFs, the $x$-bin boundaries in small $x$ regions are arranged in constant distances of $\log_{10} x$ for the grid and mapped onto this logarithmic axis for the interpolation. This is done in order to account for large gradients of PDFs for $x\rightarrow0$ as it is the case for example for gluon PDFs in $pp$ collisions.
Conversely, this can be transferred to the binning and interpolation for the electron PDFs with singular $(1-x)^{\epsilon-1}$ behaviour close to $x=1$ as elaborated in Sec. <ref> which is present in both, the LL and NLL approximation, and enhanced in the latter at least by factors $\log(1-x)$ and $\log^2(1-x)$ according to Eq. (<ref>).

NLL electron PDF $f(x)$ (displayed as $\log_{10} f(x)$) as a function of $-\log_{10}{(1-x)}$ using exact values of the analytic function <cit.> depicted in red (`exact'), interpolated values with grid point mapping $(x,y)\rightarrow (\log_{10} x, y)$ in blue (`logx_interp.'), $(x,y)\rightarrow (\log_{10} x, \log_{10} y)$ in orange (`logy_logx_interp.') and $(x,y)\rightarrow (\log_{10} (1-x), \log_{10} y)$ in green (`logy_log1-x_interp.') with $y\equiv f(x)$ (upper plot); ratio $f(x)/f_{\text{exact}}(x)$ of the interpolated electron PDF function over the `exact' function depending on $-\log_{10}(1-x)$ (lower plot)

Using the code [The code is publicly available under .] <cit.> for numerical PDF values resulting from the evaluation of $\Gamma^{\text{NLL}}_{\text{mtc}}(x)$ in Eq. (<ref>) the appropriate binning and interpolation for the argument $x$ can be investigated.
Evaluating the electron PDFs for  $x$ sampled linearly within a range $0<x\le 0.9$ and logarithmically for $1\le -\log(1-x)\le 13$ and for a fixed scale $\mu$ yields the grid points between which the PDF values in the end are supposed to be obtained through interpolation. We can now probe several interpolation approaches by sampling $x$ between the fixed grid points and evaluate the PDFs for these in the two ways, once exactly and once by interpolation.
Bicubic interpolation routines according to those used for MSTW PDFs <cit.> are applied to interpolate PDF values between the grid points. The reliability of interpolated PDF values thereby highly depends on how the grid points are mapped which is apparent from the following observation. In Fig. <ref> comparisons of PDF values $f(x)\equiv \Gamma^{\text{NLL}}_{\text{mtc}}(x) $ obtained by exact evaluation (in red) and by interpolation with a specific mapping of the grid points $(x,y)\equiv(x,f(x))$ are shown. Absolute values are depicted in the upper plot and the ratio $f(x)/f_{\text{exact}}(x)$ for interpolated over exact values in the lower plot, respectively.
The interpolated PDF with mapping $(x,y)\rightarrow (\log_{10} x, y)$ which is used for example by  as explained above and the corresponding ratio shown as blue curves yield a deviation from the exact values of more than $1\%$. The interpolation is improved to $\lesssim 1\%$ deviation if a mapping $(x,y)\rightarrow (\log_{10} x, \log_{10} y)$ is applied (in orange). A mapping $(x,y)\rightarrow (\log_{10} (1-x), \log_{10} y)$ (green curves) renders agreement of PDF values for an exact and an interpolated evaluation well-below $0.01\%$ in the full range of $x$ fixed by the grid points with the exception of $x$ values extremely close to zero where the deviation is still only around $0.1\%$.
For MC predictions this accuracy can be considered as sufficiently high[Alternatively, one may interpolate using a grid point mapping $(x,y) {\rightarrow (\log_{10} (1-x), (1-x)y)}$. The accuracy of this method is not explicitly analysed here. However, due to the nearly linear shape of the function $f(x)(1-x)$ in Fig. <ref> a well performance of a bicubic interpolation is expected.].

Function $f(x)(1-x)$ representing the electron NLL PDF times the phase-space element $\Delta x$ as an approximation of the integrand in the vicinity of $x=1$.

Left to be discussed is the treatment of PDF values for $x$ values in the extreme region $x>x_{\text{cut}}$, where $x_{\text{cut}}=1-10^{-13}$ represents the upper boundary in $x$ for grid points of the analysis in this section. First of all, note that a fixed boundary $x_{\text{cut}}$, which naively would be assumed due to machine precision restrictions of $x$, is not required. This is because nothing prevents one from using a variable transformation, e. g. $\bar{x}=1-x$, and evaluate the inverted region in the limit $\bar{x}\to 0$.
This redefinition, however, must be substituted in the analytical expression for the electron PDF in order to obtain the grid points in this region. For this reason extrapolation with exemplary results for the different mappings of the grid points shown in Fig. <ref> is not essentially necessary since there is no restriction on the upper limit of $x$ for the grid points.
From the behaviour of the curves which largely deviate from the exact one for that specific region it furthermore becomes apparent that if proceeded in the traditional way extrapolation is rather useless for the singular enhanced PDFs.

Since the singularity of the electron PDF is integrable we can consider the following fact from the MC integration perspective.
Approximately, for high-energy collider experiments matrix elements can be considered numerically constant varying $x$ in the region extremely close to one. Broadly speaking, the contributions to the integral in this region for one beam axis thus come from

\begin{align}
&\lim_{x\rightarrow1}dI(x)=dxf(x)c, &c=\text{const.},
\end{align}

i. e. the PDF times the phase-space element $dx$ and a constant number $c$. We can now investigate at which $x$ these contributions eventually become negligible: Varying a technical cut in the range $x<x_{\text{cut}}<1$, above which the PDFs are set constant to $f(x_{\text{cut}})$, would not change the integrated result within a statistical error which is phenomenologically tolerable. This can be understood by estimating the contribution of the right-most bin of the $x$ range in a MC integration as negligible. Functionally, it can be approached by

\begin{align}
\Delta x f(x)c\simeq(1-x)f(x)c\quad.
\end{align}

The function $(1-x) f(x)$ depending on $-\log_{10}(1-x)$ is plotted in Fig. <ref> with respect to an exact evaluation of the electron PDFs.
It can be seen that the area of the right-most bin is shrinking the closer $x$ approaches one which by following the trend of the `exact' curve below $1-x=10^{-16}$ yields the prerequisite for imposing a technical cut in general. In particular, it can be observed that the area $(1-x)f(x)$ for $1-x\simeq10^{-16}$ would be about $10\%$ of the contribution to the integral with respect to that for $1-x\simeq10^{-4}$ which is too large for a technical cut set to $1-x_{\text{cut}}=10^{-16}$. The investigation of the range in the vicinity of $x=1$ for which a technical cut would not affect the integrated MC cross section result within sufficiently high statistics will be deferred to future works.

For the interpolation of the PDFs between grid points with distinct scale values $\mu$ and fixed $x$ we note that the scale dependence of electron PDFs is merely due to logarithms $\log (\mu^2/\mu_0^2)$, e. g. through variables $t$ or $\xi_1$ in Eq. (<ref>). Since $\mu\geq\mu_0$, a polynomial approximation of the functional dependence of the PDFs on the scale is thus accurate enough such that the standard bicubic interpolation between linearly distributed grid points along the range $\mu_0^2\leq \mu^2 \leq \mu^2_{\text{max}}$ can be considered reliable.

§.§ Mapping of FKS phase-space for divergent electron PDFs

The peaked structure of electron PDFs at both, LL and NLL accuracy, imposes the need for a revised MC sampling of phase-space points in order to account for the full area of extremely small sections of the integration variable $x$ for the beam energy fraction.
The standard mapping applied in  in the context of ISR structure functions for an improved integration efficiency for LO calculations proceeds as follows.

We define mapping variables $\{p_1,p_2\}$ with $p_1,p_2\in[0,1]$ depending on the random number variables $\{r_1,r_2\}$ with $r_1,r_2\in[0,1]$ and using $\bar{r}_1=1-r_1$ as

\begin{align}
p_1 &= 1 - (1-r_1)^{1/\epsilon}=1- \bar{r}_1^{1/\epsilon} \label{p1 random}\\
p_2 &=
\begin{cases}
1 - {(2r_2)^{1/\epsilon}}/{2}, & u>0 \\
{(2r_2)^{1/\epsilon}}/{2}, & u<0\\
1/2, & u=0
\end{cases}
& u = 2r_2-1
\end{align}

where $\epsilon$ is an arbitrary small number. In  this number is set by the customary choice

\begin{align}
\epsilon=\frac{\alpha}{\pi}q_e^2\left(\ln \frac{s}{m^2}-1\right)\quad,
\end{align}

with $q_e$ the charge and $m$ the mass of the initial-state particle.
Integrating over $p_1$ and $p_2$ results in the Jacobian

\begin{align}
f_1 &= f_{p_1}f_{p_2}
\end{align}

\begin{align}
f_{p_1} &= \frac{1}{\epsilon}(1-r_1)^{1/\epsilon-1}=\frac{1}{\epsilon}\bar{r}_1^{1/\epsilon-1}\\
f_{p_2} &=
\begin{cases}
\frac{(2-2r_2)^{1/\epsilon-1}}{\epsilon}, & u>0 \\
\frac{(2r_2)^{1/\epsilon-1}}{\epsilon}, & u<0\\
1/\epsilon, & u=0
\end{cases}\quad.
\end{align}

This mapping gives an enhanced sampling of points at the endpoints of the mapping variables, i. e. for $p_1$ in the vicinity of $1$ and for $p_2$ approaching both, $0$ and $1$.

The dependence of $\{p_1,p_2\}$ on the physical variables $\{x_1,x_2\}$ with $x_1,x_2\in[0,1]$, the fractions of the initial state momenta, is further defined as

\begin{align}
&p_1 = x_1 \cdot x_2  &p_2 = \frac{\log x_1}{\log (x_1\cdot x_2)}
\end{align}

such that the physical variables $\{x_1, x_2\}$ to be integrated over, in turn, depend on the mapping variables as

\begin{align}
&x_1 = p_1^{p2} &x_2 = p_1^{1-p2}.
\label{p variables}
\end{align}

On the one hand, this mapping parametrisation fulfills the criterion of mapping the unit square such that the parton centre-of-mass energy is constant with respect to the centre-of-mass frame, i. e. $x_1 \cdot x_2 $ is independent from $p_2$. On the other hand, integrating over $\{x_1,x_2\}$ results in the Jacobian

\begin{equation}
f_2 = |\log p_1 |
\end{equation}

which flattens the integrand while sampling $p_1$ close to one which follows from Eq. (<ref>), and thus sampling $f_2$ close to zero. The overall Jacobian for Born or LO computations then reads

\begin{align}
f_{\text{Born}}=f_1f_2 = f_{p_1}f_{p_2}|\log p_1|.
\label{jacobianborn}
\end{align}

In this way, the integration and adaption in MC simulations is well-behaved and yield useful predictions as it is shown by simulations with  <cit.>. Crucially, the Jacobian in Eq. (<ref>) depends on the factor `$(1-r_1)^{1/\epsilon-1}|\log p_1|$' which flattens the integrand for $r_1\to 1$. Consequently, this factor can be used paradigmatically for similar parametrisations for NLO computations proposed in the following.

The condition that the total invariant mass is conserved if a parton is radiated from the system induces a rescaling of $x_1$ and $x_2$ for initial-state emissions according to Eq. (<ref>). If applying the standard parametrisation in random number variables, this leads to a breakdown of the mapping benefits, which are described above, for non-Born initial-state momenta. Essentially, for the real-emission and DGLAP-remnant component the rescaled beam momenta (before photon emission off the initial-state) are shifted closer to $1$ with respect to the unrescaled momenta. By denoting the rescaled momenta as $x^{\prime}=x+\delta x$ one can find that even for a small variation $\delta x > 0$ which is in line with the condition $x^{\prime}\le 1$ the ratio of the rescaled over the unrescaled PDFs diverges as $x^{\prime}$ approaches $1$, i. e.

\begin{align}
\lim_{x^{\prime}\rightarrow 1} \frac{\Gamma(x^{\prime})}{\Gamma(x)}=\lim_{x\rightarrow 1-\delta x}\frac{\Gamma(x+\delta x)}{\Gamma(x)}\rightarrow \infty\quad.
\end{align}

This means that in addition to the parametrisation of the unrescaled Born momenta $x_1$ and $x_2$ according to Eq. (<ref>) an improved mapping parametrisation of the variation $\delta x$ is essential.
Thus, in order to recover the previously explained integration efficiency due to the mapping we have to revise the parametrisation explicitly for rescaled initial state momenta $x_1^{\prime}$ and $x_2^{\prime}$ and thereby retaining the dependence on the FKS variables.

For simplicity, we start here with the formulation of the additional NLO integration parameters for the remnant of the subtraction of collinear ISR singularities in integrated form, i. e. the degenerate $n+1$ contribution (or DGLAP remnant) of Eq. (<ref>). For this contribution one has to integrate over exactly one additional degree of freedom compared to those of the Born cross section. The momentum dependence of the PDFs rescales as

\begin{align}
\Gamma(x_j,\mu) \longrightarrow \Gamma(x_j/z_j,\mu)
\end{align}

with $x_j < x_j/z_j < 1$ and emitter $j\in\{1,2\}$.
The standard parametrisation of ISR remnant integration variables $z_j$ following the restriction of Eq. (<ref>) is defined in terms of random number variable $r_z\in [0,1]$ as

\begin{equation}
z_j = 1 - r_z (1-x_j)
\end{equation}

with Jacobian

\begin{equation}
f_z= 1-x_j.
\end{equation}

Now, we rewrite these formulas for an improved mapping adequate for using NLL electron PDFs in NLO computations.
In order to enhance the Monte-Carlo sampling of points for $z_j$ close to $1$, it is useful to apply a bijective mapping of the random variable $r_z \in [0,1]$ defining, similar to Eq. (<ref>),

\begin{align}
p_z = 1 - (1-r_z)^{1/\epsilon}
\end{align}

with Jacobian factor

\begin{align}
f_{p_z} &= \frac{1}{\epsilon}(1-r_z)^{1/\epsilon-1}.
\end{align}

Regarding the condition $[0,1] \longmapsto [x_j,1]$ mapping $p_z \longrightarrow z_j$ we can find the parametrisation

\begin{align}
z_j = 1- p_z (1-\log p_z) (1-x_j) %+x_j
\label{improvedz}
\end{align}

which results in the factor

\begin{align}
f_{z_j}=(1-x_j)\log p_z
\end{align}

for the overall Jacobian per emitter

\begin{align}
f_{\text{DGLAP,j}}= f_{p_z}f_{z_j}=\frac{1}{\epsilon}(1-r_z)^{1/\epsilon-1}(1-x_j)\log p_z.
\label{DglapJacobian}
\end{align}

In this way, the integrand of the DGLAP remnant gets flattened which is useful with respect to the enhancement of the PDF close to $x_j=1$ for the corresponding emitter $j$.

Analogously, we proceed for the real emission phase-space parametrisation and for the IS collinear subtraction terms of the real component which both require rescaled IS momenta entering the PDFs. For the first, it is less trivial to find an adequate mapping for the rescaled momenta. It requires to get Jacobian factors analogous to Eq. (<ref>) resulting in the full integration efficiency effect as stated above and, simultaneously, retaining the actual FKS parametrisation in both, $\xi$ and $y$, defined in Eq. (<ref>). According to these definitions, for the real component the momentum dependencies of the PDFs rescale as

\begin{align}
&\Gamma(x_1,\mu) \longrightarrow\Gamma(x_1^{\prime})= \Gamma\left(\frac{x_1}{\sqrt{1-\xi}}\sqrt{\frac{2-\xi(1-y)}{2-\xi(1+y)}}\right)\\
&\Gamma(x_2,\mu) \longrightarrow\Gamma(x_2^{\prime})= \Gamma\left(\frac{x_2}{\sqrt{1-\xi}}\sqrt{\frac{2-\xi(1+y)}{2-\xi(1-y)}}\right)
\end{align}

which in the collinear limits $y=\pm 1$ have to coincide with

\begin{align}
\end{align}

entering the collinear subtraction terms. Because of the diverging gradient of the PDFs at $x_1\rightarrow1$ and $x_2\rightarrow1$ a proper local cancellation in the collinear limit formally guaranteed by

\begin{align}
&\lim_{y\rightarrow 1}\Gamma(x_1^{\prime})\Gamma(x_2^{\prime})\hat{\mathcal{R}}(\xi,y)-\Gamma\left(\frac{x_1}{1-\xi}\right)\Gamma\left({x_2}\right)\hat{\mathcal{R}}(\xi,1)\sim 0\\
&\lim_{y\rightarrow -1}\Gamma(x_1^{\prime})\Gamma(x_2^{\prime})\hat{\mathcal{R}}(\xi,y)-\Gamma\left({x_1}\right)\Gamma\left(\frac{x_2}{1-\xi}\right)\hat{\mathcal{R}}(\xi,-1)\sim 0
\end{align}

gets numerically affected for a traditional MC sampling of $\xi$ and $y$ which is used for hadron collisions.
Due to the variable $\xi$ appearing in both, the real emission amplitudes and the collinear subtraction terms, one is tempted to merely reparametrise $\xi$ analogously to the variable $z$ for the DGLAP remnant. However, one has to keep in mind that the upper boundary $\xi_{\text{max}}$, and in this turn $\xi$ itself, as well as the rescaled variables $x_1^{\prime}$ and $x_2^{\prime}$ are correlated with $y$. In this sense, it is mandatory to find a two-dimensional phase-space mapping for an improved efficiency.

An idea to technically realise such a real phase-space mapping is based on the following consideration. Note that we have exactly the same number of degrees of freedom for integration variables defining the initial-state real momenta, i. e. $\{x_1,x_2,\xi,y\}$, as the number of variables which needs to be adapted for an efficient phase-space evaluation, i. e. $\{x_1,x_2,x_1^{\prime},x_2^{\prime}\}$. In the same way as $x_1^{\prime}$ and $x_2^{\prime}$ is defined by Eq. (<ref>) through randomly determined $x_1$, $x_2$, $\xi$ and $y$ values, we can conversely derive $\xi$ and $y$ from the variables $x_1$, $x_2$, $x_1^{\prime}$ and $x_2^{\prime}$ defined through random numbers with the condition

\begin{align}
&x_1\leq x_1^{\prime}<1 &x_2\leq x_2^{\prime}<1\quad.
\end{align}

Hence, we define the rescaled variables with mapping analogously to Eq. (<ref>) as

\begin{align}
&x^{\prime}_j=1- \hat{p}_j (1-\log \hat{p}_j) (1-x_j) &j=1,2
\label{improvedmappingreal}
\end{align}

with $\hat{p}_j\in [0,1]$ which is constructed from random numbers $\hat{r}_j\in [0,1]$ as

\begin{align}
\hat{p}_j = 1 - (1-\hat{r}_j)^{1/\epsilon}
\end{align}

and leads to Jacobians

\begin{align}
f_{x^{\prime},j}=\frac{1}{\epsilon}(1-\hat{r}_j)^{1/\epsilon-1}(1-x_j)\log \hat{p}_j
\end{align}

In order to find the expressions of $\xi$ and $y$ in terms of $x_1^{\prime}$, $x_2^{\prime}$, $x_1$ and $x_2$ regarding the conditions of massless ISR construction according to Sec. <ref> we define auxiliary quantities

\begin{align}
&A\equiv\frac{x_1x_2^{\prime}}{x_2x_1^{\prime}}=\frac{2-\xi(1+y)}{2-\xi(1-y)} &B\equiv\frac{x_1x_2}{x_1^{\prime}x_2^{\prime}}=1-\xi
\end{align}

such that $\xi$ and $y$ can be derived yielding

\begin{align}
&\xi = 1-B &y=\left(\frac{1+B}{1-B}\right)\left(\frac{1-A}{1+A}\right)\quad.
\end{align}

\begin{align}
d\xi dy=\mathcal{J}_{1}(A,B)~dAdB=\mathcal{J}_{1}(A,B)\mathcal{J}_{2}(x_1^{\prime},x_2^{\prime})~dx_1^{\prime}dx_1^{\prime}
\end{align}

\begin{align}
&\mathcal{J}_{1}(A,B)=2\left(\frac{1+B}{1-B}\right)\frac{1}{(1+A)^2} &\mathcal{J}_{2}(x_1^{\prime},x_2^{\prime})=2\frac{x_1^2}{x_1^{\prime 3}x_2^{\prime}}
\end{align}

we get the final Jacobian factor for $\xi$ and $y$ parametrised in random numbers $\hat{r}_{1/2}$,

\begin{align}
\label{jacobianrealmapping}
\end{align}

This parametrisation of $\xi$ and $y$ for the radiated momentum in terms of random number variables is supposed to yield the desired improved integration efficiency due to `$(1-\hat{r}_j)^{1/\epsilon-1}\log \hat{p}_j$' factors in the Jabobian. These factors flatten the integrand if $\hat{r}_j\to 1$ which is analogous to similar factors in the Jacobians of the Born and DGLAP remnant, i. e. in Eq. (<ref>) and (<ref>). The numerical proof is deferred to future works. However, an overview on implemented elements and first successful consistency checks achieved on the way to the complete proof is given in Sec. <ref>. The desired effect for an efficient integration is per definition constrained to ISR. For FSR the radiated momentum is constructed from the same FKS variables $\xi$ and $y$ as used for the ISR singular regions. However, as the beam momentum fractions entering the PDFs receive no rescaling for FSR no additional efficiency improvement for the integration of the FKS variables is needed. In this sense, the parametrisation of $\xi$ and $y$ as shown above is redundant for FSR, but eventually has no effect on the numerical result of the MC integration.

CHAPTER: AUTOMATED NLO EW CALCULATIONS IN THE MC FRAMEWORK WHIZARD

For simulating physics with a predictive power at the level of exclusive data Monte-Carlo event generators in general represent indispensable tools. The multi-purpose program  offers the possibility of producing cross sections and simulated event samples for generic processes. Its general capabilities are introduced in Sec. <ref>. This framework is used for an automation of precision calculations at NLO in SM couplings for which the extension to automated NLO EW and NLO mixed corrections is the particular aim of this thesis. The technical details of the NLO automation in  are presented in Sec. <ref>. As formally described in Sec. <ref>, NLO contributions to arbitrary coupling powers might require the subtraction of both, QCD and EW singularities. How both correction types are treated for one NLO contribution at fixed order is explained in Sec. <ref>. Moreover, in this section challenges for IR-safe observables imposed by EW and, in particular, mixed corrections will be discussed as well.

§ THE EVENT GENERATOR WHIZARD

[The program is publicly available under .], an acronym for `$W$, $H$iggs, $Z$ And Respective Decays', is a multi-purpose Monte-Carlo event generator providing the complete system of particle-physics calculations, from setting up a Lagrangian of a chosen model up to simulated unweighted event samples.

A wide range of physics models, the SM as well as simplified and extended versions of it, are implemented as internal models which is in addition to the interface to  <cit.> and  <cit.> providing arbitrary models in the Universal Feynman Output (UFO) format <cit.>. Scattering amplitudes for arbitrary tree-level processes of a specified model are constructed and evaluated automatically by the intrinsic matrix element generator  <cit.> based on the colour-flow formalism <cit.>.
For NLO computations virtual loop and colour-/spin-correlated matrix elements in general are accessed by one-loop providers (OLPs)  <cit.>,  <cit.> and  <cit.>. The integration and event generation in  happens via  (Vegas AMPlified) <cit.> using  grid and weight adaptive integration methods based on the  algorithm <cit.> combined with multi-channel principles <cit.>. The generation of multiple phase-space channels with the heuristic algorithm implemented in  increases the efficiency of an adaptive integration even beyond the factorisability of effective integrands. For the phase-space evaluation on multi-core architectures the  algorithm has been superseded by  <cit.> by which the integration is parallelised via the message-passing interface (MPI). This capability is essential especially in view of high multiplicity processes requiring a large number of phase-space dimensions to be integrated over or CPU-intensive NLO computations. Arbitrary phase-space cuts can be set in the  steering () files via commands in the scripting language .

As event samples are of usage for comparison with experimental data only if they are unweighted the phase-space generator involves a rejection algorithm for the conversion of weighted into unweighted samples. The showering and hadronisation of events which converts partonic into hadronic events proceeds via internal algorithms <cit.> or interfaces to external tools as  <cit.>, by default shipped with , or  <cit.>. For jet clustering algorithms in order to render observables of processes involving jets IR-safe the package  <cit.> is interfaced. The output formats of event files which can be written out directly by  are , , ,  and  <cit.>. Further common formats can be used by linking  to external packages such as  <cit.>,  <cit.> or  <cit.> with an interface of the latter for exporting events in -trees <cit.>.

Concerning NLO computations, the subtraction scheme for the regularisation of IR singularities chosen for  is the FKS scheme due to its advantages mentioned already in Sec. <ref>. The formalism to achieve cross sections and observables at fixed NLO in couplings $\alpha$ or $\alpha_s$ is explained in detail in the previous chapter. For the generation of hadronic events at the precision level of NLO perturbative expansions the following caveat has to be considered.
In order to circumvent a double counting of terms describing the radiation of partons in both, the matrix elements and the parton shower, a proper matching scheme has to be applied. The method used in  is the POWHEG matching scheme <cit.> which perfectly interplays with the FKS subtraction scheme. This is
because the algorithm is designed for the generation of NLO distributions relying on unique Born phase-space points due to its description with Sudakov form factors based on ratios of real emission over Born squared matrix elements. Furthermore, the POWHEG method – related to its name which is an acronym for `Positive Weight Hardest Emission Generator' – beats other matching schemes with respect to avoiding negative event weights which is an inherent problem in MC@NLO <cit.> for example. The POWHEG matching scheme in  first has been applied in the context of $t\bar{t}(H)$ production at an $e^+e^-$ collider <cit.> and has recently been extended to $pp$ collisions <cit.>. In order to increase the precision of event samples, merging schemes are applied which combine matrix elements with different final-state parton multiplicities and parton showers.  uses the MLM scheme (of M. L. Mangano) <cit.> which is based on slicing the radiation phase-space into domains for the matrix elements and the parton shower.

The framework of  has been designed with a focus on physics simulation at lepton collisions. As one of the most established generators in this field it induced many studies and design reports for future colliders as ILC, CLIC, CEPC or FCC-ee <cit.>.
This is owing to its broad range of beam structure features. It covers effects as beamstrahlung, a strong electromagnetic field induced by dense beam bunches, which is taken into account by means of the packages  <cit.> and  by fitting the energy spectra obtained from  <cit.>. Additionally, there is the option for simulations to use Gaussian beam energy spreads.
ISR structure functions for lepton initial states are implemented to the LL accuracy with the formulas of Refs. <cit.> as presented in Sec. <ref> and (recently) to the NLL accuracy with the analytical form of Refs. <cit.>. Photon initial-states, e. g. inducing background processes at lepton colliders, can be simulated with  by the so called equivalent photon approximation (EPA) based on the concept of Weizsäcker-Williams approximation <cit.>.
Another feature with respect to the lepton beams is the possibility to include polarisation of beams for arbitrary setups, e. g. choosing polarisation beam fractions, axes or types, i. e. circular, longitudinal or transversal polarisation. Further beam structure features for lepton collisions provided by  are asymmetric beams and crossing angles, for example.

For the evaluation of PDFs in hadronic collisions there is an interface to the PDF set library  <cit.> and additionally built-in grids for specific PDF sets.

From the technical side,  is a modular written  code using the modern  standard.

§ METHODOLOGY OF AUTOMATED NLO CORRECTIONS

The first NLO calculations achieved with  were with respect to chargino production at the ILC using special-tailored NLO-EW amplitudes <cit.> and the process $pp\to b\bar{b}b\bar{b}$ at NLO QCD <cit.> based on CS subtraction. The work on fully-automated NLO computations using the FKS scheme in  started by considering pure QCD corrections with a first proof of concept for off-shell top quark production at a lepton collider given by Refs. <cit.>. Further applications of the NLO setup related to top quark physics followed by Refs. <cit.>. The capability of  computing NLO QCD corrections to lepton and hadron collision processes with arbitrary external states in a fully-automated way has been demonstrated by Refs. <cit.>. The extension of this framework accounting for NLO EW corrections for generic processes is the aim of this work.

Basic ingredients of an NLO computation, the storage of matrix elements and the one-loop provision by external tools, are presented below. Furthermore, a brief overview on the workflow of an NLO calculation is given in Sec. <ref> and on basic EW code supplements in Sec. <ref>.

§.§ Bookkeeping of matrix elements

The bookkeeping of matrix elements is an important ingredient of an automated framework computing observables for arbitrary processes, in particular for those with external states of a multi-dimensional flavour space represented by protons or jets. In , the bookkeeping is realised by state matrices implemented as tries, i. e. ordered tree data structures, with nodes representing sets of quantum numbers each. In this way a subprocess is fully described by the chain of nodes from the root to one of the last nodes. A set $i$ consists of the following elements: the particle identity $f_i$, also called `flavour', the colour flow indices $c_i$ which are assigned according to the colour flow formalism <cit.>, helicity indices $h_i$ and the momentum $\{k\}_i$.
The index $i$ thereby represents the position of a particle in Born and real flavour structures  and  with particles ordered in a code-optimising way, e. g. with the last index of  referring to the radiated parton. In accordance with the FKS scheme their entries are accessed and linked in the code based on the meaning of positions of the process particle structures $f_b$ and $f_r$ defined in Eq. (<ref>).
The leaves of the trie data structure finally correspond to the stored values of the (squared) matrix elements which are assigned to each possible branch of the trie, hence to each subprocess. The way of storing the data in tries in general relies on few memory space and allows for optimised searching algorithms.

As already mentioned above, tree-level matrix elements can be obtained by the internal matrix element generator . For NLO computations it is however possible to rely fully on external tools for tree-level, colour-/spin-correlated and virtual loop matrix elements. In some cases this is even required considering the status of the implementation at the time this thesis is written, e. g. with respect to fixed coupling powers which will be commented in more detail in Sec. <ref>.

§.§ One-loop providers

For the construction and evaluation of one-loop amplitudes the OLPs which can be linked to  rely on the following methods.

The tool  uses algebraic methods generating analytical expressions for Feynman diagrams which in a next step are simplified by $D$-dimensional reduction at integrand level <cit.> using  <cit.> or  <cit.>. Alternatively, tensor reduction <cit.> can be applied using  <cit.>. The evaluation of scalar one-loop integrals can be performed with libraries such as  <cit.>.

Contrary to this method,  relies on a numerical algorithm for constructing one-loop amplitudes in a recursive approach. The numerical recursion based on a similar algorithm as the one presented in Ref. <cit.> is thereby applied to factorised colour-stripped cut-open loop diagrams. The amplitudes are evaluated using tensorial reduction by an `on-the-fly' method <cit.>, supplied with scalar one-loop integrals from  <cit.> or . Alternatively, for the tensor reduction  or  <cit.> using OPP reduction can be used. Due to its numerical approach  represents a fast-evaluating and efficient OLP for MC integration purposes.

uses a complete recursive approach based on tree-level off-shell currents as basic building blocks. It therefore avoids the bottleneck of recomputing identical subgraphs appearing in several different Feynman diagrams, as other OLPs do, which can get severely inefficient especially for high multiplicity processes. The evaluation of one-loop scalar and tensor integrals proceeds via the  library.

While the interface to  is implemented as an own module, for  and  the communication with  relies on the Binoth Les Houches Accord (BLHA) <cit.>, a standard interface between one-loop providers and MC generators. This BLHA interface in the  code is complete in the sense that the interplay with the system of storing the matrix elements values in terms of state matrices is guaranteed. Taking this and the facts described above into account, at the moment  is considered as the preferred tool for NLO calculations to processes consisting of multiple subprocesses such as $pp$ induced processes. While all OLPs presented here are in general capable of computing NLO EW amplitudes,  is the most advanced tool considering the numerically efficient computation of EW one-loop amplitudes for arbitrary processes, including also mass dependencies of all fermions at higher orders in EW.

§.§ Workflow of NLO calculations

An NLO calculation in  in general proceeds in the following way, sketched illustratively in Fig. <ref>.

After all process flavour structures at the Born tree-level  are generated for a user-defined process the module  automatically constructs all real flavour structures  based on possible SM vertices. This proceeds by inserting partons firstly in an arbitrary way in place of the proton and jets[The process specified by the user possibly contains particle container objects, e. g. for multiple partons in protons and jets, defined as  in the  steering file which in general increases the number of subprocesses.], `' and `', if included in the process definition, and as an additional radiated parton attached to the process. The allowed partons to be inserted thereby are represented by elements of the set $\{\texttt{pr},\texttt{j},\texttt{x}\}$, for which `' and `' are otherwise excluded if not used in the process definition. Depending on the correction type of the NLO calculation, `' is defined by either the gluon `' in the case of QCD corrections or by the photon `' in the case of EW corrections. All combinations of particle sequences thus obtained are then sorted such that only those corresponding to valid SM vertices are kept.

In the next step the FKS singular regions are initialised by applying the routine . Each of the real flavour structures  is first analysed for FKS pairs possibly inducing a collinear and/or soft singularity. The singular regions are constructed for each FKS pair of  which are found as follows: For all pairs of legs which represent partons the emitter parton is traced back by either QCD or QED splittings such that a `test' underlying Born structure is constructed.
The pair of legs is then recognised as FKS pair if the test Born structure is contained in the list of pre-constructed Born flavour structures .
The index of each thus obtained singular region is linked to the position of the emitter particle  in , the  representing the positions of the partons of the FKS pair and the real flavour structure  itself.
Subsequently, the preliminary singular regions are sorted according to all prescriptions necessary to define an FKS region as described in Sec. <ref>.
Eventually, FKS regions are initialised with index  to which additionally all  with distinct emitter but same  are assigned in order to compute the complete set of factors $d_{ij}$ for $S_{(i,j)}$ according to Eq. (<ref>).
Furthermore, for the later construction of subtraction terms each  is associated with an underlying Born flavour structure  and a dynamic correction type , which will be discussed in more detail in Sec. <ref>. A complete set of FKS regions assembled in the so-called `FKS table' for an exemplary process is presented in App. <ref>.

Describing the procedure of NLO calculations for processes with multiple subprocesses it will be referred to the tool  in the following, for which the interface is the most extensively tested one in . Arrays consisting of all  and  which are part of the initialised FKS regions are the input of the routine . This routine controls all input given to OLPs using the BLHA interface by writing out  files with prefixes containing the name for a separate component of the NLO computation in .
By this, the squared matrix elements of the Born contribution are requested for in the  file. For the real part the file   contains requests for the real squared matrix elements, and for the subtraction terms  involves those corresponding to Born, colour- and spin-correlated squared matrix elements, respectively.  requests for all loop-tree interfering, colour- and spin-correlated as well as Born squared matrix elements for the virtual finite and integrated subtraction terms. In case PDFs are used,  asks for Born squared matrix elements entering the DGLAP remnant component.
An exemplary  file is given in App. <ref>. After these files are written out by ,  is called, registers all subprocess requests and responds with  answer files which contains a numbering of all squared matrix elements. This numbering is read in by  and in a further step linked to  objects (the branches of the trie) of the state matrix described in Sec. <ref>.

[node distance=2.cm]
startstop = [rectangle, rounded corners, minimum width=3cm, minimum height=1cm,text centered, draw=black, fill=red!30]
process = [rectangle, rounded corners, minimum width=3cm, minimum height=1cm, text centered, draw=black, fill=orange!30]
OLP = [rectangle, rounded corners, minimum width=3cm, minimum height=1cm, text centered, draw=black, fill=blue!30]
kinematics = [rectangle, rounded corners, minimum width=2cm, minimum height=1cm, text centered, draw=black, fill=green!30]
arrow = [thick,->,>=stealth]

at (0,0) ( );
at (3,0) (start) [startstop] $\begin{array}{c}
\texttt{radiation\_generator}\\
\text{Construct real flavour structures}
\end{array}$;
(pro1) [process, below of=start] $\begin{array}{c}
\texttt{region\_data\_init}\\
\text{Initialise FKS regions}
\end{array}$;
(pro2) [OLP, right of=pro1, xshift=7cm] $\begin{array}{c}
\texttt{blha\_configuration\_write}\\
\text{Write \texttt{.olp} file}
\end{array}$;
(pro3) [OLP, below of=pro2] $\begin{array}{c}
\texttt{read\_contract\_file}\\
\text{Read \texttt{.olc} file}
\end{array}$;
(pro4) [OLP, below of=pro3] $\begin{array}{c}
\texttt{evaluate\_interaction}\\
\text{Evaluate sq. matrix elements}
\end{array}$;
(pro5) [process, below of=pro4, xshift=-7.5cm] $\begin{array}{c}
\texttt{process\_instance\_evaluate\_trace}\\
\text{Evaluate FKS terms}
\end{array}$;
(mci) [kinematics, below of=pro1, xshift=4cm,yshift=1.3cm]
(pro6) [kinematics, below of=pro1, xshift=3.5cm,yshift=-0.2cm] $\begin{array}{c}
\texttt{compute\_seed\_kinematics}\\
\text{Construct phase-space}
\end{array}$
(pro7) [kinematics, below of=pro6] $\begin{array}{c}
\text{Compute }S_{\tilde{\alpha}},\\
\text{cut expressions, etc.}
\end{array}$
[arrow] (start) – (pro1);
[arrow] (pro1) – (pro5);
[arrow] (pro1) – (pro2);
[arrow] (pro1) – (pro6);
[arrow] (pro6) – (pro7);
[arrow] (pro7) – (pro5);
[arrow] (pro2) – (pro3);
[arrow] (pro3) – (pro4);
[arrow] (pro4) – (pro5);
[arrow] (pro6) – (pro4);
[arrow] (mci) – (pro6);

Sketched workflow for NLO calculations in

Outgoing from the found emitter particles by initialising the FKS regions the radiation phase-space will be constructed. From the multi-channel integration objects  all ingredients for the construction of the phase-space such as random numbers are provided. In these random number variables the Born phase space $\bar{\Phi}_n$ as well as the FKS variables, $\xi$ and $y$, are parametrised as  is called. The NLO phase-space is thereby constructed for each index  associated with phase-space configurations distinguishing the Born-like points from those of the real emission and real phase-space configurations for distinct emitter positions [If radiation is possible from both,  and , the emitter position for an FKS region is replaced by . This is necessary for the phase-space construction in order to treat the initial state emitter symmetrically according to Sec. <ref> and <ref>].

The evaluation of matrix elements proceeds in the routine  with the input on the phase-space points per  and indices of the list of squared matrix elements requests from the  files linked to the  objects of the state matrix. After the  routines are called for the computation of squared amplitudes, their numerical values are stored as  array entries each linked to the  objects. Note, that for the case that  is linked to  the complete part of the code evaluation indicated in blue in Fig. <ref> is done in one step.

From the phase-space configurations all kinematical FKS factors entering the NLO calculation such as $S_{\tilde{\alpha}}$ as well as phase-space cut expressions are evaluated. These in addition to the previously stored squared matrix elements values are then used to evaluate the FKS terms by means of the constructed FKS singular regions. Precisely, in  all integrands represented by contributions of the Born in Eq. (<ref>), the real-subtracted in Eq. (<ref>), the virtual plus integrated subtraction terms in Eq. (<ref>) and the DGLAP remnant in Eq. (<ref>) are computed.

§.§ Basic EW extensions

An NLO EW calculation, either for integrated cross sections or event simulation, with respect to Eq. (<ref>) can be steered with  by running a SINDARIN file with commands of the usual NLO settings and additionally

\begin{align}
\texttt{\$nlo\_correction\_type = "EW"}
\label{SindarinEW}
\end{align}

For EW corrections all group theoretical factors entering the FKS terms have to correspond to those of QED. Due to a switch to  as correction type these are called in place of those corresponding to QCD, which is the default for NLO calculations in . This concerns the Altarelli-Parisi splitting kernels from Eqs. (<ref>) to (<ref>) and (<ref>) to (<ref>) and App. <ref>, Casimir operators $C^{\text{QED}}$ and definitions $\gamma^{\text{QED}}$ and $\gamma^{\prime \text{QED}}$ defined in Table <ref> for the collinear terms and charge correlations according to Eq. (<ref>) for the soft terms, respectively. In addition to these numerical factors the prescriptions of photons entering SM vertices have to be regarded throughout the construction of subtraction terms concerning the NLO EW automated corrections.