Solving DGLAP

We are solving the DGLAP equations [AP77, Dok77, GL72] given in x-space by

\[\frac{d}{d\ln(\mu_F^2)} \mathbf{f}(x,\mu_F^2) = \int\limits_x^1\!\frac{dy}{y}\, \mathbf{P}(x/y,a_s) \cdot \mathbf{f}(y,\mu_F^2)\]

with \(\mathbf P\) the Altarelli-Parisi splitting functions (see pQCD ingredients). In Mellin space the DGLAP equations are just differential equations:

\[\frac{d}{d\ln(\mu_F^2)} \tilde{\mathbf{f}}(\mu_F^2) = -\gamma(a_s) \cdot \tilde{\mathbf{f}}(\mu_F^2)\]

(Note the additional minus in the definition for \(\gamma\)).

We change the evolution variable to the (monotonic) Strong Coupling \(a_s(\mu_F^2)\) and the equations to solve become

\[\frac{d}{da_s} \tilde{\mathbf{f}}(a_s) = \frac{d\ln(\mu_F^2)}{da_s} \cdot \frac{d \tilde{\mathbf{f}}(\mu_F^2)}{d\ln(\mu_F^2)} = -\frac{\gamma(a_s)}{\beta(a_s)} \cdot \tilde{\mathbf{f}}(a_s)\]

This assumes the factorization scale \(\mu_F^2\) (the inherit scale of the PDF) and the renormalization scale \(\mu_R^2\) (the inherit scale for the strong coupling) to be equal.

The (formal) solution can then be written in terms of an EKO \(\mathbf E\) [Bon12]

\[\begin{split}\tilde{\mathbf{f}}(a_s) &= \tilde{\mathbf{E}}(a_s \leftarrow a_s^0) \cdot \tilde{\mathbf{f}}(a_s^0)\\ \tilde{\mathbf{E}}(a_s \leftarrow a_s^0) &= \mathcal P \exp\left[-\int\limits_{a_s^0}^{a_s} \frac{\gamma(a_s')}{\beta(a_s')} da_s' \right]\end{split}\]

with \(\mathcal P\) the path-ordering operator. In the non-singlet sector the equations decouple and we do not need to worry about neither matrices nor the path-ordering.

Using Interpolation on both the initial and final PDF, we can then discretize the EKO in x-space and define \({\mathbf{E}}_{k,j}\) (represented by Operator) by

\[{\mathbf{E}}_{k,j}(a_s \leftarrow a_s^0) = \mathcal{M}^{-1}\left[\tilde{\mathbf{E}}(a_s \leftarrow a_s^0)\tilde p_j\right](x_k)\]

Now, we can write the solution to DGLAP in a true matrix operator scheme and find

\[\mathbf{f}(x_k,a_s) = {\mathbf{E}}_{k,j}(a_s \leftarrow a_s^0) \mathbf{f}(x_j,a_s^0)\]

so the EKO is a rank-4 operator acting both in flavor and momentum fraction space.

The issue of matching conditions when crossing flavor thresholds is discussed in a separate document

LO evolution

Expanding the anomalous dimension \(\gamma(a_s)\) and the beta function \(\beta(a_s)\) to LO we obtain the (exact) EKO:

\[\begin{split}\ln \tilde {\mathbf E}^{(0)}(a_s \leftarrow a_s^0) &= \gamma^{(0)}\int\limits_{a_s^0}^{a_s} \frac{da_s'}{\beta_0 a_s'} = \gamma^{(0)} \cdot j^{(0,0)}(a_s,a_s^0)\\ j^{(0,0)}(a_s,a_s^0) &= \int\limits_{a_s^0}^{a_s} \frac{da_s'}{\beta_0 a_s'} = \frac{\ln(a_s/a_s^0)}{\beta_0}\end{split}\]

In LO we always use the exact solution.

Non-Singlet

We find

\[\frac{d}{da_s} \tilde f_{ns}^{(0)}(a_s) = \frac{\gamma_{ns}^{(0)}}{\beta_0 a_s} \cdot \tilde f_{ns}^{(0)}(a_s)\]

with \(\gamma_{ns}^{(0)} = \gamma_{ns,+}^{(0)} = \gamma_{ns,-}^{(0)} = \gamma_{ns,v}^{(0)} = \gamma_{qq}^{(0)}\).

The EKO is then given by a simple exponential [Vog05]

\[\tilde E^{(0)}_{ns}(a_s \leftarrow a_s^0) = \exp\left[\gamma_{ns}^{(0)} \ln(a_s/a_s^0)/\beta_0 \right]\]

Singlet

We find

\[\begin{split}\frac{d}{da_s} \dSV{0}{a_s} = \frac{\gamma_S^{(0)}}{\beta_0 a_s} \cdot \dSV{0}{a_s}\,, \qquad \gamma_S^{(0)} = \begin{pmatrix} \gamma_{qq}^{(0)} & \gamma_{qg}^{(0)}\\ \gamma_{gq}^{(0)} & \gamma_{gg}^{(0)} \end{pmatrix}\end{split}\]

In order to exponentiate the EKO, we decompose it \(\ln \mathbf{\tilde E}^{(0)}_S = \lambda_+ {\mathbf e}_+ + \lambda_- {\mathbf e}_-\) with the eigenvalues \(\lambda_{\pm}\) and the projectors \(\mathbf e_{\pm}\) given by [Vog05]

\[\begin{split}\lambda_{\pm} &= \frac 1 {2} \left( \ln \tilde E_{qq}^{(0)} + \ln \tilde E_{gg}^{(0)} \pm \sqrt{(\ln \tilde E_{qq}^{(0)}-\ln \tilde E_{gg}^{(0)})^2 + 4\ln \tilde E_{qg}^{(0)}\ln \tilde E_{gq}^{(0)}} \right)\\ {\mathbf e}_{\pm} &= \frac{1}{\lambda_{\pm} - \lambda_{\mp}} \left( \ln \mathbf{\tilde E}^{(0)}_S - \lambda_{\mp} \mathbf I \right)\end{split}\]

with \(\mathbf I\) the 2x2 identity matrix in flavor space and, e.g., \(\ln \tilde E_{qq}^{(0)} = \gamma_{qq}^{(0)}j^{(0,0)}(a_s,a_s^0)\).

The projectors obey the usual properties, i.e.

\[{\mathbf e}_{\pm} \cdot {\mathbf e}_{\pm} = {\mathbf e}_{\pm}\,,\quad {\mathbf e}_{\pm} \cdot {\mathbf e}_{\mp} = 0\,,\quad \ep + \em = \mathbf I\]

and thus the exponentiation becomes easier again.

The EKO is then given by

\[\ESk{0}{a_s}{a_s^0} = \ep \exp(\lambda_{+}) + \em \exp(\lambda_{-})\]

NLO evolution

Non-Singlet

We find

\[\frac{d}{da_s} \tilde f_{ns}^{(1)}(a_s) = \frac{\gamma_{ns}^{(0)} a_s + \gamma_{ns}^{(1)} a_s^2}{\beta_0 a_s^2 + \beta_1 a_s^3} \cdot \tilde f_{ns}^{(1)}(a_s)\]

with \(\gamma_{ns} \in \{\gamma_{ns,+},\gamma_{ns,-}=\gamma_{ns,v}\}\).

We obtain the (exact) EKO [Bon12, RA98, Vog05]:

\[\begin{split}\ln \tilde E^{(1)}_{ns}(a_s \leftarrow a_s^0) &= \gamma^{(0)} \cdot j^{(0,1)}(a_s,a_s^0) + \gamma^{(1)} \cdot j^{(1,1)}(a_s,a_s^0)\\ j^{(1,1)}(a_s,a_s^0) &= \int\limits_{a_s^0}^{a_s}\!da_s'\,\frac{a_s'^2}{\beta_0 a_s'^2 + \beta_1 a_s'^3} = \frac{1}{\beta_1}\ln\left(\frac{1+b_1 a_s}{1+b_1 a_s^0}\right)\\ j^{(0,1)}(a_s,a_s^0) &= \int\limits_{a_s^0}^{a_s}\!da_s'\,\frac{a_s'}{\beta_0 a_s'^2 + \beta_1 a_s'^3} = j^{(0,0)}(a_s,a_s^0) - b_1 j^{(1,1)}(a_s,a_s^0)\end{split}\]

