diff --git a/crates/eko/src/bib.rs b/crates/eko/src/bib.rs
new file mode 120000
index 000000000..1e713280e
--- /dev/null
+++ b/crates/eko/src/bib.rs
@@ -0,0 +1 @@
+../../ekore/src/bib.rs
\ No newline at end of file
diff --git a/crates/eko/src/lib.rs b/crates/eko/src/lib.rs
index a5dc15a6a..67bf8956a 100644
--- a/crates/eko/src/lib.rs
+++ b/crates/eko/src/lib.rs
@@ -3,7 +3,8 @@
 use ekore::harmonics::cache::Cache;
 use std::ffi::c_void;
 
-mod mellin;
+pub mod bib;
+pub mod mellin;
 
 /// QCD intergration kernel inside quad.
 ///
diff --git a/crates/eko/src/mellin.rs b/crates/eko/src/mellin.rs
index f2f968fe4..1322bbe08 100644
--- a/crates/eko/src/mellin.rs
+++ b/crates/eko/src/mellin.rs
@@ -9,10 +9,13 @@ use std::f64::consts::PI;
 /// Talbot inversion path.
 ///
 /// Implements the algorithm presented in [\[Abate\]](crate::bib::Abate).
-/// $p_{\text{Talbot}}(t) =  o + r \cdot ( \theta \cot(\theta) + i\theta)$ with $\theta = \pi(2t-1)$
-/// The default values for the parameters $r,o$ are given by $r = 1/2, o = 0$ for
-/// the non-singlet integrals and by $r = \frac{2}{5} \frac{16}{1 - \ln(x)}, o = 1$
-/// for the singlet sector. Note that the non-singlet kernels evolve poles only up to
+///
+/// $$p_{\text{Talbot}}(t) =  o + r \cdot ( \theta \cot(\theta) + i\theta) ~ \text{with}~ \theta = \pi(2t-1)$$
+///
+/// The default values for the parameters $r,o$ are given by
+/// $r = \frac{2}{5} \frac{16}{0.1 - \ln(x)}$ and $o = 0$ for
+/// the non-singlet integrals and $o = 1$ for the singlet sector.
+/// Note that the non-singlet kernels evolve poles only up to
 /// $N=0$ whereas the singlet kernels have poles up to $N=1$.
 pub struct TalbotPath {
     /// integration variable
@@ -26,16 +29,20 @@ pub struct TalbotPath {
 }
 
 impl TalbotPath {
-    /// Auxilary angle.
+    /// Auxiliary angle.
     fn theta(&self) -> f64 {
         PI * (2.0 * self.t - 1.0)
     }
 
     /// Constructor from parameters.
     pub fn new(t: f64, logx: f64, is_singlet: bool) -> Self {
+        // The prescription suggested by Abate for r is 0.4 * M / ( - logx)
+        // with M the number of accurate digits; Maria Ubiali suggested in her thesis M = 16.
+        // However, this seems to yield unstable results for the OME in the large x region
+        // so we add an additional regularization, which makes the path less "edgy" there.
         Self {
             t,
-            r: 0.4 * 16.0 / (1.0 - logx),
+            r: 0.4 * 16.0 / (0.1 - logx),
             o: if is_singlet { 1. } else { 0. },
         }
     }
diff --git a/crates/ekore/refs.bib b/crates/ekore/refs.bib
index 3c63b69e7..e8ef02c2b 100644
--- a/crates/ekore/refs.bib
+++ b/crates/ekore/refs.bib
@@ -50,3 +50,22 @@ @article{Abate
       journal = {International Journal for Numerical Methods in Engineering},
       doi = {10.1002/nme.995}
 }
+@phdthesis{MuselliPhD,
+    author = "Muselli, Claudio",
+    title = "{Resummations of Transverse Momentum Distributions}",
+    doi = "10.13130/muselli-claudio_phd2017-12-01",
+    school = "Università degli studi di Milano, Dipartimento di Fisica",
+    year = "2017"
+}
+@article{Vogt2004ns,
+    author = "Vogt, A.",
+    archivePrefix = "arXiv",
+    doi = "10.1016/j.cpc.2005.03.103",
+    eprint = "hep-ph/0408244",
+    journal = "Comput. Phys. Commun.",
+    pages = "65--92",
+    reportNumber = "NIKHEF-04-011",
+    title = "{Efficient evolution of unpolarized and polarized parton distributions with QCD-PEGASUS}",
+    volume = "170",
+    year = "2005"
+}
diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs
index c2d0fb0d7..dec73bf04 100644
--- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs
+++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs
@@ -1,14 +1,37 @@
 //! The unpolarized, space-like anomalous dimensions at various couplings power.
 
+use crate::constants::{PID_NSM, PID_NSP, PID_NSV};
 use crate::harmonics::cache::Cache;
 use num::complex::Complex;
 use num::Zero;
 pub mod as1;
+pub mod as2;
+pub mod as3;
 
 /// Compute the tower of the non-singlet anomalous dimensions.
-pub fn gamma_ns_qcd(order_qcd: usize, _mode: u16, c: &mut Cache, nf: u8) -> Vec<Complex<f64>> {
+pub fn gamma_ns_qcd(order_qcd: usize, mode: u16, c: &mut Cache, nf: u8) -> Vec<Complex<f64>> {
     let mut gamma_ns = vec![Complex::<f64>::zero(); order_qcd];
     gamma_ns[0] = as1::gamma_ns(c, nf);
+    // NLO and beyond
+    if order_qcd >= 2 {
+        let gamma_ns_1 = match mode {
+            PID_NSP => as2::gamma_nsp(c, nf),
+            // To fill the full valence vector in NNLO we need to add gamma_ns^1 explicitly here
+            PID_NSM | PID_NSV => as2::gamma_nsm(c, nf),
+            _ => panic!("Unkown non-singlet sector element"),
+        };
+        gamma_ns[1] = gamma_ns_1
+    }
+    // NNLO and beyond
+    if order_qcd >= 3 {
+        let gamma_ns_2 = match mode {
+            PID_NSP => as3::gamma_nsp(c, nf),
+            PID_NSM => as3::gamma_nsm(c, nf),
+            PID_NSV => as3::gamma_nsv(c, nf),
+            _ => panic!("Unkown non-singlet sector element"),
+        };
+        gamma_ns[2] = gamma_ns_2
+    }
     gamma_ns
 }
 
@@ -22,5 +45,13 @@ pub fn gamma_singlet_qcd(order_qcd: usize, c: &mut Cache, nf: u8) -> Vec<[[Compl
         order_qcd
     ];
     gamma_S[0] = as1::gamma_singlet(c, nf);
+    // NLO and beyond
+    if order_qcd >= 2 {
+        gamma_S[1] = as2::gamma_singlet(c, nf);
+    }
+    // NNLO and beyond
+    if order_qcd >= 3 {
+        gamma_S[2] = as3::gamma_singlet(c, nf);
+    }
     gamma_S
 }
diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as2.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as2.rs
new file mode 100644
index 000000000..b2b16459c
--- /dev/null
+++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as2.rs
@@ -0,0 +1,338 @@
+//! NLO QCD
+
+use ::num::complex::Complex;
+use std::f64::consts::LN_2;
+
+use crate::constants::{CA, CF, TR, ZETA2, ZETA3};
+use crate::harmonics::cache::{Cache, K};
+
+/// Compute the valence-like non-singlet anomalous dimension.
+///
+/// Implements Eq. (3.6) of [\[Moch:2004pa\]][crate::bib::Moch2004pa].
+pub fn gamma_nsm(c: &mut Cache, nf: u8) -> Complex<f64> {
+    let N = c.n;
+    let S1 = c.get(K::S1);
+    let S2 = c.get(K::S2);
+    let Sp1m = c.get(K::S1mh);
+    let Sp2m = c.get(K::S2mh);
+    let Sp3m = c.get(K::S3mh);
+    let g3n = c.get(K::G3);
+    #[rustfmt::skip]
+    let gqq1m_cfca = 16.0 * g3n
+        - (144.0 + N * (1.0 + N) * (156.0 + N * (340.0 + N * (655.0 + 51.0 * N * (2.0 + N)))))
+            / (18.0 * N.powu(3) * (1. + N).powu(3))
+        + (-14.666666666666666 + 8.0 / N - 8.0 / (1.0 + N)) * S2
+        - (4.0 * Sp2m) / (N + N.powu(2))
+        + S1 * (29.77777777777778 + 16.0 / N.powu(2) - 16.0 * S2 + 8.0 * Sp2m)
+        + 2.0 * Sp3m
+        + 10.0 * ZETA3
+        + ZETA2 * (16.0 * S1 - 16.0 * Sp1m - (16.0 * (1. + N * LN_2)) / N);
+    #[rustfmt::skip]
+    let gqq1m_cfcf = -32. * g3n
+        + (24. - N * (-32. + 3. * N * (-8. + N * (3. + N) * (3. + N.powu(2)))))
+            / (2. * N.powu(3) * (1. + N).powu(3))
+        + (12. - 8. / N + 8. / (1. + N)) * S2
+        + S1 * (-24. / N.powu(2) - 8. / (1. + N).powu(2) + 16. * S2 - 16. * Sp2m)
+        + (8. * Sp2m) / (N + N.powu(2))
+        - 4. * Sp3m
+        - 20. * ZETA3
+        + ZETA2 * (-32. * S1 + 32. * Sp1m + 32. * (1. / N + LN_2));
+    #[rustfmt::skip]
+    let gqq1m_cfnf = (-12. + N * (20. + N * (47. + 3. * N * (2. + N))))
+        / (9. * N.powu(2) * (1. + N).powu(2))
+        - (40. * S1) / 9.
+        + (8. * S2) / 3.;
+
+    CF * (CA * gqq1m_cfca + CF * gqq1m_cfcf + 2.0 * TR * (nf as f64) * gqq1m_cfnf)
+}
+
+/// Compute the singlet-like non-singlet anomalous dimension.
+///
+/// Implements Eq. (3.5) of [\[Moch:2004pa\]][crate::bib::Moch2004pa].
+pub fn gamma_nsp(c: &mut Cache, nf: u8) -> Complex<f64> {
+    let N = c.n;
+    let S1 = c.get(K::S1);
+    let S2 = c.get(K::S2);
+    let Sp1p = c.get(K::S1h);
+    let Sp2p = c.get(K::S2h);
+    let Sp3p = c.get(K::S3h);
+    let g3n = c.get(K::G3);
+    #[rustfmt::skip]
+    let gqq1p_cfca = -16. * g3n
+        + (132. - N * (340. + N * (655. + 51. * N * (2. + N))))
+            / (18. * N.powu(2) * (1. + N).powu(2))
+        + (-14.666666666666666 + 8. / N - 8. / (1. + N)) * S2
+        - (4. * Sp2p) / (N + N.powu(2))
+        + S1 * (29.77777777777778 - 16. / N.powu(2) - 16. * S2 + 8. * Sp2p)
+        + 2. * Sp3p
+        + 10. * ZETA3
+        + ZETA2 * (16. * S1 - 16. * Sp1p + 16. * (1. / N - LN_2));
+    #[rustfmt::skip]
+    let gqq1p_cfcf = 32. * g3n
+        - (8. + N * (32. + N * (40. + 3. * N * (3. + N) * (3. + N.powu(2)))))
+            / (2. * N.powu(3) * (1. + N).powu(3))
+        + (12. - 8. / N + 8. / (1. + N)) * S2
+        + S1 * (40. / N.powu(2) - 8. / (1. + N).powu(2) + 16. * S2 - 16. * Sp2p)
+        + (8. * Sp2p) / (N + N.powu(2))
+        - 4. * Sp3p
+        - 20. * ZETA3
+        + ZETA2 * (-32. * S1 + 32. * Sp1p + 32. * (-(1. / N) + LN_2));
+    #[rustfmt::skip]
+    let gqq1p_cfnf = (-12. + N * (20. + N * (47. + 3. * N * (2. + N))))
+        / (9. * N.powu(2) * (1. + N).powu(2))
+        - (40. * S1) / 9.
+        + (8. * S2) / 3.;
+
+    CF * (CA * gqq1p_cfca + CF * gqq1p_cfcf + 2.0 * TR * (nf as f64) * gqq1p_cfnf)
+}
+
+/// Compute the pure-singlet quark-quark anomalous dimension.
+///
+/// Implements Eq. (3.6) of [\[Vogt:2004mw\]][crate::bib::Vogt2004mw].
+pub fn gamma_ps(c: &mut Cache, nf: u8) -> Complex<f64> {
+    let N = c.n;
+    let gqqps1_nfcf = (-4. * (2. + N * (5. + N)) * (4. + N * (4. + N * (7. + 5. * N))))
+        / ((-1. + N) * N.powu(3) * (1. + N).powu(3) * (2. + N).powu(2));
+    2.0 * TR * (nf as f64) * CF * gqqps1_nfcf
+}
+
+/// Compute the quark-gluon singlet anomalous dimension.
+///
+/// Implements Eq. (3.7) of [\[Vogt:2004mw\]][crate::bib::Vogt2004mw].
+pub fn gamma_qg(c: &mut Cache, nf: u8) -> Complex<f64> {
+    let N = c.n;
+    let S1 = c.get(K::S1);
+    let S2 = c.get(K::S2);
+    let Sp2p = c.get(K::S2h);
+    #[rustfmt::skip]
+    let gqg1_nfca = (-4.
+        * (16.
+            + N * (64.
+                + N * (104.
+                    + N * (128. + N * (85. + N * (36. + N * (25. + N * (15. + N * (6. + N))))))))))
+        / ((-1. + N) * N.powu(3) * (1. + N).powu(3) * (2. + N).powu(3))
+        - (16. * (3. + 2. * N) * S1) / (2. + 3. * N + N.powu(2)).powu(2)
+        + (4. * (2. + N + N.powu(2)) * S1.powu(2)) / (N * (2. + 3. * N + N.powu(2)))
+        - (4. * (2. + N + N.powu(2)) * S2) / (N * (2. + 3. * N + N.powu(2)))
+        + (4. * (2. + N + N.powu(2)) * Sp2p) / (N * (2. + 3. * N + N.powu(2)));
+    #[rustfmt::skip]
+    let gqg1_nfcf = (-2. * (4. + N * (8. + N * (1. + N) * (25. + N * (26. + 5. * N * (2. + N))))))
+        / (N.powu(3) * (1. + N).powu(3) * (2. + N))
+        + (8. * S1) / N.powu(2)
+        - (4. * (2. + N + N.powu(2)) * S1.powu(2)) / (N * (2. + 3. * N + N.powu(2)))
+        + (4. * (2. + N + N.powu(2)) * S2) / (N * (2. + 3. * N + N.powu(2)));
+    2.0 * TR * (nf as f64) * (CA * gqg1_nfca + CF * gqg1_nfcf)
+}
+
+/// Compute the gluon-quark singlet anomalous dimension.
+///
+/// Implements Eq. (3.8) of [\[Vogt:2004mw\]][crate::bib::Vogt2004mw].
+pub fn gamma_gq(c: &mut Cache, nf: u8) -> Complex<f64> {
+    let N = c.n;
+    let S1 = c.get(K::S1);
+    let S2 = c.get(K::S2);
+    let Sp2p = c.get(K::S2h);
+    #[rustfmt::skip]
+    let ggq1_cfcf = (-8.
+        + 2. * N * (-12. + N * (-1. + N * (28. + N * (43. + 6. * N * (5. + 2. * N))))))
+        / ((-1. + N) * N.powu(3) * (1. + N).powu(3))
+        - (4. * (10. + N * (17. + N * (8. + 5. * N))) * S1) / ((-1. + N) * N * (1. + N).powu(2))
+        + (4. * (2. + N + N.powu(2)) * S1.powu(2)) / (N * (-1. + N.powu(2)))
+        + (4. * (2. + N + N.powu(2)) * S2) / (N * (-1. + N.powu(2)));
+    #[rustfmt::skip]
+    let ggq1_cfca = (-4.
+        * (144.
+            + N * (432.
+                + N * (-152.
+                    + N * (-1304.
+                        + N * (-1031.
+                            + N * (695. + N * (1678. + N * (1400. + N * (621. + 109. * N))))))))))
+        / (9. * N.powu(3) * (1. + N).powu(3) * (-2. + N + N.powu(2)).powu(2))
+        + (4. * (-12. + N * (-22. + 41. * N + 17. * N.powu(3))) * S1)
+            / (3. * (-1. + N).powu(2) * N.powu(2) * (1. + N))
+        + ((8. + 4. * N + 4. * N.powu(2)) * S1.powu(2)) / (N - N.powu(3))
+        + ((8. + 4. * N + 4. * N.powu(2)) * S2) / (N - N.powu(3))
+        + (4. * (2. + N + N.powu(2)) * Sp2p) / (N * (-1. + N.powu(2)));
+    #[rustfmt::skip]
+    let ggq1_cfnf = (8. * (16. + N * (27. + N * (13. + 8. * N))))
+        / (9. * (-1. + N) * N * (1. + N).powu(2))
+        - (8. * (2. + N + N.powu(2)) * S1) / (3. * N * (-1. + N.powu(2)));
+    CF * (CA * ggq1_cfca + CF * ggq1_cfcf + 2.0 * TR * (nf as f64) * ggq1_cfnf)
+}
+
+/// Compute the gluon-gluon singlet anomalous dimension.
+///
+/// Implements Eq. (3.9) of [\[Vogt:2004mw\]][crate::bib::Vogt2004mw].
+pub fn gamma_gg(c: &mut Cache, nf: u8) -> Complex<f64> {
+    let N = c.n;
+    let S1 = c.get(K::S1);
+    let Sp1p = c.get(K::S1h);
+    let Sp2p = c.get(K::S2h);
+    let Sp3p = c.get(K::S3h);
+    let g3n = c.get(K::G3);
+    #[rustfmt::skip]
+    let ggg1_caca = 16. * g3n
+        - (2.
+            * (576.
+                + N * (1488.
+                    + N * (560.
+                        + N * (-1248.
+                            + N * (-1384.
+                                + N * (1663.
+                                    + N * (4514.
+                                        + N * (4744.
+                                            + N * (3030.
+                                                + N * (1225. + 48. * N * (7. + N))))))))))))
+            / (9. * (-1. + N).powu(2) * N.powu(3) * (1. + N).powu(3) * (2. + N).powu(3))
+        + S1 * (29.77777777777778 + 16. / (-1. + N).powu(2) + 16. / (1. + N).powu(2)
+            - 16. / (2. + N).powu(2)
+            - 8. * Sp2p)
+        + (16. * (1. + N + N.powu(2)) * Sp2p) / (N * (1. + N) * (-2. + N + N.powu(2)))
+        - 2. * Sp3p
+        - 10. * ZETA3
+        + ZETA2 * (-16. * S1 + 16. * Sp1p + 16. * (-(1. / N) + LN_2));
+    #[rustfmt::skip]
+    let ggg1_canf = (8. * (6. + N * (1. + N) * (28. + N * (1. + N) * (13. + 3. * N * (1. + N)))))
+        / (9. * N.powu(2) * (1. + N).powu(2) * (-2. + N + N.powu(2)))
+        - (40. * S1) / 9.;
+    #[rustfmt::skip]
+    let ggg1_cfnf = (2.
+        * (-8.
+            + N * (-8.
+                + N * (-10. + N * (-22. + N * (-3. + N * (6. + N * (8. + N * (4. + N)))))))))
+        / (N.powu(3) * (1. + N).powu(3) * (-2. + N + N.powu(2)));
+
+    CA * CA * ggg1_caca + 2.0 * TR * (nf as f64) * (CA * ggg1_canf + CF * ggg1_cfnf)
+}
+
+/// Compute the singlet anomalous dimension matrix.
+pub fn gamma_singlet(c: &mut Cache, nf: u8) -> [[Complex<f64>; 2]; 2] {
+    let gamma_qq = gamma_nsp(c, nf) + gamma_ps(c, nf);
+    [
+        [gamma_qq, gamma_qg(c, nf)],
+        [gamma_gq(c, nf), gamma_gg(c, nf)],
+    ]
+}
+
+#[cfg(test)]
+mod tests {
+    use crate::cmplx;
+    use crate::{anomalous_dimensions::unpolarized::spacelike::as2::*, harmonics::cache::Cache};
+    use float_cmp::assert_approx_eq;
+    use num::complex::Complex;
+    use num::traits::Pow;
+    use std::f64::consts::PI;
+
+    const NF: u8 = 5;
+
+    #[test]
+    fn physical_constraints() {
+        // number conservation
+        let mut c = Cache::new(cmplx![1., 0.]);
+        assert_approx_eq!(f64, gamma_nsm(&mut c, NF).re, 0.0, epsilon = 2e-6);
+
+        // momentum conservation
+        let mut c = Cache::new(cmplx![2., 0.]);
+        let gS1 = gamma_singlet(&mut c, NF);
+
+        // gluon momentum conservation
+        assert_approx_eq!(f64, (gS1[0][1] + gS1[1][1]).re, 0.0, epsilon = 4e-5);
+        // quark momentum conservation
+        assert_approx_eq!(f64, (gS1[0][0] + gS1[1][0]).re, 0.0, epsilon = 2e-6);
+    }
+
+    #[test]
+    fn N2() {
+        // reference values are obtained from MMa
+        let mut c = Cache::new(cmplx![2., 0.]);
+
+        // ns+
+        assert_approx_eq!(
+            f64,
+            gamma_nsp(&mut c, NF).re,
+            (-112.0 * CF + 376.0 * CA - 64.0 * (NF as f64)) * CF / 27.0,
+            epsilon = 2e-6
+        );
+
+        // ns-
+        let check = ((34.0 / 27.0 * (-47.0 + 6. * PI.pow(2)) - 16.0 * ZETA3) * CF
+            + (373.0 / 9.0 - 34.0 * PI.pow(2) / 9.0 + 8.0 * ZETA3) * CA
+            - 64.0 * (NF as f64) / 27.0)
+            * CF;
+        assert_approx_eq!(f64, gamma_nsm(&mut c, NF).re, check, epsilon = 2e-6);
+
+        // singlet sector
+        let gS1 = gamma_singlet(&mut c, NF);
+        // ps
+        assert_approx_eq!(
+            f64,
+            gamma_ps(&mut c, NF).re,
+            -40.0 * CF * (NF as f64) / 27.0
+        );
+        // qg
+        assert_approx_eq!(
+            f64,
+            gS1[0][1].re,
+            (-74.0 * CF - 35.0 * CA) * (NF as f64) / 27.0
+        );
+        // gq
+        assert_approx_eq!(
+            f64,
+            gS1[1][0].re,
+            (112.0 * CF - 376.0 * CA + 104.0 * (NF as f64)) * CF / 27.0,
+            epsilon = 1e-13
+        );
+    }
+
+    #[test]
+    fn N3() {
+        let mut c = Cache::new(cmplx![3., 0.]);
+        // ns+
+        let check = ((-34487.0 / 432.0 + 86.0 * PI.pow(2) / 9.0 - 16.0 * ZETA3) * CF
+            + (459.0 / 8.0 - 43.0 * PI.pow(2) / 9.0 + 8.0 * ZETA3) * CA
+            - 415.0 * (NF as f64) / 108.0)
+            * CF;
+        assert_approx_eq!(f64, gamma_nsp(&mut c, NF).re, check, epsilon = 2e-6);
+
+        // singlet sector
+        let gS1 = gamma_singlet(&mut c, NF);
+        // ps
+        assert_approx_eq!(
+            f64,
+            gamma_ps(&mut c, NF).re,
+            -1391.0 * CF * (NF as f64) / 5400.0
+        );
+        // gq
+        assert_approx_eq!(
+            f64,
+            gS1[1][0].re,
+            (973.0 / 432.0 * CF
+                + (2801.0 / 5400.0 - 7.0 * PI.pow(2) / 9.0) * CA
+                + 61.0 / 54.0 * (NF as f64))
+                * CF
+        );
+        // gg
+        assert_approx_eq!(
+            f64,
+            gS1[1][1].re,
+            (-79909.0 / 3375.0 + 194.0 * PI.pow(2) / 45.0 - 8.0 * ZETA3) * CA.pow(2)
+                - 967.0 / 270.0 * CA * (NF as f64)
+                + 541.0 / 216.0 * CF * (NF as f64),
+            epsilon = 3e-5
+        );
+    }
+
+    #[test]
+    fn N4() {
+        let mut c = Cache::new(cmplx![4., 0.]);
+        // singlet sector
+        let gS1 = gamma_singlet(&mut c, NF);
+        // qg
+        assert_approx_eq!(
+            f64,
+            gS1[0][1].re,
+            (-56317.0 / 18000.0 * CF + 16387.0 / 9000.0 * CA) * (NF as f64),
+            epsilon = 1e-14
+        )
+    }
+}
diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as3.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as3.rs
new file mode 100644
index 000000000..b43608e01
--- /dev/null
+++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as3.rs
@@ -0,0 +1,456 @@
+//! NNLO QCD
+
+use ::num::complex::Complex;
+use num::traits::Pow;
+
+use crate::cmplx;
+use crate::constants::{ZETA2, ZETA3};
+use crate::harmonics::cache::{Cache, K};
+
+/// Compute the valence-like non-singlet anomalous dimension.
+///
+/// Implements Eq. (3.8) of [\[Moch:2004pa\]][crate::bib::Moch2004pa].
+pub fn gamma_nsm(c: &mut Cache, nf: u8) -> Complex<f64> {
+    let N = c.n;
+    let S1 = c.get(K::S1);
+    let S2 = c.get(K::S2);
+    let S3 = c.get(K::S3);
+    let E1 = S1 / N.powu(2) + (S2 - ZETA2) / N;
+    let E2 = 2.0 * (-S1 / N.powu(3) + (ZETA2 - S2) / N.powu(2) - (S3 - ZETA3) / N);
+
+    #[rustfmt::skip]
+    let pm2 =
+        -1174.898 * (S1 - 1.0 / N)
+        + 1295.470
+        - 714.1 * S1 / N
+        - 433.2 / (N + 3.0)
+        + 297.0 / (N + 2.0)
+        - 3505.0 / (N + 1.0)
+        + 1860.2 / N
+        - 1465.2 / N.powu(2)
+        + 399.2 * 2.0 / N.powu(3)
+        - 320.0 / 9.0 * 6.0 / N.powu(4)
+        + 116.0 / 81.0 * 24.0 / N.powu(5)
+        + 684.0 * E1
+        + 251.2 * E2;
+
+    #[rustfmt::skip]
+    let pm2_nf =
+        183.187 * (S1 - 1.0 / N)
+        - 173.933
+        + 5120./ 81.0 * S1 / N
+        + 34.76 / (N + 3.0)
+        + 77.89 / (N + 2.0)
+        + 406.5 / (N + 1.0)
+        - 216.62 / N
+        + 172.69 / N.powu(2)
+        - 3216.0 / 81.0 * 2.0 / N.powu(3)
+        + 256.0 / 81.0 * 6.0 / N.powu(4)
+        - 65.43 * E1
+        + 1.136 * 6.0 / (N + 1.0).powu(4);
+
+    #[rustfmt::skip]
+    let pf2_nfnf =
+        -(
+            17.0 / 72.0
+            - 2.0 / 27.0 * S1
+            - 10.0 / 27.0 * S2
+            + 2.0 / 9.0 * S3
+            - (12.0 * N.powu(4) + 2.0 * N.powu(3) - 12.0 * N.powu(2) - 2.0 * N + 3.0)
+            / (27.0 * N.powu(3) * (N + 1.0).powu(3))
+        )* 32.0 / 3.0;
+
+    let result = pm2 + (nf as f64) * pm2_nf + (nf as f64).pow(2) * pf2_nfnf;
+    -1. * result
+}
+
+/// Compute the singlet-like non-singlet anomalous dimension.
+///
+/// Implements Eq. (3.7) of [\[Moch:2004pa\]][crate::bib::Moch2004pa].
+pub fn gamma_nsp(c: &mut Cache, nf: u8) -> Complex<f64> {
+    let N = c.n;
+    let S1 = c.get(K::S1);
+    let S2 = c.get(K::S2);
+    let S3 = c.get(K::S3);
+    let E1 = S1 / N.powu(2) + (S2 - ZETA2) / N;
+    let E2 = 2.0 * (-S1 / N.powu(3) + (ZETA2 - S2) / N.powu(2) - (S3 - ZETA3) / N);
+
+    #[rustfmt::skip]
+    let pp2 =
+        -1174.898 * (S1 - 1.0 / N)
+        + 1295.384
+        - 714.1 * S1 / N
+        - 522.1 / (N + 3.0)
+        + 243.6 / (N + 2.0)
+        - 3135.0 / (N + 1.0)
+        + 1641.1 / N
+        - 1258.0 / N.powu(2)
+        + 294.9 * 2.0 / N.powu(3)
+        - 800. / 27.0 * 6.0 / N.powu(4)
+        + 128. / 81.0 * 24.0 / N.powu(5)
+        + 563.9 * E1
+        + 256.8 * E2;
+
+    #[rustfmt::skip]
+    let pp2_nf =
+        183.187 * (S1 - 1.0 / N)
+        - 173.924
+        + 5120. / 81.0 * S1 / N
+        + 44.79 / (N + 3.0)
+        + 72.94 / (N + 2.0)
+        + 381.1 / (N + 1.0)
+        - 197.0 / N
+        + 152.6 / N.powu(2)
+        - 2608.0 / 81.0 * 2.0 / N.powu(3)
+        + 192.0 / 81.0 * 6.0 / N.powu(4)
+        - 56.66 * E1
+        + 1.497 * 6.0 / (N + 1.0).powu(4);
+
+    #[rustfmt::skip]
+    let pf2_nfnf =
+        -(
+            17.0 / 72.0
+            - 2.0 / 27.0 * S1
+            - 10.0 / 27.0 * S2
+            + 2.0 / 9.0 * S3
+            - (12.0 * N.powu(4) + 2.0 * N.powu(3) - 12.0 * N.powu(2) - 2.0 * N + 3.0)
+            / (27.0 * N.powu(3) * (N + 1.0).powu(3))
+        )* 32.0/ 3.0;
+
+    let result = pp2 + (nf as f64) * pp2_nf + (nf as f64).pow(2) * pf2_nfnf;
+    -1. * result
+}
+
+/// Compute the valence non-singlet anomalous dimension.
+///
+/// Implements Eq. (3.9) of [\[Moch:2004pa\]][crate::bib::Moch2004pa].
+pub fn gamma_nsv(c: &mut Cache, nf: u8) -> Complex<f64> {
+    let N = c.n;
+    let S1 = c.get(K::S1);
+    let S2 = c.get(K::S2);
+    let S3 = c.get(K::S3);
+    let E1 = S1 / N.powu(2) + (S2 - ZETA2) / N;
+    let E2 = 2.0 * (-S1 / N.powu(3) + (ZETA2 - S2) / N.powu(2) - (S3 - ZETA3) / N);
+    let B11 = -(S1 + 1.0 / (N + 1.0)) / (N + 1.0);
+    let B12 = -(S1 + 1.0 / (N + 1.0) + 1.0 / (N + 2.0)) / (N + 2.0);
+
+    let B1M = if N.im.abs() < 1.0e-5 && (N - 1.0).re.abs() < 1.0e-5 {
+        cmplx![-ZETA2, 0.]
+    } else {
+        -(S1 - 1.0 / N) / (N - 1.0)
+    };
+
+    #[rustfmt::skip]
+    let ps2 = -(
+        -163.9 * (B1M + S1 / N)
+        - 7.208 * (B11 - B12)
+        + 4.82 * (1.0 / (N + 3.0) - 1.0 / (N + 4.0))
+        - 43.12 * (1.0 / (N + 2.0) - 1.0 / (N + 3.0))
+        + 44.51 * (1.0 / (N + 1.0) - 1.0 / (N + 2.0))
+        + 151.49 * (1.0 / N - 1.0 / (N + 1.0))
+        - 178.04 / N.powu(2)
+        + 6.892 * 2.0 / N.powu(3)
+        - 40.0 / 27.0 * (-2.0 * 6.0 / N.powu(4) - 24.0 / N.powu(5))
+        - 173.1 * E1
+        + 46.18 * E2
+    );
+
+    let result = gamma_nsm(c, nf) + (nf as f64) * ps2;
+    result
+}
+
+/// Compute the pure-singlet quark-quark anomalous dimension.
+///
+/// Implements Eq. (3.10) of [\[Vogt:2004mw\]][crate::bib::Vogt2004mw].
+pub fn gamma_ps(c: &mut Cache, nf: u8) -> Complex<f64> {
+    let N = c.n;
+    let S1 = c.get(K::S1);
+    let S2 = c.get(K::S2);
+    let S3 = c.get(K::S3);
+
+    let E1 = S1 / N.powu(2) + (S2 - ZETA2) / N;
+    let E11 = (S1 + 1.0 / (N + 1.0)) / (N + 1.0).powu(2)
+        + (S2 + 1.0 / (N + 1.0).powu(2) - ZETA2) / (N + 1.0);
+    let B21 = ((S1 + 1.0 / (N + 1.0)).powu(2) + S2 + 1.0 / (N + 1.0).powu(2)) / (N + 1.0);
+    let B3 = -(S1.powu(3) + 3.0 * S1 * S2 + 2.0 * S3) / N;
+
+    #[rustfmt::skip]
+    let B31 = -(
+        (S1 + 1.0 / (N + 1.0)).powu(3)
+        + 3.0 * (S1 + 1.0 / (N + 1.0)) * (S2 + 1.0 / (N + 1.0).powu(2))
+        + 2.0 * (S3 + 1.0 / (N + 1.0).powu(3))
+    ) / (N + 1.0);
+
+    #[rustfmt::skip]
+    let ps1 =
+        -3584.0 / 27.0 * (-1.0 / (N - 1.0).powu(2) + 1.0 / N.powu(2))
+        - 506.0 * (1.0 / (N - 1.0) - 1.0 / N)
+        + 160.0 / 27.0 * (24.0 / N.powu(5) - 24.0 / (N + 1.0).powu(5))
+        - 400.0 / 9.0 * (-6.0 / N.powu(4) + 6.0 / (N + 1.0).powu(4))
+        + 131.4 * (2.0 / N.powu(3) - 2.0 / (N + 1.0).powu(3))
+        - 661.6 * (-1.0 / N.powu(2) + 1.0 / (N + 1.0).powu(2))
+        - 5.926 * (B3 - B31)
+        - 9.751 * ((S1.powu(2) + S2) / N - B21)
+        - 72.11 * (-S1 / N + (S1 + 1.0 / (N + 1.0)) / (N + 1.0))
+        + 177.4 * (1.0 / N - 1.0 / (N + 1.0))
+        + 392.9 * (1.0 / (N + 1.0) - 1.0 / (N + 2.0))
+        - 101.4 * (1.0 / (N + 2.0) - 1.0 / (N + 3.0))
+        - 57.04 * (E1 - E11);
+
+    #[rustfmt::skip]
+    let ps2 =
+        256.0 / 81.0 * (1.0 / (N - 1.0) - 1.0 / N)
+        + 32.0 / 27.0 * (-6.0 / N.powu(4) + 6.0 / (N + 1.0).powu(4))
+        + 17.89 * (2.0 / N.powu(3) - 2.0 / (N + 1.0).powu(3))
+        + 61.75 * (-1.0 / N.powu(2) + 1.0 / (N + 1.0).powu(2))
+        + 1.778 * ((S1.powu(2) + S2) / N - B21)
+        + 5.944 * (-S1 / N + (S1 + 1.0 / (N + 1.0)) / (N + 1.0))
+        + 100.1 * (1.0 / N - 1.0 / (N + 1.0))
+        - 125.2 * (1.0 / (N + 1.0) - 1.0 / (N + 2.0))
+        + 49.26 * (1.0 / (N + 2.0) - 1.0 / (N + 3.0))
+        - 12.59 * (1.0 / (N + 3.0) - 1.0 / (N + 4.0))
+        - 1.889 * (E1 - E11);
+
+    let result = (nf as f64) * ps1 + (nf as f64).pow(2) * ps2;
+    -1.0 * result
+}
+
+/// Compute the quark-gluon singlet anomalous dimension.
+///
+/// Implements Eq. (3.11) of [\[Vogt:2004mw\]][crate::bib::Vogt2004mw].
+pub fn gamma_qg(c: &mut Cache, nf: u8) -> Complex<f64> {
+    let N = c.n;
+    let S1 = c.get(K::S1);
+    let S2 = c.get(K::S2);
+    let S3 = c.get(K::S3);
+    let S4 = c.get(K::S4);
+
+    let E1 = S1 / N.powu(2) + (S2 - ZETA2) / N;
+    let E2 = 2.0 * (-S1 / N.powu(3) + (ZETA2 - S2) / N.powu(2) - (S3 - ZETA3) / N);
+    let B3 = -(S1.powu(3) + 3.0 * S1 * S2 + 2.0 * S3) / N;
+    let B4 = (S1.powu(4) + 6.0 * S1.powu(2) * S2 + 8.0 * S1 * S3 + 3.0 * S2.powu(2) + 6.0 * S4) / N;
+
+    #[rustfmt::skip]
+    let qg1 =
+        896.0 / 3.0 / (N - 1.0).powu(2)
+        - 1268.3 / (N - 1.0)
+        + 536.0 / 27.0 * 24.0 / N.powu(5)
+        + 44.0 / 3.0 * 6.0 / N.powu(4)
+        + 881.5 * 2.0 / N.powu(3)
+        - 424.9 / N.powu(2)
+        + 100.0 / 27.0 * B4
+        - 70.0 / 9.0 * B3
+        - 120.5 * (S1.powu(2) + S2) / N
+        - 104.42 * S1 / N
+        + 2522.0 / N
+        - 3316.0 / (N + 1.0)
+        + 2126.0 / (N + 2.0)
+        + 1823.0 * E1
+        - 25.22 * E2
+        + 252.5 * 6.0 / (N + 1.0).powu(4);
+
+    #[rustfmt::skip]
+    let qg2 =
+        1112.0 / 243.0 / (N - 1.0)
+        - 16.0 / 9.0 * 24.0 / N.powu(5)
+        + 376.0 / 27.0 * 6.0 / N.powu(4)
+        - 90.8 * 2.0 / N.powu(3)
+        + 254.0 / N.powu(2)
+        + 20.0 / 27.0 * B3
+        + 200.0 / 27.0 * (S1.powu(2) + S2) / N
+        + 5.496 * S1 / N
+        - 252.0 / N
+        + 158.0 / (N + 1.0)
+        + 145.4 / (N + 2.0)
+        - 139.28 / (N + 3.0)
+        - 53.09 * E1
+        - 80.616 * E2
+        - 98.07 * 2.0 / (N + 1.0).powu(3)
+        - 11.70 * 6.0 / (N + 1.0).powu(4);
+
+    let result = (nf as f64) * qg1 + (nf as f64).pow(2) * qg2;
+    -1.0 * result
+}
+
+/// Compute the gluon-quark singlet anomalous dimension.
+///
+/// Implements Eq. (3.12) of [\[Vogt:2004mw\]][crate::bib::Vogt2004mw].
+pub fn gamma_gq(c: &mut Cache, nf: u8) -> Complex<f64> {
+    let N = c.n;
+    let S1 = c.get(K::S1);
+    let S2 = c.get(K::S2);
+    let S3 = c.get(K::S3);
+    let S4 = c.get(K::S4);
+
+    let E1 = S1 / N.powu(2) + (S2 - ZETA2) / N;
+    let E2 = 2.0 * (-S1 / N.powu(3) + (ZETA2 - S2) / N.powu(2) - (S3 - ZETA3) / N);
+    let B21 = ((S1 + 1.0 / (N + 1.0)).powu(2) + S2 + 1.0 / (N + 1.0).powu(2)) / (N + 1.0);
+    let B3 = -(S1.powu(3) + 3.0 * S1 * S2 + 2.0 * S3) / N;
+    let B4 = (S1.powu(4) + 6.0 * S1.powu(2) * S2 + 8.0 * S1 * S3 + 3.0 * S2.powu(2) + 6.0 * S4) / N;
+
+    #[rustfmt::skip]
+    let gq0 =
+        -1189.3 * 1.0 / (N - 1.0).powu(2)
+        + 6163.1 / (N - 1.0)
+        - 4288.0 / 81.0 * 24.0 / N.powu(5)
+        - 1568.0 / 9.0 * 6.0 / N.powu(4)
+        - 1794.0 * 2.0 / N.powu(3)
+        - 4033.0 * 1.0 / N.powu(2)
+        + 400.0 / 81.0 * B4
+        + 2200.0 / 27.0 * B3
+        + 606.3 * (S1.powu(2) + S2) / N
+        - 2193.0 * S1 / N
+        - 4307.0 / N
+        + 489.3 / (N + 1.0)
+        + 1452.0 / (N + 2.0)
+        + 146.0 / (N + 3.0)
+        - 447.3 * E2
+        - 972.9 * 2.0 / (N + 1.0).powu(3);
+
+    #[rustfmt::skip]
+    let gq1 =
+        -71.082 / (N - 1.0).powu(2)
+        - 46.41 / (N - 1.0)
+        + 128.0 / 27.0 * 24.0 / N.powu(5)
+        - 704. / 81.0 * 6.0 / N.powu(4)
+        + 20.39 * 2.0 / N.powu(3)
+        - 174.8 * 1.0 / N.powu(2)
+        - 400.0 / 81.0 * B3
+        - 68.069 * (S1.powu(2) + S2) / N
+        + 296.7 * S1 / N
+        - 183.8 / N
+        + 33.35 / (N + 1.0)
+        - 277.9 / (N + 2.0)
+        + 108.6 * 2.0 / (N + 1.0).powu(3)
+        - 49.68 * E1;
+
+    #[rustfmt::skip]
+    let gq2 = (
+        64.0 * (-1.0 / (N - 1.0) + 1.0 / N + 2.0 / (N + 1.0))
+        + 320.0
+        * (
+            -(S1 - 1.0 / N) / (N - 1.0)
+            + S1 / N
+            - 0.8 * (S1 + 1.0 / (N + 1.0)) / (N + 1.0)
+        )
+        + 96.0
+        * (
+            ((S1 - 1.0 / N).powu(2) + S2 - 1.0 / N.powu(2)) / (N - 1.0)
+            - (S1.powu(2) + S2) / N
+            + 0.5 * B21
+        )
+    ) / 27.0;
+
+    let result = gq0 + (nf as f64) * gq1 + (nf as f64).pow(2) * gq2;
+    -1.0 * result
+}
+
+/// Compute the gluon-quark singlet anomalous dimension.
+///
+/// Implements Eq. (3.13) of [\[Vogt:2004mw\]][crate::bib::Vogt2004mw].
+pub fn gamma_gg(c: &mut Cache, nf: u8) -> Complex<f64> {
+    let N = c.n;
+    let S1 = c.get(K::S1);
+    let S2 = c.get(K::S2);
+    let S3 = c.get(K::S3);
+
+    let E1 = S1 / N.powu(2) + (S2 - ZETA2) / N;
+    let E11 = (S1 + 1.0 / (N + 1.0)) / (N + 1.0).powu(2)
+        + (S2 + 1.0 / (N + 1.0).powu(2) - ZETA2) / (N + 1.0);
+    let E2 = 2.0 * (-S1 / N.powu(3) + (ZETA2 - S2) / N.powu(2) - (S3 - ZETA3) / N);
+
+    #[rustfmt::skip]
+    let gg0 =
+        -2675.8 / (N - 1.0).powu(2)
+        + 14214.0 / (N - 1.0)
+        - 144.0 * 24.0 / N.powu(5)
+        - 72.0 * 6.0 / N.powu(4)
+        - 7471.0 * 2.0 / N.powu(3)
+        - 274.4 / N.powu(2)
+        - 20852.0 / N
+        + 3968.0 / (N + 1.0)
+        - 3363.0 / (N + 2.0)
+        + 4848.0 / (N + 3.0)
+        + 7305.0 * E1
+        + 8757.0 * E2
+        - 3589.0 * S1 / N
+        + 4425.894
+        - 2643.521 * (S1 - 1.0 / N);
+
+    #[rustfmt::skip]
+    let gg1 =
+        -157.27 / (N - 1.0).powu(2)
+        + 182.96 / (N - 1.0)
+        + 512.0 / 27.0 * 24.0 / N.powu(5)
+        - 832.0 / 9.0 * 6.0 / N.powu(4)
+        + 491.3 * 2.0 / N.powu(3)
+        - 1541.0 / N.powu(2)
+        - 350.2 / N
+        + 755.7 / (N + 1.0)
+        - 713.8 / (N + 2.0)
+        + 559.3 / (N + 3.0)
+        + 26.15 * E1
+        - 808.7 * E2
+        + 320.0 * S1 / N
+        - 528.723
+        + 412.172 * (S1 - 1.0 / N);
+
+    #[rustfmt::skip]
+    let gg2 =
+        -680.0 / 243.0 / (N - 1.0)
+        + 32.0 / 27.0 * 6.0 / N.powu(4)
+        + 9.680 * 2.0 / N.powu(3)
+        + 3.422 / N.powu(2)
+        - 13.878 / N
+        + 153.4 / (N + 1.0)
+        - 187.7 / (N + 2.0)
+        + 52.75 / (N + 3.0)
+        - 115.6 * E1
+        + 85.25 * E11
+        - 63.23 * E2
+        + 6.4630
+        + 16.0 / 9.0 * (S1 - 1.0 / N);
+
+    let result = gg0 + (nf as f64) * gg1 + (nf as f64).pow(2) * gg2;
+    -1.0 * result
+}
+
+/// Compute the singlet anomalous dimension matrix.
+pub fn gamma_singlet(c: &mut Cache, nf: u8) -> [[Complex<f64>; 2]; 2] {
+    let gamma_qq = gamma_nsp(c, nf) + gamma_ps(c, nf);
+    [
+        [gamma_qq, gamma_qg(c, nf)],
+        [gamma_gq(c, nf), gamma_gg(c, nf)],
+    ]
+}
+
+#[cfg(test)]
+mod tests {
+    use crate::cmplx;
+    use crate::{anomalous_dimensions::unpolarized::spacelike::as3::*, harmonics::cache::Cache};
+    use float_cmp::assert_approx_eq;
+    use num::complex::Complex;
+
+    const NF: u8 = 5;
+
+    #[test]
+    fn physical_constraints() {
+        // number conservation
+        let mut c = Cache::new(cmplx![1., 0.]);
+        assert_approx_eq!(f64, gamma_nsv(&mut c, NF).re, -0.000960586, epsilon = 3e-7);
+        assert_approx_eq!(f64, gamma_nsm(&mut c, NF).re, 0.000594225, epsilon = 6e-7);
+
+        let mut c = Cache::new(cmplx![2., 0.]);
+        let gS2 = gamma_singlet(&mut c, NF);
+        // gluon momentum conservation
+        assert_approx_eq!(f64, (gS2[0][1] + gS2[1][1]).re, -0.00388726, epsilon = 2e-6);
+        // quark momentum conservation
+        assert_approx_eq!(f64, (gS2[0][0] + gS2[1][0]).re, 0.00169375, epsilon = 2e-6);
+    }
+
+    #[test]
+    fn N2() {
+        let mut c = Cache::new(cmplx![2., 0.]);
+        assert_approx_eq!(f64, gamma_nsv(&mut c, NF).re, 188.325593, epsilon = 3e-7);
+    }
+}
diff --git a/crates/ekore/src/bib.rs b/crates/ekore/src/bib.rs
index c62f16d02..167294b41 100644
--- a/crates/ekore/src/bib.rs
+++ b/crates/ekore/src/bib.rs
@@ -1,4 +1,4 @@
-//! List of References (autogenerated on 2023-11-29T18:07:41.518255).
+//! List of References (autogenerated on 2024-04-25T18:22:43.797877).
 
 /// The Three loop splitting functions in QCD: The Nonsinglet case
 ///
@@ -43,3 +43,25 @@ pub fn KOLBIG1972221() {}
 ///
 /// DOI: [10.1002/nme.995](https:dx.doi.org/10.1002/nme.995)
 pub fn Abate() {}
+
+/// Resummations of Transverse Momentum Distributions
+///
+/// Muselli, Claudio
+///
+/// Published as PhD thesis at Università degli studi di Milano, Dipartimento di Fisica (2017)
+///
+///
+///
+/// DOI: [10.13130/muselli-claudio_phd2017-12-01](https:dx.doi.org/10.13130/muselli-claudio_phd2017-12-01)
+pub fn MuselliPhD() {}
+
+/// Efficient evolution of unpolarized and polarized parton distributions with QCD-PEGASUS
+///
+/// Vogt, A.
+///
+/// Published in: Comput. Phys. Commun. 170 (2005), 65--92
+///
+/// e-Print: [hep-ph/0408244](https://arxiv.org/abs/hep-ph/0408244)
+///
+/// DOI: [10.1016/j.cpc.2005.03.103](https:dx.doi.org/10.1016/j.cpc.2005.03.103)
+pub fn Vogt2004ns() {}
diff --git a/crates/ekore/src/constants.rs b/crates/ekore/src/constants.rs
index bfe44fe9c..c824bdf5a 100644
--- a/crates/ekore/src/constants.rs
+++ b/crates/ekore/src/constants.rs
@@ -19,3 +19,25 @@ pub const CA: f64 = NC as f64;
 ///
 /// Defaults to $C_F = \frac{N_C^2-1}{2N_C} = 4/3$.
 pub const CF: f64 = ((NC * NC - 1) as f64) / ((2 * NC) as f64);
+
+/// Riemann zeta function at z = 2.
+///
+/// $\zeta(2) = \pi^2 / 6$.
+pub const ZETA2: f64 = 1.6449340668482264;
+
+/// Riemann zeta function at z = 3.
+pub const ZETA3: f64 = 1.2020569031595942;
+
+/// Riemann zeta function at z = 4.
+///
+/// $\zeta(4) = \pi^4 / 90$.
+pub const ZETA4: f64 = 1.082323233711138;
+
+/// singlet-like non-singlet PID
+pub const PID_NSP: u16 = 10101;
+
+/// valence-like non-singlet PID
+pub const PID_NSM: u16 = 10201;
+
+/// non-singlet all-valence PID
+pub const PID_NSV: u16 = 10200;
diff --git a/crates/ekore/src/harmonics.rs b/crates/ekore/src/harmonics.rs
index b829060fc..a8a3b3241 100644
--- a/crates/ekore/src/harmonics.rs
+++ b/crates/ekore/src/harmonics.rs
@@ -1,5 +1,9 @@
 //! Tools to compute harmonic sums and related special functions.
 
 pub mod cache;
+pub mod g_functions;
 pub mod polygamma;
-mod w1;
+pub mod w1;
+pub mod w2;
+pub mod w3;
+pub mod w4;
diff --git a/crates/ekore/src/harmonics/cache.rs b/crates/ekore/src/harmonics/cache.rs
index b026e4f75..0fed95258 100644
--- a/crates/ekore/src/harmonics/cache.rs
+++ b/crates/ekore/src/harmonics/cache.rs
@@ -1,15 +1,35 @@
 //! Cache harmonic sums for given Mellin N.
 
 use hashbrown::HashMap;
-use num::complex::Complex;
+use num::{complex::Complex, Zero};
 
-use crate::harmonics::w1;
+use crate::harmonics::{g_functions, w1, w2, w3, w4};
 
 /// List of available elements.
 #[derive(Debug, PartialEq, Eq, Hash)]
 pub enum K {
     /// $S_1(N)$
     S1,
+    /// $S_2(N)$
+    S2,
+    /// $S_3(N)$
+    S3,
+    /// $S_4(N)$
+    S4,
+    /// $S_1(N/2)$
+    S1h,
+    /// $S_2(N/2)$
+    S2h,
+    /// $S_3(N/2)$
+    S3h,
+    /// $S_1((N-1)/2)$
+    S1mh,
+    /// $S_2((N-1)/2)$
+    S2mh,
+    /// $S_3((N-1)/2)$
+    S3mh,
+    /// $g_3(N)$
+    G3,
 }
 
 /// Hold all cached values.
@@ -39,9 +59,62 @@ impl Cache {
         // compute new
         let val = match k {
             K::S1 => w1::S1(self.n),
+            K::S2 => w2::S2(self.n),
+            K::S3 => w3::S3(self.n),
+            K::S4 => w4::S4(self.n),
+            K::S1h => w1::S1(self.n / 2.),
+            K::S2h => w2::S2(self.n / 2.),
+            K::S3h => w3::S3(self.n / 2.),
+            K::S1mh => w1::S1((self.n - 1.) / 2.),
+            K::S2mh => w2::S2((self.n - 1.) / 2.),
+            K::S3mh => w3::S3((self.n - 1.) / 2.),
+            K::G3 => g_functions::g3(self.n, self.get(K::S1)),
         };
         // insert
         self.m.insert(k, val);
         val
     }
 }
+
+/// Recursive computation of harmonic sums.
+///
+/// Compute the harmonic sum $S_{w}(N+k)$ stating from the value $S_{w}(N)$ via the recurrence relations.
+pub fn recursive_harmonic_sum(
+    base_value: Complex<f64>,
+    n: Complex<f64>,
+    iterations: usize,
+    weight: u32,
+) -> Complex<f64> {
+    let mut fact = Complex::zero();
+    for i in 1..iterations + 1 {
+        fact = fact + (1.0 / (n + (i as f64))).powu(weight);
+    }
+    base_value + fact
+}
+
+#[cfg(test)]
+mod tests {
+    use crate::harmonics::cache::recursive_harmonic_sum;
+    use crate::harmonics::{w1, w2, w3, w4};
+    use crate::{assert_approx_eq_cmplx, cmplx};
+    use num::complex::Complex;
+
+    #[test]
+    fn test_recursive_harmonic_sum() {
+        const SX: [fn(Complex<f64>) -> Complex<f64>; 4] = [w1::S1, w2::S2, w3::S3, w4::S4];
+        const NS: [Complex<f64>; 2] = [cmplx![1.0, 0.0], cmplx![2.34, 3.45]];
+        const ITERS: [usize; 2] = [1, 2];
+        for sit in SX.iter().enumerate() {
+            for nit in NS.iter().enumerate() {
+                let n = *nit.1;
+                for iit in ITERS.iter().enumerate() {
+                    let iterations = *iit.1;
+                    let s_base = sit.1(n);
+                    let s_test = sit.1(n + (iterations as f64));
+                    let s_ref = recursive_harmonic_sum(s_base, n, iterations, 1 + (sit.0 as u32));
+                    assert_approx_eq_cmplx!(f64, s_test, s_ref);
+                }
+            }
+        }
+    }
+}
diff --git a/crates/ekore/src/harmonics/g_functions.rs b/crates/ekore/src/harmonics/g_functions.rs
new file mode 100644
index 000000000..ff0f8789f
--- /dev/null
+++ b/crates/ekore/src/harmonics/g_functions.rs
@@ -0,0 +1,50 @@
+//! Auxilary functions for harmonics sums of weight = 3,4.
+
+use crate::constants::ZETA2;
+use crate::harmonics::cache::recursive_harmonic_sum as s;
+use num::{complex::Complex, Zero};
+
+/// Compute the Mellin transform of $\text{Li}_2(x)/(1+x)$.
+///
+/// This function appears in the analytic continuation of the harmonic sum
+/// $S_{-2,1}(N)$ which in turn appears in the NLO anomalous dimension.
+///
+/// We use the name from [\[MuselliPhD\]](crate::bib::MuselliPhD), but not his implementation - rather we use the
+/// Pegasus [\[Vogt:2004ns\]](crate::bib::Vogt2004ns) implementation.
+pub fn g3(N: Complex<f64>, S1: Complex<f64>) -> Complex<f64> {
+    const CS: [f64; 7] = [
+        1.0000e0, -0.9992e0, 0.9851e0, -0.9005e0, 0.6621e0, -0.3174e0, 0.0699e0,
+    ];
+    let mut g3 = Complex::zero();
+    for cit in CS.iter().enumerate() {
+        let Nj = N + (cit.0 as f64);
+        g3 = g3 + (*cit.1) * (ZETA2 - s(S1, N, cit.0, 1) / Nj) / Nj;
+    }
+    g3
+}
+
+#[cfg(test)]
+mod tests {
+    use crate::harmonics::g_functions::g3;
+    use crate::harmonics::w1;
+    use crate::{assert_approx_eq_cmplx, cmplx};
+    use num::complex::Complex;
+
+    #[test]
+    fn test_mellin_g3() {
+        const NS: [Complex<f64>; 3] = [cmplx![1.0, 0.0], cmplx![2.0, 0.0], cmplx![1.0, 1.0]];
+        // NIntegrate[x^({1, 2, 1 + I} - 1) PolyLog[2, x]/(1 + x), {x, 0, 1}]
+        const REFVALS: [Complex<f64>; 3] = [
+            cmplx![0.3888958462, 0.],
+            cmplx![0.2560382207, 0.],
+            cmplx![0.3049381491, -0.1589060625],
+        ];
+        for it in NS.iter().enumerate() {
+            let n = *it.1;
+            let s1 = w1::S1(n);
+            let refval = REFVALS[it.0];
+            let g3 = g3(n, s1);
+            assert_approx_eq_cmplx![f64, g3, refval, epsilon = 1e-6];
+        }
+    }
+}
diff --git a/crates/ekore/src/harmonics/polygamma.rs b/crates/ekore/src/harmonics/polygamma.rs
index f5236d8e9..7415db88e 100644
--- a/crates/ekore/src/harmonics/polygamma.rs
+++ b/crates/ekore/src/harmonics/polygamma.rs
@@ -6,7 +6,8 @@ use std::f64::consts::PI;
 #[allow(clippy::excessive_precision, clippy::assign_op_pattern)]
 /// Compute the polygamma functions $\psi_k(z)$.
 ///
-/// Reimplementation of ``WPSIPG`` (C317) in [CERNlib](http://cernlib.web.cern.ch/cernlib/) given by [[KOLBIG1972221]][crate::bib::KOLBIG1972221].
+/// Reimplementation of ``WPSIPG`` (C317) in [CERNlib](http://cernlib.web.cern.ch/cernlib/)
+/// given by [\[KOLBIG1972221\]](crate::bib::KOLBIG1972221).
 ///
 /// TODO: introduce back errors
 pub fn cern_polygamma(Z: Complex<f64>, K: usize) -> Complex<f64> {
@@ -83,9 +84,7 @@ pub fn cern_polygamma(Z: Complex<f64>, K: usize) -> Complex<f64> {
     }
     let mut R = 1. / V.powu(2);
     let mut P = R * C[K][6 - 1];
-    for i in (1..=5).rev()
-    // (int i = 5; i>1-1; i--)
-    {
+    for i in (1..=5).rev() {
         P = R * (C[K][i - 1] + P);
     }
     H = (SGN[K] as f64)
@@ -121,9 +120,8 @@ pub fn cern_polygamma(Z: Complex<f64>, K: usize) -> Complex<f64> {
 
 #[cfg(test)]
 mod tests {
-    use crate::cmplx;
     use crate::harmonics::polygamma::cern_polygamma;
-    use float_cmp::assert_approx_eq;
+    use crate::{assert_approx_eq_cmplx, cmplx};
     use num::complex::Complex;
 
     #[test]
@@ -202,8 +200,7 @@ mod tests {
             for zit in ZS.iter().enumerate() {
                 let fref = FORTRAN_REF[kit.0][zit.0];
                 let me = cern_polygamma(*zit.1, *kit.1);
-                assert_approx_eq!(f64, me.re, fref.re, ulps = 32);
-                assert_approx_eq!(f64, me.im, fref.im, ulps = 32);
+                assert_approx_eq_cmplx!(f64, me, fref, ulps = 32);
             }
         }
     }
diff --git a/crates/ekore/src/harmonics/w1.rs b/crates/ekore/src/harmonics/w1.rs
index 7de50bf0b..a1aa893d9 100644
--- a/crates/ekore/src/harmonics/w1.rs
+++ b/crates/ekore/src/harmonics/w1.rs
@@ -1,3 +1,4 @@
+//! Harmonic sums of weight 1.
 use num::complex::Complex;
 
 use crate::harmonics::polygamma::cern_polygamma;
diff --git a/crates/ekore/src/harmonics/w2.rs b/crates/ekore/src/harmonics/w2.rs
new file mode 100644
index 000000000..d55a32679
--- /dev/null
+++ b/crates/ekore/src/harmonics/w2.rs
@@ -0,0 +1,13 @@
+//! Harmonic sums of weight 2.
+use num::complex::Complex;
+
+use crate::constants::ZETA2;
+use crate::harmonics::polygamma::cern_polygamma;
+
+/// Compute the harmonic sum $S_2(N)$.
+///
+/// $$S_2(N) = \sum\limits_{j=1}^N \frac 1 {j^2} = -\psi_1(N+1)+\zeta(2)$$
+/// with $\psi_1(N)$ the trigamma function and $\zeta$ the Riemann zeta function.
+pub fn S2(N: Complex<f64>) -> Complex<f64> {
+    -cern_polygamma(N + 1.0, 1) + ZETA2
+}
diff --git a/crates/ekore/src/harmonics/w3.rs b/crates/ekore/src/harmonics/w3.rs
new file mode 100644
index 000000000..7850c678d
--- /dev/null
+++ b/crates/ekore/src/harmonics/w3.rs
@@ -0,0 +1,13 @@
+//! Harmonic sums of weight 3.
+use num::complex::Complex;
+
+use crate::constants::ZETA3;
+use crate::harmonics::polygamma::cern_polygamma;
+
+/// Compute the harmonic sum $S_3(N)$.
+///
+/// $$S_3(N) = \sum\limits_{j=1}^N \frac 1 {j^3} = \frac 1 2 \psi_2(N+1)+\zeta(3)$$
+/// with $\psi_2(N)$ the 2nd polygamma function and $\zeta$ the Riemann zeta function.
+pub fn S3(N: Complex<f64>) -> Complex<f64> {
+    0.5 * cern_polygamma(N + 1.0, 2) + ZETA3
+}
diff --git a/crates/ekore/src/harmonics/w4.rs b/crates/ekore/src/harmonics/w4.rs
new file mode 100644
index 000000000..66835ac5f
--- /dev/null
+++ b/crates/ekore/src/harmonics/w4.rs
@@ -0,0 +1,13 @@
+//! Harmonic sums of weight 4.
+use num::complex::Complex;
+
+use crate::constants::ZETA4;
+use crate::harmonics::polygamma::cern_polygamma;
+
+/// Compute the harmonic sum $S_4(N)$.
+///
+/// $$S_4(N) = \sum\limits_{j=1}^N \frac 1 {j^4} = - \frac 1 6 \psi_3(N+1)+\zeta(4)$$
+/// with $\psi_3(N)$ the 3rd polygamma function and $\zeta$ the Riemann zeta function.
+pub fn S4(N: Complex<f64>) -> Complex<f64> {
+    ZETA4 - 1.0 / 6.0 * cern_polygamma(N + 1.0, 3)
+}
diff --git a/crates/ekore/src/util.rs b/crates/ekore/src/util.rs
index 5c4e34f4d..55b1ab2bd 100644
--- a/crates/ekore/src/util.rs
+++ b/crates/ekore/src/util.rs
@@ -7,3 +7,23 @@ macro_rules! cmplx {
         Complex::new($re, $im)
     };
 }
+
+/// Shorthand complex number contructor.
+#[macro_export]
+macro_rules! assert_approx_eq_cmplx {
+    ($size:ty, $ref:expr, $target:expr) => {
+        use float_cmp::assert_approx_eq;
+        assert_approx_eq!($size, $ref.re, $target.re);
+        assert_approx_eq!($size, $ref.im, $target.im);
+    };
+    ($size:ty, $ref:expr, $target:expr, ulps=$ulps:expr) => {
+        use float_cmp::assert_approx_eq;
+        assert_approx_eq!($size, $ref.re, $target.re, ulps = $ulps);
+        assert_approx_eq!($size, $ref.im, $target.im, ulps = $ulps);
+    };
+    ($size:ty, $ref:expr, $target:expr, epsilon=$epsilon:expr) => {
+        use float_cmp::assert_approx_eq;
+        assert_approx_eq!($size, $ref.re, $target.re, epsilon = $epsilon);
+        assert_approx_eq!($size, $ref.im, $target.im, epsilon = $epsilon);
+    };
+}
diff --git a/crates/make_bib.py b/crates/make_bib.py
index b085414c9..df5437a04 100644
--- a/crates/make_bib.py
+++ b/crates/make_bib.py
@@ -21,6 +21,7 @@
 /// {doi}"""
 # Combine publication information
 PUB = """Published in: {journal} {volume} ({year}), {pages}"""
+PHD = """Published as PhD thesis at {school} ({year})"""
 
 
 def clean_nl(t: str) -> str:
@@ -39,12 +40,18 @@ def clean_nl(t: str) -> str:
 for el in bib_database.entries:
     title = re.sub(r"^\{(.+)\}$", r"\1", clean_nl(el.fields_dict["title"].value))
     author = el.fields_dict["author"].value
-    publication = PUB.format(
-        journal=el.fields_dict["journal"].value,
-        volume=el.fields_dict["volume"].value,
-        year=el.fields_dict["year"].value,
-        pages=el.fields_dict["pages"].value,
-    )
+    if el.entry_type == "phdthesis":
+        publication = PHD.format(
+            school=el.fields_dict["school"].value,
+            year=el.fields_dict["year"].value,
+        )
+    else:
+        publication = PUB.format(
+            journal=el.fields_dict["journal"].value,
+            volume=el.fields_dict["volume"].value,
+            year=el.fields_dict["year"].value,
+            pages=el.fields_dict["pages"].value,
+        )
     eprint = ""
     if (
         "eprint" in el.fields_dict
diff --git a/flake.lock b/flake.lock
index d66f4ed08..d243cec36 100644
--- a/flake.lock
+++ b/flake.lock
@@ -197,6 +197,22 @@
         "type": "github"
       }
     },
+    "nixpkgs_2": {
+      "locked": {
+        "lastModified": 1716330097,
+        "narHash": "sha256-8BO3B7e3BiyIDsaKA0tY8O88rClYRTjvAp66y+VBUeU=",
+        "owner": "NixOS",
+        "repo": "nixpkgs",
+        "rev": "5710852ba686cc1fd0d3b8e22b3117d43ba374c2",
+        "type": "github"
+      },
+      "original": {
+        "owner": "NixOS",
+        "ref": "nixos-unstable",
+        "repo": "nixpkgs",
+        "type": "github"
+      }
+    },
     "pre-commit-hooks": {
       "inputs": {
         "flake-compat": [
@@ -229,7 +245,7 @@
       "inputs": {
         "devenv": "devenv",
         "flake-parts": "flake-parts",
-        "nixpkgs": "nixpkgs"
+        "nixpkgs": "nixpkgs_2"
       }
     },
     "systems": {
diff --git a/flake.nix b/flake.nix
index ffed59900..28cacfc7c 100644
--- a/flake.nix
+++ b/flake.nix
@@ -21,15 +21,18 @@
 
       perSystem = {pkgs, ...}: {
         devenv.shells.default = {
-          packages = with pkgs; [pre-commit poethepoet maturin];
+          packages = with pkgs; [maturin poethepoet pre-commit stdenv.cc.cc.lib];
 
           languages = {
             python = {
               enable = true;
               poetry = {
                 enable = true;
-                install.enable = true;
-                install.allExtras = true;
+                install = {
+                  enable = true;
+                  allExtras = true;
+                  groups = ["dev" "test"];
+                };
               };
             };
             rust.enable = true;
diff --git a/rustify.sh b/rustify.sh
index 7e99996f0..297afdb3d 100755
--- a/rustify.sh
+++ b/rustify.sh
@@ -1,4 +1,4 @@
-#!/usr/bin/bash
+#!/usr/bin/env bash
 
 # git diff --merge-base master pyproject.toml > pyproject.toml.patch
 patch -p1 <pyproject.toml.patch
diff --git a/src/eko/evolution_operator/quad_ker.py b/src/eko/evolution_operator/quad_ker.py
index d7b5858b0..798cebe6b 100644
--- a/src/eko/evolution_operator/quad_ker.py
+++ b/src/eko/evolution_operator/quad_ker.py
@@ -68,6 +68,7 @@ def select_singlet_element(ker, mode0, mode1):
         nb.types.bool_,  # is_threshold
     ),
     cache=True,
+    nopython=True,
 )
 def cb_quad_ker_qcd(
     re_gamma_raw,
@@ -103,7 +104,7 @@ def cb_quad_ker_qcd(
     areas = nb.carray(areas_raw, (areas_x, areas_y))
     pj = interpolation.evaluate_grid(n, is_log, logx, areas)
     # TODO recover parameters
-    method = "iterate-expanded"
+    method = "iterate-exact"
     sv_mode = sv.Modes.exponentiated
     order = (order_qcd, 0)
     ev_op_max_order = (ev_op_max_order_qcd, 0)
diff --git a/src/ekore/anomalous_dimensions/unpolarized/space_like/as2.py b/src/ekore/anomalous_dimensions/unpolarized/space_like/as2.py
index 01a115c5b..28f8533f4 100644
--- a/src/ekore/anomalous_dimensions/unpolarized/space_like/as2.py
+++ b/src/ekore/anomalous_dimensions/unpolarized/space_like/as2.py
@@ -18,7 +18,7 @@
 def gamma_nsm(n, nf, cache):
     r"""Compute the |NLO| valence-like non-singlet anomalous dimension.
 
-    Implements Eq. (3.5) of :cite:`Moch:2004pa`.
+    Implements Eq. (3.6) of :cite:`Moch:2004pa`.
 
     Parameters
     ----------