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, but tis constraint can however be lifted (see Scale Variations).

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

Leading Order

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.

LO Non-Singlet Evolution

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]\]

LO Singlet Evolution

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_{-})\]

Next-to-Leading Order

NLO Non-Singlet Evolution

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 \(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)\) and \(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)\)

  • 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]\]

NLO Singlet Evolution

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} 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\]

Next-to-Next-to-Leading Order

NNLO Non-Singlet Evolution

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^{(0)} \cdot j^{(0,2)}(a_s,a_s^0) + \gamma^{(1)} \cdot j^{(1,2)}(a_s,a_s^0) + \gamma^{(2)} \cdot 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 &= atan \left( \frac{b_1 + 2 a_s b_2 }{ \Delta} \right) - atan \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] \\ &= atan \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 pertubative 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^{(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]\]

NNLO Singlet Evolution

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-1} \,\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}\]

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.