Note that we recover the LO solution:

\[\ln \tilde E^{(1)}_{ns}(a_s \leftarrow a_s^0) = \ln \tilde E^{(0)}_{ns}(a_s \leftarrow a_s^0) + j^{(1,1)}(a_s,a_s^0)(\gamma^{(1)} - b_1 \gamma^{(0)})\]

In NLO we provide different strategies to define the EKO:

  • method in ['iterate-exact', 'decompose-exact', 'perturbative-exact']: use the exact solution as defined above

  • method in ['iterate-expanded', 'decompose-expanded', 'perturbative-expanded']: use the exact LO solution and substitute:

    \[\begin{split}j^{(1,1)}(a_s,a_s^0) \to j^{(1,1)}_{exp}(a_s,a_s^0) &= \frac 1 {\beta_0}(a_s - a_s^0) \\ j^{(0,1)}(a_s,a_s^0) \to j^{(0,1)}_{exp}(a_s,a_s^0) &= j^{(0,0)}(a_s,a_s^0) - b_1 j^{(1,1)}_{exp}(a_s,a_s^0) \\\end{split}\]

  • method = 'ordered-truncated': expanding the argument of the exponential of the new term but keeping the order we obtain:

\[\tilde E^{(1)}_{ns}(a_s \leftarrow a_s^0) = \tilde E^{(0)}_{ns}(a_s \leftarrow a_s^0) \frac{1 + a_s/\beta_0 (\gamma_{ns}^{(1)} - b_1 \gamma_{ns}^{(0)})}{1 + a_s^0/\beta_0 (\gamma_{ns}^{(1)} - b_1 \gamma_{ns}^{(0)})}\]

  • method = 'truncated': expanding the whole exponential of the new term we obtain:

\[\tilde E^{(1)}_{ns}(a_s \leftarrow a_s^0) = \tilde E^{(0)}_{ns}(a_s \leftarrow a_s^0) \left[1 + (a_s - a_s^0)/\beta_0 (\gamma_{ns}^{(1)} - b_1 \gamma_{ns}^{(0)}) \right]\]

Singlet

We find

\[\frac{d}{da_s} \dSV{1}{a_s} = \frac{\gamma_{S}^{(0)} a_s + \gamma_{S}^{(1)} a_s^2}{\beta_0 a_s^2 + \beta_1 a_s^3} \cdot \dSV{1}{a_s}\]

with \(\gamma_{S}^{(0)} \gamma_{S}^{(1)} \neq \gamma_{S}^{(1)} \gamma_{S}^{(0)}\).

Here the strategies are:

  • for method in ['iterate-exact', 'iterate-expanded'] we use a discretized path-ordering [Bon12]:

\[\ESk{1}{a_s}{a_s^0} = \prod\limits_{k=n}^{0} \ESk{1}{a_s^{k+1}}{a_s^{k}}\quad \text{with}\quad a_s^{n+1} = a_s\]

where the order of the product is such that later EKO are to the left and

\[\begin{split}\ESk{1}{a_s^{k+1}}{a_s^{k}} &= \exp\left(-\frac{\gamma(a_s^{k+1/2})}{\beta(a_s^{k+1/2})} \Delta a_s \right) \\ a_s^{k+1/2} &= a_0 + \left(k+ \frac 1 2\right) \Delta a_s\\ \Delta a_s &= \frac{a_s - a_s^0}{n + 1}\end{split}\]

using the projector algebra from LO to exponentiate the single steps.

  • for method in ['decompose-exact', 'decompose-expanded']: use the exact or the approximate exact integrals from the non-singlet sector and then decompose \(\ln \tilde{\mathbf E}^{(1)}\) - this will neglect the non-commutativity of the singlet matrices.

  • for method in ['perturbative-exact', 'perturbative-expanded', 'ordered-truncated', 'truncated'] we seek for an perturbative solution around the (exact) leading order operator:

    We set [Vog05]

    \[\frac{d}{da_s} \dSV{1}{a_s} = \frac{\mathbf R (a_s)}{a_s} \cdot \dSV{1}{a_s}\,, \quad \mathbf R (a_s) = \sum\limits_{k=0} a_s^k \mathbf R_{k}\]

    where in NLO we find

    \[\mathbf R_0 = \gamma_{S}^{(0)}/\beta_0\,,\quad \mathbf R_1 = \gamma_{S}^{(1)}/\beta_0 - b_1 \gamma_{S}^{(0)} /\beta_0\]

    and for the higher coefficients

    • method = 'perturbative-exact': \(\mathbf R_k = - b_1 \mathbf R_{k-1}\,\text{for}\,k>1\)

    • method = 'perturbative-expanded': \(\mathbf R_k = 0\,\text{for}\,k>1\)

    We make an ansatz for the solution

    \[\ESk{1}{a_s}{a_s^0} = \mathbf U (a_s) \ESk{0}{a_s}{a_s^0} {\mathbf U}^{-1} (a_s^0), \quad \mathbf U (a_s) = \mathbf I + \sum\limits_{k=1} a_s^k \mathbf U_k\]

    Inserting this ansatz into the differential equation and sorting by powers of \(a_s\), we obtain a recursive set of commutator relations for the evolution operator coefficients \(\mathbf U_k\):

    \[\begin{split}[\mathbf U_1, \mathbf R_0] &= \mathbf R_1 - \mathbf U_1\\ [\mathbf U_k, \mathbf R_0] &= \mathbf R_k + \sum\limits_{j=1}^{k-1} \mathbf R_{k-j} \mathbf U_j - k \mathbf U_k = \mathbf{R}_k' - k \mathbf U_k\,,k>1\end{split}\]

    Multiplying these equations with \(\mathbf e_{\pm}\) from left and right and using the identity

    \[\mathbf U_k = \em \mathbf U_k \em + \em \mathbf U_k \ep + \ep \mathbf U_k \em + \ep \mathbf U_k \ep\]

    we obtain the \(\mathbf U_k\)

    \[\mathbf U_k = \frac{ \em \mathbf{R}_k' \em + \ep \mathbf{R}_k' \ep } k + \frac{\ep \mathbf{R}_k' \em}{r_- - r_+ + k} + \frac{\em \mathbf{R}_k' \ep}{r_+ - r_- + k}\]

    with \(r_{\pm} =\frac 1 {2\beta_0} \left( \gamma_{qq}^{(0)} + \gamma_{gg}^{(0)} \pm \sqrt{(\gamma_{qq}^{(0)}-\gamma_{gg}^{(0)})^2 + 4\gamma_{qg}^{(0)}\gamma_{gq}^{(0)}} \right)\).

    So the strategies are

    • method in ['perturbative-exact', 'perturbative-expanded']: approximate the full evolution operator \(\mathbf U(a_s)\) with an expansion up to ev_op_max_order

    • method in ['ordered-truncated', 'truncated']: truncate the evolution operator \(\mathbf U(a_s)\) and use

    \[\ESk{1}{a_s}{a_s^0} = \ESk{0}{a_s}{a_s^0} + a_s \mathbf U_1 \ESk{0}{a_s}{a_s^0} - a_s^0 \ESk{0}{a_s}{a_s^0} \mathbf U_1\]

NNLO evolution

Non-Singlet

We find

\[\frac{d}{da_s} \tilde f_{ns}^{(2)}(a_s) = \frac{\gamma_{ns}^{(0)} a_s + \gamma_{ns}^{(1)} a_s^2 + \gamma_{ns}^{(2)} a_s^3 }{\beta_0 a_s^2 + \beta_1 a_s^3 + \beta_2 a_s^4} \cdot \tilde f_{ns}^{(2)}(a_s)\]

with \(\gamma_{ns} \in \{\gamma_{ns,+},\gamma_{ns,-},\gamma_{ns,v}\}\).

We obtain the (exact) EKO [CCG08, Vog05]:

\[\begin{split}\ln \tilde E^{(2)}_{ns}(a_s \leftarrow a_s^0) &= \gamma_{ns}^{(0)} j^{(0,2)}(a_s,a_s^0) + \gamma_{ns}^{(1)} j^{(1,2)}(a_s,a_s^0) + \gamma_{ns}^{(2)} j^{(2,2)}(a_s,a_s^0)\\\end{split}\]

with:

\[\begin{split}j^{(2,2)}(a_s,a_s^0) &= \int\limits_{a_s^0}^{a_s}\!da_s'\,\frac{a_s'^3}{\beta_0 a_s'^2 + \beta_1 a_s'^3 + \beta_2 a_s'^4} = \frac{1}{\beta_2}\ln\left(\frac{1 + a_s ( b_1 + b_2 a_s ) }{ 1 + a_s^0 ( b_1 + b_2 a_s^0 )}\right) - \frac{b_1}{ \beta_2 \Delta} \delta \\ \delta &= \arctan \left( \frac{b_1 + 2 a_s b_2 }{ \Delta} \right) - \arctan \left( \frac{b_1 + 2 a_s^0 b_2 }{ \Delta} \right) \\ &= \frac{i}{2} \left[ ln \left( \frac{ \Delta - i (b_1 + 2a_s b_2)}{ \Delta + i (b_1 + 2a_s b_2)}\right) - ln \left( \frac{ \Delta - i (b_1 + 2a_s^0 b_2)}{ \Delta + i (b_1 + 2a_s^0 b_2)}\right) \right] \\ &= \arctan \left( \frac{\Delta ( a_s - a_s^0 )}{ 2 + b_1 (a_s + a_s^0) + 2 a_s a_s^0 b_2 } \right) \\ \Delta &= \sqrt{4 b_2 - b_1^2 }\end{split}\]

and:

\[\begin{split}j^{(1,2)}(a_s,a_s^0) &= \int\limits_{a_s^0}^{a_s}\!da_s'\,\frac{a_s'^2}{\beta_0 a_s'^2 + \beta_1 a_s'^3 + \beta_2 a_s'^4} = \frac{2}{\beta_0 \Delta} \delta \\ j^{(0,2)}(a_s,a_s^0) &= \int\limits_{a_s^0}^{a_s}\!da_s'\,\frac{a_s'}{\beta_0 a_s'^2 + \beta_1 a_s'^3 + \beta_2 a_s'^4} = j^{(0,0)}(a_s,a_s^0) - b_1 j^{(1,2)}(a_s,a_s^0) - b_2 j^{(2,2)}(a_s,a_s^0)\end{split}\]

Note, plugging the numerical values of \(\beta_i\) we find that the \(\Delta \in \mathbb{R}\) if \(n_f < 6\). However you can notice that \(\Delta\) appears always with \(\delta\) and the fraction \(\frac{\delta}{\Delta} \in \mathbb{R}, \forall n_f\).

We can recover the LO solution:

\[\ln \tilde E^{(2)}_{ns}(a_s \leftarrow a_s^0) = \ln \tilde E^{(0)}_{ns}(a_s \leftarrow a_s^0) + j^{(1,2)}(a_s,a_s^0)(\gamma^{(1)} - b_1 \gamma^{(0)}) + j^{(2,2)}(a_s,a_s^0)(\gamma^{(2)} - b_2 \gamma^{(0)})\]

And thus the NLO solution:

\[\begin{split}\ln \tilde E^{(2)}_{ns}(a_s \leftarrow a_s^0) &= \ln \tilde E^{(1)}_{ns}(a_s \leftarrow a_s^0) + j^{(1,2)'}(a_s,a_s^0)(\gamma^{(1)} - b_1 \gamma^{(0)}) + j^{(2,2)}(a_s,a_s^0)(\gamma^{(2)} - b_2 \gamma^{(0)}) \\ j^{(1,2)'}(a_s,a_s^0) &= \int\limits_{a_s^0}^{a_s}\!da_s'\,\frac{ \beta_2 a_s'^2}{( \beta_0 + \beta_1 a_s' + \beta_2 a_s'^2 ) (\beta_0 + \beta_1 a_s')}\end{split}\]

In NNLO we provide different strategies to define the EKO:

  • method in ['iterate-exact', 'decompose-exact', 'perturbative-exact']: use the exact solution as defined above

  • method in ['iterate-expanded', 'decompose-expanded', 'perturbative-expanded']: use the exact LO solution and expand all functions \(j^{(n,m)}(a_s,a_s^0)\) to the order \(\mathcal o(a_s^3)\). We find:

\[\begin{split}j^{(2,2)}(a_s,a_s^0) \approx j^{(2,2)}_{exp}(a_s,a_s^0) &= \frac{1}{2\beta_0} (a_s^2 - (a_s^0)^{2}) \\ j^{(1,2)}(a_s,a_s^0) \approx j^{(1,2)}_{exp}(a_s,a_s^0) &= \frac{1}{\beta_0} [ (a_s - a_s^0) - \frac{b_1}{2} (a_s^2 - (a_s^0)^{2})] \\ j^{(0,2)}(a_s,a_s^0) \approx j^{(0,2)}_{exp}(a_s,a_s^0) &= j^{(0,0)}(a_s,a_s^0) - b_1 j^{(1,2)}_{exp}(a_s,a_s^0) - b_2 j^{(2,2)}_{exp}(a_s,a_s^0) \\ &= j^{(0,0)}(a_s,a_s^0) - \frac{1}{\beta_0} [ b_1 (a_s - a_s^0) + \frac{b_1^2-b_2}{2} (a_s^2 - (a_s^0)^{2}) ] \\\end{split}\]

This method corresponds to IMODEV=2 of [Vog05].

  • method = 'ordered-truncated': for this method we follow the prescription from [Vog05] and we get:

\[\tilde E^{(2)}_{ns}(a_s \leftarrow a_s^0) = \tilde E^{(0)}_{ns}(a_s \leftarrow a_s^0) \frac{ 1 + a_s U_1 + a_s^2 U_2 }{ 1 + a_s^0 U_1 + (a_s^0)^{2} U_2 }\]

with the unitary matrices defined consistently with the method perturbative adopted for NLO singlet evolution:

\[\begin{split}U_1 &= R_1 = \frac{1}{\beta_0}[ \gamma^{(1)} - b_1 \gamma^{(0)}] \\ U_2 &= \frac{1}{2}[ R_1^2 + R_2 ] \\ R_2 &= \gamma_{ns}^{(2)}/\beta_0 - b_1 R_1 - b_2 R_0 \\\end{split}\]

This method corresponds to IMODEV=3 of [Vog05].

  • method = 'truncated': we expand the whole exponential and keeping terms within \(\mathcal o(a_s^3)\). This method is the fastest among the ones provided by our program. We obtain:

\[\tilde E^{(2)}_{ns}(a_s \leftarrow a_s^0) = \tilde E^{(0)}_{ns}(a_s \leftarrow a_s^0) \left [ 1 + U_1 (a_s - a_s^0) + a_s^2 U_2 - a_s a_s^0 U_1^2 + (a_s^0)^{2} ( U_1^2 - U_2 ) \right]\]

Singlet

For the singlet evolution we find:

\[\frac{d}{da_s} \dSV{2}{a_s} = \frac{\gamma_{S}^{(0)} a_s + \gamma_{S}^{(1)} a_s^2 + \gamma_{S}^{(2)} a_s^3}{\beta_0 a_s^2 + \beta_1 a_s^3 + \beta_2 a_s^4} \cdot \dSV{2}{a_s}\]

with \(\gamma_{S}^{(i)} \gamma_{S}^{(j)} \neq \gamma_{S}^{(j)} \gamma_{S}^{(i)}, \quad i,j=0,1,2\).

In analogy to NLO we define the following strategies :

  • for method in ['iterate-exact', 'iterate-expanded'] we use a discretized path-ordering [Bon12]:

\[\ESk{2}{a_s}{a_s^0} = \prod\limits_{k=n}^{0} \ESk{2}{a_s^{k+1}}{a_s^{k}} \quad \text{with} \quad a_s^{n+1} = a_s\]

All the procedure is identical to NLO, simply the beat function is now expanded until \(\mathcal o(a_s^4)\)

  • for method in ['decompose-exact', 'decompose-expanded']: use the exact or the approximate exact integrals from the non-singlet sector and then decompose \(\ln \tilde{\mathbf E}^{(2)}\) - this will neglect the non-commutativity of the singlet matrices.

  • for method in ['perturbative-exact', 'perturbative-expanded', 'ordered-truncated', 'truncated'] we seek for an perturbative solution around the (exact) leading order operator. We set [Vog05]

    \[\frac{d}{da_s} \dSV{2}{a_s} = \frac{\mathbf R (a_s)}{a_s} \cdot \dSV{2}{a_s}\,, \quad \mathbf R (a_s) = \sum\limits_{k=0} a_s^k \mathbf R_{k}\]

    Finding one additional term compared to NLO:

    \[\begin{split}\mathbf R_2 & = \gamma_{S}^{(2)}/\beta_0 - b_1 \mathbf R_1 - b_2 \mathbf R_0 \\ & = \frac{1}{\beta_0} [ \gamma_{S}^{(2)} - b_1 \gamma_{S}^{(1)} - \gamma_{S}^{(0)} ( b_2 - b_1^2 ) ]\end{split}\]

    and for the higher coefficients

    • method = 'perturbative-exact': \(\mathbf R_k = - b_1 \mathbf R_{k-1} - b_2 \mathbf R_{k-2} \,\text{for}\,k>2\)

    • method = 'perturbative-expanded': \(\mathbf R_k = 0\,\text{for}\,k>2\)

    The solution ansatz becomes:

    \[\ESk{2}{a_s}{a_s^0} = \mathbf U (a_s) \ESk{0}{a_s}{a_s^0} {\mathbf U}^{-1} (a_s^0), \quad \mathbf U (a_s) = \mathbf I + \sum\limits_{k=1} a_s^k \mathbf U_k\]

    with:

    \[\begin{split}[\mathbf U_2, \mathbf R_0] &= \mathbf R_2 + \mathbf R_1 \mathbf U_1 - 2 \mathbf U_2\\\end{split}\]

    So the strategies are:

    • method in ['perturbative-exact', 'perturbative-expanded']: approximate the full evolution operator \(\mathbf U(a_s)\) with an expansion up to ev_op_max_order

    • method in ['ordered-truncated', 'truncated']: truncate the evolution operator \(\mathbf U(a_s)\) and use

    \[\begin{split}\ESk{2}{a_s}{a_s^0} &= \ESk{0}{a_s}{a_s^0} + a_s \mathbf U_1 \ESk{0}{a_s}{a_s^0} - a_s^0 \ESk{0}{a_s}{a_s^0} \mathbf U_1 \\ &\hspace{20pt} + a_s^2 \mathbf U_2 \ESk{0}{a_s}{a_s^0} \\ &\hspace{20pt} + a_s a_s^0 \mathbf U_1 \ESk{0}{a_s}{a_s^0} \mathbf U_1 \\ &\hspace{20pt}- (a_s^0)^{2} \ESk{0}{a_s}{a_s^0} ( \mathbf U_1^2 - \mathbf U_2 )\end{split}\]

N3LO evolution

Non-Singlet

At N3LO the DGLAP expansion reads:

\[\frac{d}{da_s} \tilde f_{ns}^{(2)}(a_s) = \frac{\gamma_{ns}^{(0)} a_s + \gamma_{ns}^{(1)} a_s^2 + \gamma_{ns}^{(2)} a_s^3 + \gamma_{ns}^{(3)} a_s^4 }{\beta_0 a_s^2 + \beta_1 a_s^3 + \beta_2 a_s^4 + \beta_3 a_s^5 } \cdot \tilde f_{ns}^{(2)}(a_s)\]

with \(\gamma_{ns} \in \{\gamma_{ns,+},\gamma_{ns,-},\gamma_{ns,v}\}\).

We obtain the (exact) EKO in analogy of the previous orders and recovering the LO solution:

\[\begin{split}\ln \tilde E^{(2)}_{ns}(a_s \leftarrow a_s^0) &= \gamma_{ns}^{(0)} j^{(0,3)}(a_s,a_s^0) + \gamma_{ns}^{(1)} j^{(1,3)}(a_s,a_s^0) + \gamma_{ns}^{(2)} j^{(2,3)}(a_s,a_s^0) + \gamma_{ns}^{(3)} j^{(3,3)}(a_s,a_s^0)\\\end{split}\]

with:

\[\begin{split}j^{(3,3)}(a_s,a_s^0) &= \frac{1}{\beta_0} \sum_{r=r_1}^{r_3} \frac{ r^2 [\ln(a_s-r) - \ln(a_s^0 - r)]}{b_1 + 2 b_2 r + 3 b_3 r^2} \\ j^{(2,3)}(a_s,a_s^0) &= \frac{1}{\beta_0} \sum_{r=r_1}^{r_3} \frac{r [\ln(a_s-r) - \ln(a_s^0 - r)]}{b_1 + 2 b_2 r + 3 b_3 r^2} \\ j^{(1,3)}(a_s,a_s^0) &= \frac{1}{\beta_0} \sum_{r=r_1}^{r_3} \frac{\ln(a_s-r) - \ln(a_s^0 - r)}{b_1 + 2 b_2 r + 3 b_3 r^2} \\ j^{(0,3)}(a_s,a_s^0) &= j^{(0,0)}(a_s,a_s^0) - b_1 j^{(1,3)}(a_s,a_s^0) - b_2 j^{(2,3)}(a_s,a_s^0)- b_3 j^{(3,3)}(a_s,a_s^0)\end{split}\]

where the sum is carried on the complex roots of the beta function expansion:

\[a_s \in \{r_1, r_2, r_3 \} | \quad 1 + b_1 a_s + b_2 a_s^2 + b_3 a_s^3 = 0\]

You can notice that in the denominator of the integrals appears always the derivative of this expansion. We remark that even though the roots are complex the total integral is real.

Also in this case we provide a we provide different strategies to define the EKO:

  • method in ['iterate-exact', 'decompose-exact', 'perturbative-exact']: use the exact solution as defined above

  • method in ['iterate-expanded', 'decompose-expanded', 'perturbative-expanded']: use the exact LO solution and expand all functions \(j^{(n,m)}(a_s,a_s^0)\) to the order \(\mathcal o(a_s^3)\). We find:

\[\begin{split}j^{(3,3)}(a_s,a_s^0) &\approx j^{(3,3)}_{exp}(a_s,a_s^0) = \frac{1}{3 \beta_0} (a_s^3 - (a_s^0)^3) \\ j^{(2,3)}(a_s,a_s^0) &\approx j^{(2,3)}_{exp}(a_s,a_s^0) = \frac{1}{\beta_0} [ \frac{1}{2} (a_s^2 - (a_s^0)^2) - \frac{b_1}{3} (a_s^3 - (a_s^0)^3) ]\\ j^{(1,3)}(a_s,a_s^0) &\approx j^{(1,3)}_{exp}(a_s,a_s^0) = \frac{1}{\beta_0} [ (a_s - a_s^0) - \frac{b_1}{2} (a_s^2 - (a_s^0)^2) + \frac{b_1^2-b_2}{3} (a_s^3 - (a_s^0)^3) ]\\ j^{(0,2)}(a_s,a_s^0) &\approx j^{(0,3)}_{exp}(a_s,a_s^0) = j^{(0,0)}(a_s,a_s^0) - b_1 j^{(1,3)}_{exp}(a_s,a_s^0) - b_2 j^{(2,3)}_{exp}(a_s,a_s^0)- b_3 j^{(3,3)}_{exp}(a_s,a_s^0) \\\end{split}\]

  • method = 'ordered-truncated': performing the expansion one order higher wrt to NNLO we get:

\[\tilde E^{(2)}_{ns}(a_s \leftarrow a_s^0) = \tilde E^{(0)}_{ns}(a_s \leftarrow a_s^0) \frac{ 1 + a_s U_1 + a_s^2 U_2 + a_s^3 U_3 }{ 1 + a_s^0 U_1 + (a_s^0)^{2} U_2 + (a_s^0)^{3} U_3 }\]

with the new unitary matrices defined:

\[\begin{split}U_3 &= \frac{1}{3} [R_3 + R_2 U_1 + R_1 U_2] \\ R_3 &= \gamma_{ns}^{(3)}/\beta_0 - b_1 R_2 - b_2 R_1 - b_3 R_0 \\\end{split}\]

  • method = 'truncated':

\[\begin{split}\tilde E^{(2)}_{ns}(a_s \leftarrow a_s^0) = \tilde E^{(0)}_{ns}(a_s \leftarrow a_s^0) & \left[ \right. 1 \\ & + U_1 (a_s - a_s^0) \\ & + a_s^2 U_2 - a_s a_s^0 U_1^2 + (a_s^0)^{2} ( U_1^2 - U_2 ) \\ & + a_s^3 U_3 - a_s^2 a_s^0 U_2 U_1 + a_s (a_s^0)^{2} U_1 ( U_1^2 - U_2 ) - (a_s^0)^{3} (U_1^3 - 2 U_1 U_2 + U_3) \left. \right] \\\end{split}\]

Singlet

For the singlet evolution we find:

\[\frac{d}{da_s} \dSV{2}{a_s} = \frac{\gamma_{S}^{(0)} a_s + \gamma_{S}^{(1)} a_s^2 + \gamma_{S}^{(2)} a_s^3 + \gamma_{S}^{(3)} a_s^3}{\beta_0 a_s^2 + \beta_1 a_s^3 + \beta_2 a_s^4 + \beta_3 a_s^5} \cdot \dSV{2}{a_s}\]

with \(\gamma_{S}^{(i)} \gamma_{S}^{(j)} \neq \gamma_{S}^{(j)} \gamma_{S}^{(i)}, \quad i,j=0,1,2,3\).

In analogy to NLO we define the following strategies :

  • for method in ['iterate-exact', 'iterate-expanded']: the solution strategies is exactly the same as in NLO and NNLO simply the beat function is now expanded until \(\mathcal o(a_s^5)\)

  • for method in ['decompose-exact', 'decompose-expanded']: use the exact or the approximate exact integrals from the non-singlet sector and then decompose \(\ln \tilde{\mathbf E}^{(3)}\) - this will neglect the non-commutativity of the singlet matrices.

  • for method in ['perturbative-exact', 'perturbative-expanded', 'ordered-truncated', 'truncated'] we seek for an perturbative solution around the (exact) leading order operator. Following the notation used for previous orders we have:

    \[\begin{split}\mathbf R_2 & = \gamma_{S}^{(3)}/\beta_0 - b_1 \mathbf R_2 - b_2 \mathbf R_1 - b_3 \mathbf R_0 \\\end{split}\]

    and for the higher coefficients:

    • method = 'perturbative-exact': \(\mathbf R_k = - b_1 \mathbf R_{k-1} - b_2 \mathbf R_{k-2} - b_3 \mathbf R_{k-3} \,\text{for}\,k>3\)

    • method = 'perturbative-expanded': \(\mathbf R_k = 0\,\text{for}\,k>3\)

    The new unitary matrix entering in the evolution ansatz follows the commutation relation:

    \[\begin{split}[\mathbf U_3, \mathbf R_0] &= \mathbf R_3 + \mathbf R_2 \mathbf U_1 + \mathbf R_1 \mathbf U_2 - 3 \mathbf U_3\\\end{split}\]

    So the strategies are:

    • method in ['perturbative-exact', 'perturbative-expanded']: approximate the full evolution operator \(\mathbf U(a_s)\) with an expansion up to ev_op_max_order

    • method in ['ordered-truncated', 'truncated']: truncate the evolution operator \(\mathbf U(a_s)\) and use

    \[\begin{split}\ESk{3}{a_s}{a_s^0} = \mathbf E^{0} &+ a_s \mathbf U_1 \mathbf E^{0} - a_s^0 \mathbf E^{0} \mathbf U_1 \\ & + a_s^2 \mathbf U_2 \mathbf E^{0} + a_s a_s^0 \mathbf U_1 \mathbf E^{0} \mathbf U_1 - (a_s^0)^{2} \mathbf E^{0} ( \mathbf U_1^2 - \mathbf U_2 ) \\ & + a_s^3 \mathbf U_3 \mathbf E^{0} - a_s^2 a_s^0 \mathbf U_2 \mathbf E^{0} \mathbf U_1 + a_s (a_s^0)^2 \mathbf U_1 E^{0} (\mathbf U_1^2 - \mathbf U_2) \\ & - (a_s^0)^3 \mathbf E^{0} (\mathbf U_1^3 - \mathbf U_1 \mathbf U_2 - \mathbf U_2 \mathbf U_1 + \mathbf U_3) \\\end{split}\]
    \[\mathbf E^{0} = \ESk{0}{a_s}{a_s^0}\]

Intrinsic evolution

We also consider the evolution of intrinsic heavy PDF. Since these are massive partons they can not split any collinear particles and thus they do not participate in the DGLAP evolution. Instead, their evolution is simply an identity operation: e.g. for an intrinsic charm distribution we get for \(m_c^2 > Q_1^2 > Q_0^2\):

\[\begin{split}\tilde c(Q_1^2) &= \tilde c(Q_0^2)\\ \tilde {\bar c}(Q_1^2) &= \tilde{\bar c}(Q_0^2)\end{split}\]

After crossing the mass threshold (charm in this example) the PDF can not be considered intrinsic any longer and hence, they have to be rejoined with their evolution basis elements and take then again part in the ordinary collinear evolution.

Mixed QCD \(\otimes\) QED evolution

For the moment in this case only the exact evolution is implemented.

Singlet

The evolution is obtained in the same way of the pure QCD case, with the only difference that now both \(\gamma\) and \(\beta_{qcd}\) contain the QED corrections.

In the case in which \(\alpha_{em}\) is running, at every step of the iteration the corresponding value of \(a_{em}(a_s)\) is used.

Non singlet

For the non singlet, being it diagonal, the solution is straightforward. When \(\alpha_{em}\) is fixed, the terms proportional to it are just a constant in the splitting functions, and therefore they can be integrated directly. For example at order=(1,1) we have

\[\begin{split}\tilde E^{(1,1)}_{ns}(a_s \leftarrow a_s^0) &= \exp \Bigl( -\int_{\log \mu_0^2}^{\log \mu^2}d\log\mu^2 \gamma_{ns}^{(1,0)} a_s(\log\mu^2) + \gamma_{ns}^{(1,1)} a_s(\log\mu^2) a_{em} + \gamma_{ns}^{(0,1)} a_em \Bigr) \\ & = \exp \Bigl( \int_{a_s^0}^{a_s}da_s\frac{\gamma_{ns}^{(1,0)} a_s + \gamma_{ns}^{(1,1)} a_s a_{em} + \gamma_{ns}^{(0,1)} a_em}{a_s^2(\beta_0 + \beta_0^{mix} a_{em})} -\int_{\log \mu_0^2}^{\log \mu^2}d\log\mu^2 \gamma_{ns}^{(0,1)} a_em\Bigr)\end{split}\]

In the last expression, the first term can be integrated with the \(j^{(n,m)\) functions, while the second term is trivial.

In the case of \(\alpha_{em}\) running, the \(a_s\) integration integral is divided in steps, such that in every step \(\alpha_{em}\) is considered constant. In this way the solution will be the product of the solutions of every integration step:

\[\tilde E^{(1,1)}_{ns}(a_s \leftarrow a_s^0) = \prod\limits_{k=n}^{0} E^{(1,1)}_{ns}(a_s^{k+1} \leftarrow a_s^k, a_{em}^k)\]