diff --git a/CHANGELOG.md b/CHANGELOG.md
index 0830525be..ccad7bb74 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -16,6 +16,7 @@ Items without prefix refer to a global change.
 ## [Unreleased](https://github.com/NNPDF/eko/compare/v0.15.4...HEAD)
 
 ### Changed
+- rust: Visibility of functions in `ekore` crate changed, `constants.rs` split into 3 files ([#543](https://github.com/NNPDF/eko/pull/543))
 - rust: Bumped `rust-version` to 1.85.0, `edition` to 2024 ([#545](https://github.com/NNPDF/eko/pull/545))
 
 ## [0.15.4](https://github.com/NNPDF/eko/compare/v0.15.3...v0.15.4) - 2026-06-24
diff --git a/crates/ekore/src/anomalous_dimensions/polarized/spacelike.rs b/crates/ekore/src/anomalous_dimensions/polarized/spacelike.rs
index 99c5dd858..2a6dedba2 100644
--- a/crates/ekore/src/anomalous_dimensions/polarized/spacelike.rs
+++ b/crates/ekore/src/anomalous_dimensions/polarized/spacelike.rs
@@ -4,8 +4,8 @@ use crate::constants::{MAX_ORDER_QCD, PID_NSM, PID_NSP, PID_NSV};
 use crate::harmonics::cache::Cache;
 use num::Zero;
 use num::complex::Complex;
-pub mod as1;
-pub mod as2;
+mod as1;
+mod as2;
 // pub mod as3;
 
 /// Compute the tower of the non-singlet anomalous dimensions.
diff --git a/crates/ekore/src/anomalous_dimensions/polarized/spacelike/as1.rs b/crates/ekore/src/anomalous_dimensions/polarized/spacelike/as1.rs
index 1d12335cd..629044404 100644
--- a/crates/ekore/src/anomalous_dimensions/polarized/spacelike/as1.rs
+++ b/crates/ekore/src/anomalous_dimensions/polarized/spacelike/as1.rs
@@ -9,14 +9,14 @@ use crate::harmonics::cache::{Cache, K};
 /// Compute the non-singlet anomalous dimension.
 ///
 /// Identical to the unpolarized counterpart.
-pub fn gamma_ns(c: &mut Cache, nf: u8) -> Complex<f64> {
+pub(super) fn gamma_ns(c: &mut Cache, nf: u8) -> Complex<f64> {
     unpol(c, nf)
 }
 
 /// Compute the quark-gluon anomalous dimension.
 ///
 /// Implements Eq. (A.1) of [\[Gluck:1995yr\]][crate::bib::Gluck1995yr].
-pub fn gamma_qg(c: &mut Cache, nf: u8) -> Complex<f64> {
+pub(super) fn gamma_qg(c: &mut Cache, nf: u8) -> Complex<f64> {
     let N = c.n();
     let gamma = -(N - 1.) / N / (N + 1.);
     2.0 * TR * 2.0 * (nf as f64) * gamma
@@ -25,7 +25,7 @@ pub fn gamma_qg(c: &mut Cache, nf: u8) -> Complex<f64> {
 /// Compute the gluon-quark anomalous dimension.
 ///
 /// Implements Eq. (A.1) of [\[Gluck:1995yr\]][crate::bib::Gluck1995yr].
-pub fn gamma_gq(c: &mut Cache, _nf: u8) -> Complex<f64> {
+pub(super) fn gamma_gq(c: &mut Cache, _nf: u8) -> Complex<f64> {
     let N = c.n();
     let gamma = -(N + 2.) / N / (N + 1.);
     2.0 * CF * gamma
@@ -34,7 +34,7 @@ pub fn gamma_gq(c: &mut Cache, _nf: u8) -> Complex<f64> {
 /// Compute the gluon-gluon anomalous dimension.
 ///
 /// Implements Eq. (A.1) of [\[Gluck:1995yr\]][crate::bib::Gluck1995yr].
-pub fn gamma_gg(c: &mut Cache, nf: u8) -> Complex<f64> {
+pub(super) fn gamma_gg(c: &mut Cache, nf: u8) -> Complex<f64> {
     let N = c.n();
     let S1 = c.get(K::S1);
     let gamma = -S1 + 2. / N / (N + 1.);
@@ -42,7 +42,7 @@ pub fn gamma_gg(c: &mut Cache, nf: u8) -> Complex<f64> {
 }
 
 /// Compute the singlet anomalous dimension matrix.
-pub fn gamma_singlet(c: &mut Cache, nf: u8) -> [[Complex<f64>; 2]; 2] {
+pub(super) fn gamma_singlet(c: &mut Cache, nf: u8) -> [[Complex<f64>; 2]; 2] {
     let gamma_qq = gamma_ns(c, nf);
     [
         [gamma_qq, gamma_qg(c, nf)],
diff --git a/crates/ekore/src/anomalous_dimensions/polarized/spacelike/as2.rs b/crates/ekore/src/anomalous_dimensions/polarized/spacelike/as2.rs
index d29d48932..ee1714e2e 100644
--- a/crates/ekore/src/anomalous_dimensions/polarized/spacelike/as2.rs
+++ b/crates/ekore/src/anomalous_dimensions/polarized/spacelike/as2.rs
@@ -7,18 +7,18 @@ use super::super::super::unpolarized::spacelike::as2::gamma_nsp as unpol_nsp;
 use crate::constants::{CA, CF, TR, ZETA2, ZETA3};
 use crate::harmonics::cache::{Cache, K};
 
-pub fn gamma_nsm(c: &mut Cache, nf: u8) -> Complex<f64> {
+pub(super) fn gamma_nsm(c: &mut Cache, nf: u8) -> Complex<f64> {
     unpol_nsp(c, nf)
 }
 
-pub fn gamma_nsp(c: &mut Cache, nf: u8) -> Complex<f64> {
+pub(super) fn gamma_nsp(c: &mut Cache, nf: u8) -> Complex<f64> {
     unpol_nsm(c, nf)
 }
 
 /// Compute the pure-singlet quark-quark anomalous dimension.
 ///
 /// Implements Eq. (A.3) of [\[Gluck:1995yr\]][crate::bib::Gluck1995yr].
-pub fn gamma_ps(c: &mut Cache, nf: u8) -> Complex<f64> {
+pub(super) fn gamma_ps(c: &mut Cache, nf: u8) -> Complex<f64> {
     let n = c.n();
     let gqqps1_nfcf = (2. * (n + 2.) * (1. + 2. * n + n.powu(3))) / ((1. + n).powu(3) * n.powu(3));
     4.0 * TR * (nf as f64) * CF * gqqps1_nfcf
@@ -27,7 +27,7 @@ pub fn gamma_ps(c: &mut Cache, nf: u8) -> Complex<f64> {
 /// Compute the quark-gluon singlet anomalous dimension.
 ///
 /// Implements Eq. (A.4) of [\[Gluck:1995yr\]][crate::bib::Gluck1995yr].
-pub fn gamma_qg(c: &mut Cache, nf: u8) -> Complex<f64> {
+pub(super) 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);
@@ -53,7 +53,7 @@ pub fn gamma_qg(c: &mut Cache, nf: u8) -> Complex<f64> {
 /// Compute the gluon-quark singlet anomalous dimension.
 ///
 /// Implements Eq. (A.5) of [\[Gluck:1995yr\]][crate::bib::Gluck1995yr].
-pub fn gamma_gq(c: &mut Cache, nf: u8) -> Complex<f64> {
+pub(super) 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);
@@ -83,7 +83,7 @@ pub fn gamma_gq(c: &mut Cache, nf: u8) -> Complex<f64> {
 /// Compute the gluon-gluon singlet anomalous dimension.
 ///
 /// Implements Eq. (A.6) of [\[Gluck:1995yr\]][crate::bib::Gluck1995yr].
-pub fn gamma_gg(c: &mut Cache, nf: u8) -> Complex<f64> {
+pub(super) fn gamma_gg(c: &mut Cache, nf: u8) -> Complex<f64> {
     let n = c.n();
     let S1 = c.get(K::S1);
     let Sp1m = c.get(K::S1mh);
@@ -118,7 +118,7 @@ pub fn gamma_gg(c: &mut Cache, nf: u8) -> Complex<f64> {
 }
 
 /// Compute the singlet anomalous dimension matrix.
-pub fn gamma_singlet(c: &mut Cache, nf: u8) -> [[Complex<f64>; 2]; 2] {
+pub(super) 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)],
diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs
index e7e92c3e6..aab033980 100644
--- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs
+++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs
@@ -7,13 +7,13 @@ use crate::constants::{
 use crate::harmonics::cache::Cache;
 use num::Zero;
 use num::complex::Complex;
-pub mod aem1;
-pub mod aem2;
-pub mod as1;
-pub mod as1aem1;
-pub mod as2;
-pub mod as3;
-pub mod as4;
+mod aem1;
+mod aem2;
+pub(crate) mod as1;
+mod as1aem1;
+pub(crate) mod as2;
+mod as3;
+mod as4;
 
 /// Compute the tower of the non-singlet anomalous dimensions.
 ///
diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem1.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem1.rs
index 74976c938..67dcf9c55 100644
--- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem1.rs
+++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem1.rs
@@ -11,7 +11,7 @@ use crate::anomalous_dimensions::unpolarized::spacelike::as1;
 /// Compute the photon-quark anomalous dimension.
 ///
 /// Implements Eq. (2.5) of [\[Carrazza:2015dea\]][crate::bib::Carrazza2015dea].
-pub fn gamma_phq(c: &mut Cache, nf: u8) -> Complex<f64> {
+pub(super) fn gamma_phq(c: &mut Cache, nf: u8) -> Complex<f64> {
     as1::gamma_gq(c, nf) / CF
 }
 
@@ -20,14 +20,14 @@ pub fn gamma_phq(c: &mut Cache, nf: u8) -> Complex<f64> {
 /// Implements Eq. (2.5) of [\[Carrazza:2015dea\]][crate::bib::Carrazza2015dea].
 /// However, we are adding the $N_C$ and the $2n_f$ factors from $\theta$ inside the
 /// definition of $\gamma_{q \gamma}^{(0)}(N)$.
-pub fn gamma_qph(c: &mut Cache, nf: u8) -> Complex<f64> {
+pub(super) fn gamma_qph(c: &mut Cache, nf: u8) -> Complex<f64> {
     as1::gamma_qg(c, nf) / TR * (NC as f64)
 }
 
 /// Compute the photon-photon anomalous dimension.
 ///
 /// Implements Eq. (2.5) of [\[Carrazza:2015dea\]][crate::bib::Carrazza2015dea].
-pub fn gamma_phph(_c: &mut Cache, nf: u8) -> Complex<f64> {
+pub(super) fn gamma_phph(_c: &mut Cache, nf: u8) -> Complex<f64> {
     let cc = ChargeCombinations { nf };
     cmplx!(
         4.0 / 3.0 * (NC as f64) * ((cc.nu() as f64) * EU2 + (cc.nd() as f64) * ED2),
@@ -38,12 +38,12 @@ pub fn gamma_phph(_c: &mut Cache, nf: u8) -> Complex<f64> {
 /// Compute the non-singlet anomalous dimension.
 ///
 /// Implements Eq. (2.5) of [\[Carrazza:2015dea\]][crate::bib::Carrazza2015dea].
-pub fn gamma_ns(c: &mut Cache, nf: u8) -> Complex<f64> {
+pub(super) fn gamma_ns(c: &mut Cache, nf: u8) -> Complex<f64> {
     as1::gamma_ns(c, nf) / CF
 }
 
 /// Compute the singlet anomalous dimension matrix.
-pub fn gamma_singlet(c: &mut Cache, nf: u8) -> [[Complex<f64>; 4]; 4] {
+pub(super) fn gamma_singlet(c: &mut Cache, nf: u8) -> [[Complex<f64>; 4]; 4] {
     let cc = ChargeCombinations { nf };
 
     let gamma_ph_q = gamma_phq(c, nf);
@@ -79,7 +79,7 @@ pub fn gamma_singlet(c: &mut Cache, nf: u8) -> [[Complex<f64>; 4]; 4] {
 }
 
 /// Compute the valence anomalous dimension matrix.
-pub fn gamma_valence(c: &mut Cache, nf: u8) -> [[Complex<f64>; 2]; 2] {
+pub(super) fn gamma_valence(c: &mut Cache, nf: u8) -> [[Complex<f64>; 2]; 2] {
     let cc = ChargeCombinations { nf };
     let g = gamma_ns(c, nf);
     [
diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem2.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem2.rs
index 5db88cbed..bcb7bf306 100644
--- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem2.rs
+++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem2.rs
@@ -10,7 +10,7 @@ use crate::anomalous_dimensions::unpolarized::spacelike::as1aem1;
 /// Compute the photon-photon anomalous dimension.
 ///
 /// Implements Eq. (68) of [\[deFlorian:2016gvk\]][crate::bib::deFlorian2016gvk].
-pub fn gamma_phph(c: &mut Cache, nf: u8) -> Complex<f64> {
+pub(super) fn gamma_phph(c: &mut Cache, nf: u8) -> Complex<f64> {
     let cc = ChargeCombinations { nf };
     (NC as f64)
         * ((cc.nu() as f64) * EU2.powi(2) + (cc.nd() as f64) * ED2.powi(2))
@@ -20,14 +20,14 @@ pub fn gamma_phph(c: &mut Cache, nf: u8) -> Complex<f64> {
 /// Compute the up-photon anomalous dimension.
 ///
 /// Implements Eq. (55) of [\[deFlorian:2016gvk\]][crate::bib::deFlorian2016gvk] for q=u.
-pub fn gamma_uph(c: &mut Cache, nf: u8) -> Complex<f64> {
+pub(super) fn gamma_uph(c: &mut Cache, nf: u8) -> Complex<f64> {
     EU2 * as1aem1::gamma_qph(c, nf) / CF
 }
 
 /// Compute the down-photon anomalous dimension.
 ///
 /// Implements Eq. (55) of [\[deFlorian:2016gvk\]][crate::bib::deFlorian2016gvk] for q=d.
-pub fn gamma_dph(c: &mut Cache, nf: u8) -> Complex<f64> {
+pub(super) fn gamma_dph(c: &mut Cache, nf: u8) -> Complex<f64> {
     ED2 * as1aem1::gamma_qph(c, nf) / CF
 }
 
@@ -50,14 +50,14 @@ fn gamma_phq(c: &mut Cache, nf: u8) -> Complex<f64> {
 /// Compute the photon-up anomalous dimension.
 ///
 /// Implements Eq. (56) of [\[deFlorian:2016gvk\]][crate::bib::deFlorian2016gvk].
-pub fn gamma_phu(c: &mut Cache, nf: u8) -> Complex<f64> {
+pub(super) fn gamma_phu(c: &mut Cache, nf: u8) -> Complex<f64> {
     EU2 * as1aem1::gamma_phq(c, nf) / CF + gamma_phq(c, nf)
 }
 
 /// Compute the photon-down anomalous dimension.
 ///
 /// Implements Eq. (56) of [\[deFlorian:2016gvk\]][crate::bib::deFlorian2016gvk].
-pub fn gamma_phd(c: &mut Cache, nf: u8) -> Complex<f64> {
+pub(super) fn gamma_phd(c: &mut Cache, nf: u8) -> Complex<f64> {
     ED2 * as1aem1::gamma_phq(c, nf) / CF + gamma_phq(c, nf)
 }
 
@@ -80,35 +80,35 @@ fn gamma_nsq(c: &mut Cache, nf: u8) -> Complex<f64> {
 /// Compute the singlet-like non-singlet up anomalous dimension.
 ///
 /// Implements sum of Eqs. (57-58) of [\[deFlorian:2016gvk\]][crate::bib::deFlorian2016gvk] for q=u.
-pub fn gamma_nspu(c: &mut Cache, nf: u8) -> Complex<f64> {
+pub(super) fn gamma_nspu(c: &mut Cache, nf: u8) -> Complex<f64> {
     EU2 * as1aem1::gamma_nsp(c, nf) / CF / 2. + gamma_nsq(c, nf)
 }
 
 /// Compute the singlet-like non-singlet down anomalous dimension.
 ///
 /// Implements sum of Eqs. (57-58) of [\[deFlorian:2016gvk\]][crate::bib::deFlorian2016gvk] for q=d.
-pub fn gamma_nspd(c: &mut Cache, nf: u8) -> Complex<f64> {
+pub(super) fn gamma_nspd(c: &mut Cache, nf: u8) -> Complex<f64> {
     ED2 * as1aem1::gamma_nsp(c, nf) / CF / 2. + gamma_nsq(c, nf)
 }
 
 /// Compute the valence-like non-singlet up anomalous dimension.
 ///
 /// Implements difference between Eqs. (57-58) of [\[deFlorian:2016gvk\]][crate::bib::deFlorian2016gvk] for q=u.
-pub fn gamma_nsmu(c: &mut Cache, nf: u8) -> Complex<f64> {
+pub(super) fn gamma_nsmu(c: &mut Cache, nf: u8) -> Complex<f64> {
     EU2 * as1aem1::gamma_nsm(c, nf) / CF / 2. + gamma_nsq(c, nf)
 }
 
 /// Compute the valence-like non-singlet down anomalous dimension.
 ///
 /// Implements difference between Eqs. (57-58) of [\[deFlorian:2016gvk\]][crate::bib::deFlorian2016gvk] for q=d.
-pub fn gamma_nsmd(c: &mut Cache, nf: u8) -> Complex<f64> {
+pub(super) fn gamma_nsmd(c: &mut Cache, nf: u8) -> Complex<f64> {
     ED2 * as1aem1::gamma_nsm(c, nf) / CF / 2. + gamma_nsq(c, nf)
 }
 
 /// Compute the pure-singlet anomalous dimension.
 ///
 /// Implements Eq. (59) of [\[deFlorian:2016gvk\]][crate::bib::deFlorian2016gvk].
-pub fn gamma_ps(c: &mut Cache, nf: u8) -> Complex<f64> {
+pub(super) fn gamma_ps(c: &mut Cache, nf: u8) -> Complex<f64> {
     let N = c.n();
     let result = -4.0 * (2.0 + N * (5.0 + N)) * (4.0 + N * (4.0 + N * (7.0 + 5.0 * N)))
         / ((-1.0 + N) * N.powu(3) * (1.0 + N).powu(3) * (2.0 + N).powu(2));
@@ -116,7 +116,7 @@ pub fn gamma_ps(c: &mut Cache, nf: u8) -> Complex<f64> {
 }
 
 /// Compute the singlet anomalous dimension matrix.
-pub fn gamma_singlet(c: &mut Cache, nf: u8) -> [[Complex<f64>; 4]; 4] {
+pub(super) fn gamma_singlet(c: &mut Cache, nf: u8) -> [[Complex<f64>; 4]; 4] {
     let cc = ChargeCombinations { nf };
     const E2M: f64 = EU2 - ED2;
 
@@ -159,7 +159,7 @@ pub fn gamma_singlet(c: &mut Cache, nf: u8) -> [[Complex<f64>; 4]; 4] {
 }
 
 /// Compute the valence anomalous dimension matrix.
-pub fn gamma_valence(c: &mut Cache, nf: u8) -> [[Complex<f64>; 2]; 2] {
+pub(super) fn gamma_valence(c: &mut Cache, nf: u8) -> [[Complex<f64>; 2]; 2] {
     let cc = ChargeCombinations { nf };
     let gamma_nsmu = gamma_nsmu(c, nf);
     let gamma_nsmd = gamma_nsmd(c, nf);
diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1.rs
index 1b9bb5c65..f5d3a6190 100644
--- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1.rs
+++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1.rs
@@ -8,7 +8,7 @@ use crate::harmonics::cache::{Cache, K};
 /// Compute the non-singlet anomalous dimension.
 ///
 /// Implements Eq. (3.4) of [\[Moch:2004pa\]][crate::bib::Moch2004pa].
-pub fn gamma_ns(c: &mut Cache, _nf: u8) -> Complex<f64> {
+pub(crate) fn gamma_ns(c: &mut Cache, _nf: u8) -> Complex<f64> {
     let N = c.n();
     let S1 = c.get(K::S1);
     let gamma = -(3.0 - 4.0 * S1 + 2.0 / N / (N + 1.0));
@@ -18,7 +18,7 @@ pub fn gamma_ns(c: &mut Cache, _nf: u8) -> Complex<f64> {
 /// Compute the quark-gluon anomalous dimension.
 ///
 /// Implements Eq. (3.5) of [\[Vogt:2004mw\]](crate::bib::Vogt2004mw).
-pub fn gamma_qg(c: &mut Cache, nf: u8) -> Complex<f64> {
+pub(super) fn gamma_qg(c: &mut Cache, nf: u8) -> Complex<f64> {
     let N = c.n();
     let gamma = -(N.powu(2) + N + 2.0) / (N * (N + 1.0) * (N + 2.0));
     2.0 * TR * 2.0 * (nf as f64) * gamma
@@ -27,7 +27,7 @@ pub fn gamma_qg(c: &mut Cache, nf: u8) -> Complex<f64> {
 /// Compute the gluon-quark anomalous dimension.
 ///
 /// Implements Eq. (3.5) of [\[Vogt:2004mw\]](crate::bib::Vogt2004mw).
-pub fn gamma_gq(c: &mut Cache, _nf: u8) -> Complex<f64> {
+pub(super) fn gamma_gq(c: &mut Cache, _nf: u8) -> Complex<f64> {
     let N = c.n();
     let gamma = -(N.powu(2) + N + 2.0) / (N * (N + 1.0) * (N - 1.0));
     2.0 * CF * gamma
@@ -36,7 +36,7 @@ pub fn gamma_gq(c: &mut Cache, _nf: u8) -> Complex<f64> {
 /// Compute the gluon-gluon anomalous dimension.
 ///
 /// Implements Eq. (3.5) of [\[Vogt:2004mw\]](crate::bib::Vogt2004mw).
-pub fn gamma_gg(c: &mut Cache, nf: u8) -> Complex<f64> {
+pub(super) fn gamma_gg(c: &mut Cache, nf: u8) -> Complex<f64> {
     let N = c.n();
     let S1 = c.get(K::S1);
     let gamma = S1 - 1.0 / N / (N - 1.0) - 1.0 / (N + 1.0) / (N + 2.0);
@@ -44,7 +44,7 @@ pub fn gamma_gg(c: &mut Cache, nf: u8) -> Complex<f64> {
 }
 
 /// Compute the singlet anomalous dimension matrix.
-pub fn gamma_singlet(c: &mut Cache, nf: u8) -> [[Complex<f64>; 2]; 2] {
+pub(super) fn gamma_singlet(c: &mut Cache, nf: u8) -> [[Complex<f64>; 2]; 2] {
     let gamma_qq = gamma_ns(c, nf);
     [
         [gamma_qq, gamma_qg(c, nf)],
diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs
index 487f44965..2fffb3cc9 100644
--- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs
+++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs
@@ -9,7 +9,7 @@ use std::f64::consts::PI;
 /// Compute the photon-quark anomalous dimension.
 ///
 /// Implements Eq. (36) of [\[deFlorian:2015ujt\]][crate::bib::deFlorian2015ujt].
-pub fn gamma_phq(c: &mut Cache, _nf: u8) -> Complex<f64> {
+pub(super) fn gamma_phq(c: &mut Cache, _nf: u8) -> Complex<f64> {
     let N = c.n();
     let S1 = c.get(K::S1);
     let S2 = c.get(K::S2);
@@ -41,7 +41,7 @@ pub fn gamma_phq(c: &mut Cache, _nf: u8) -> Complex<f64> {
 /// Compute the quark-photon anomalous dimension.
 ///
 /// Implements Eq. (26) of [\[deFlorian:2015ujt\]][crate::bib::deFlorian2015ujt].
-pub fn gamma_qph(c: &mut Cache, nf: u8) -> Complex<f64> {
+pub(super) fn gamma_qph(c: &mut Cache, nf: u8) -> Complex<f64> {
     let N = c.n();
     let S1 = c.get(K::S1);
     let S2 = c.get(K::S2);
@@ -67,7 +67,7 @@ pub fn gamma_qph(c: &mut Cache, nf: u8) -> Complex<f64> {
 /// Compute the gluon-photon anomalous dimension.
 ///
 /// Implements Eq. (27) of [\[deFlorian:2015ujt\]][crate::bib::deFlorian2015ujt].
-pub fn gamma_gph(c: &mut Cache, _nf: u8) -> Complex<f64> {
+pub(super) fn gamma_gph(c: &mut Cache, _nf: u8) -> Complex<f64> {
     let N = c.n();
     CF * CA
         * (8.0 * (-4.0 + N * (-4.0 + N * (-5.0 + N * (-10.0 + N + 2.0 * N.powu(2) * (2.0 + N))))))
@@ -77,28 +77,28 @@ pub fn gamma_gph(c: &mut Cache, _nf: u8) -> Complex<f64> {
 /// Compute the photon-gluon anomalous dimension.
 ///
 /// Implements Eq. (30) of [\[deFlorian:2015ujt\]][crate::bib::deFlorian2015ujt].
-pub fn gamma_phg(c: &mut Cache, nf: u8) -> Complex<f64> {
+pub(super) fn gamma_phg(c: &mut Cache, nf: u8) -> Complex<f64> {
     TR / CF / CA * (NC as f64) * gamma_gph(c, nf)
 }
 
 /// Compute the quark-gluon singlet anomalous dimension.
 ///
 /// Implements Eq. (29) of [\[deFlorian:2015ujt\]][crate::bib::deFlorian2015ujt].
-pub fn gamma_qg(c: &mut Cache, nf: u8) -> Complex<f64> {
+pub(super) fn gamma_qg(c: &mut Cache, nf: u8) -> Complex<f64> {
     TR / CF / CA * (NC as f64) * gamma_qph(c, nf)
 }
 
 /// Compute the gluon-quark singlet anomalous dimension.
 ///
 /// Implements Eq. (35) of [\[deFlorian:2015ujt\]][crate::bib::deFlorian2015ujt].
-pub fn gamma_gq(c: &mut Cache, nf: u8) -> Complex<f64> {
+pub(super) fn gamma_gq(c: &mut Cache, nf: u8) -> Complex<f64> {
     gamma_phq(c, nf)
 }
 
 /// Compute the photon-photon singlet anomalous dimension.
 ///
 /// Implements Eq. (28) of [\[deFlorian:2015ujt\]][crate::bib::deFlorian2015ujt].
-pub fn gamma_phph(_c: &mut Cache, nf: u8) -> Complex<f64> {
+pub(super) fn gamma_phph(_c: &mut Cache, nf: u8) -> Complex<f64> {
     let cc = ChargeCombinations { nf };
     cmplx!(
         4.0 * CF * CA * ((cc.nu() as f64) * EU2 + (cc.nd() as f64) * ED2),
@@ -109,7 +109,7 @@ pub fn gamma_phph(_c: &mut Cache, nf: u8) -> Complex<f64> {
 /// Compute the gluon-gluon singlet anomalous dimension.
 ///
 /// Implements Eq. (31) of [\[deFlorian:2015ujt\]][crate::bib::deFlorian2015ujt].
-pub fn gamma_gg(_c: &mut Cache, _nf: u8) -> Complex<f64> {
+pub(super) fn gamma_gg(_c: &mut Cache, _nf: u8) -> Complex<f64> {
     cmplx!(4.0 * TR * (NC as f64), 0.)
 }
 
@@ -126,7 +126,7 @@ fn g3_shift(c: &mut Cache) -> Complex<f64> {
 /// Compute the singlet-like non-singlet anomalous dimension.
 ///
 /// Implements Eqs. (33-34) of [\[deFlorian:2015ujt\]][crate::bib::deFlorian2015ujt].
-pub fn gamma_nsp(c: &mut Cache, _nf: u8) -> Complex<f64> {
+pub(super) 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);
@@ -158,7 +158,7 @@ pub fn gamma_nsp(c: &mut Cache, _nf: u8) -> Complex<f64> {
 /// Compute the valence-like non-singlet anomalous dimension.
 ///
 /// Implements Eqs. (33-34) of [\[deFlorian:2015ujt\]][crate::bib::deFlorian2015ujt].
-pub fn gamma_nsm(c: &mut Cache, _nf: u8) -> Complex<f64> {
+pub(super) 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);
@@ -206,7 +206,7 @@ pub fn gamma_nsm(c: &mut Cache, _nf: u8) -> Complex<f64> {
 }
 
 /// Compute the singlet anomalous dimension matrix.
-pub fn gamma_singlet(c: &mut Cache, nf: u8) -> [[Complex<f64>; 4]; 4] {
+pub(super) fn gamma_singlet(c: &mut Cache, nf: u8) -> [[Complex<f64>; 4]; 4] {
     let cc = ChargeCombinations { nf };
     let e2_tot = nf as f64 * cc.e2avg();
     let gq = gamma_gq(c, nf);
@@ -243,7 +243,7 @@ pub fn gamma_singlet(c: &mut Cache, nf: u8) -> [[Complex<f64>; 4]; 4] {
 }
 
 /// Compute the valence anomalous dimension matrix.
-pub fn gamma_valence(c: &mut Cache, nf: u8) -> [[Complex<f64>; 2]; 2] {
+pub(super) fn gamma_valence(c: &mut Cache, nf: u8) -> [[Complex<f64>; 2]; 2] {
     let cc = ChargeCombinations { nf };
     let g = gamma_nsm(c, nf);
     [
diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as2.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as2.rs
index df17fecee..7626eb0a0 100644
--- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as2.rs
+++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as2.rs
@@ -9,7 +9,7 @@ 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> {
+pub(crate) 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);
@@ -49,7 +49,7 @@ pub fn gamma_nsm(c: &mut Cache, nf: u8) -> Complex<f64> {
 /// 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> {
+pub(crate) 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);
@@ -89,7 +89,7 @@ pub fn gamma_nsp(c: &mut Cache, nf: u8) -> Complex<f64> {
 /// 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> {
+pub(super) 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));
@@ -99,7 +99,7 @@ pub fn gamma_ps(c: &mut Cache, nf: u8) -> Complex<f64> {
 /// 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> {
+pub(super) 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);
@@ -127,7 +127,7 @@ pub fn gamma_qg(c: &mut Cache, nf: u8) -> Complex<f64> {
 /// 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> {
+pub(super) 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);
@@ -163,7 +163,7 @@ pub fn gamma_gq(c: &mut Cache, nf: u8) -> Complex<f64> {
 /// 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> {
+pub(super) 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);
@@ -206,7 +206,7 @@ pub fn gamma_gg(c: &mut Cache, nf: u8) -> Complex<f64> {
 }
 
 /// Compute the singlet anomalous dimension matrix.
-pub fn gamma_singlet(c: &mut Cache, nf: u8) -> [[Complex<f64>; 2]; 2] {
+pub(super) 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)],
diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as3.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as3.rs
index de4eac4fc..f77263e7f 100644
--- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as3.rs
+++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as3.rs
@@ -10,7 +10,7 @@ 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> {
+pub(super) 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);
@@ -67,7 +67,7 @@ pub fn gamma_nsm(c: &mut Cache, nf: u8) -> Complex<f64> {
 /// 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> {
+pub(super) 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);
@@ -124,7 +124,7 @@ pub fn gamma_nsp(c: &mut Cache, nf: u8) -> Complex<f64> {
 /// 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> {
+pub(super) 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);
@@ -161,7 +161,7 @@ pub fn gamma_nsv(c: &mut Cache, nf: u8) -> Complex<f64> {
 /// 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> {
+pub(super) 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);
@@ -217,7 +217,7 @@ pub fn gamma_ps(c: &mut Cache, nf: u8) -> Complex<f64> {
 /// 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> {
+pub(super) 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);
@@ -274,7 +274,7 @@ pub fn gamma_qg(c: &mut Cache, nf: u8) -> Complex<f64> {
 /// 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> {
+pub(super) 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);
@@ -347,7 +347,7 @@ pub fn gamma_gq(c: &mut Cache, nf: u8) -> Complex<f64> {
 /// 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> {
+pub(super) 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);
@@ -415,7 +415,7 @@ pub fn gamma_gg(c: &mut Cache, nf: u8) -> Complex<f64> {
 }
 
 /// Compute the singlet anomalous dimension matrix.
-pub fn gamma_singlet(c: &mut Cache, nf: u8) -> [[Complex<f64>; 2]; 2] {
+pub(super) 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)],
diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as4.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as4.rs
index 7f21e14e2..389b7700b 100644
--- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as4.rs
+++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as4.rs
@@ -11,16 +11,12 @@ mod gqg;
 use crate::harmonics::cache::Cache;
 use num::complex::Complex;
 
-pub use ggg::gamma_gg;
-pub use ggq::gamma_gq;
-pub use gnsm::gamma_nsm;
-pub use gnsp::gamma_nsp;
-pub use gnsv::gamma_nsv;
-pub use gps::gamma_ps;
-pub use gqg::gamma_qg;
+pub(super) use gnsm::gamma_nsm;
+pub(super) use gnsp::gamma_nsp;
+pub(super) use gnsv::gamma_nsv;
 
 /// Compute the singlet anomalous dimension matrix.
-pub fn gamma_singlet(c: &mut Cache, nf: u8, variation: [u8; 4]) -> [[Complex<f64>; 2]; 2] {
+pub(super) fn gamma_singlet(c: &mut Cache, nf: u8, variation: [u8; 4]) -> [[Complex<f64>; 2]; 2] {
     let gamma_qq = gnsp::gamma_nsp(c, nf, variation[3]) + gps::gamma_ps(c, nf, variation[3]);
     [
         [gamma_qq, gqg::gamma_qg(c, nf, variation[2])],
diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as4/ggg.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as4/ggg.rs
index 8f02ec240..cc6ec34e1 100644
--- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as4/ggg.rs
+++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as4/ggg.rs
@@ -16,7 +16,7 @@ use crate::harmonics::log_functions::{
 /// first 10 even moments together with small-x/large-x constraints.
 /// The two sets providing the error estimate are called via `variation = 1`
 /// and `variation = 2`.  Any other value of `variation` invokes their average.
-pub fn gamma_gg(c: &mut Cache, nf: u8, variation: u8) -> Complex<f64> {
+pub(super) fn gamma_gg(c: &mut Cache, nf: u8, variation: u8) -> Complex<f64> {
     let n = c.n();
     let S1 = c.get(K::S1);
     let S2 = c.get(K::S2);
diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as4/ggq.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as4/ggq.rs
index ff7ead02e..949eec4c8 100644
--- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as4/ggq.rs
+++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as4/ggq.rs
@@ -14,7 +14,7 @@ use crate::harmonics::log_functions::{lm11, lm12, lm12m1, lm13, lm14, lm14m1, lm
 /// first 10 even moments together with small-x/large-x constraints.
 /// The two sets providing the error estimate are called via `variation = 1`
 /// and `variation = 2`.  Any other value of `variation` invokes their average.
-pub fn gamma_gq(c: &mut Cache, nf: u8, variation: u8) -> Complex<f64> {
+pub(super) fn gamma_gq(c: &mut Cache, nf: u8, variation: u8) -> Complex<f64> {
     let n = c.n();
     let S1 = c.get(K::S1);
     let S2 = c.get(K::S2);
diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as4/gnsm.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as4/gnsm.rs
index 507a42af7..17d011f3d 100644
--- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as4/gnsm.rs
+++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as4/gnsm.rs
@@ -19,7 +19,7 @@ use crate::harmonics::log_functions::{lm11, lm11m1, lm12m1, lm13m1};
 /// eight even moments together with small-x and large-x constraints.
 /// The two sets spanning the error estimate are called via `variation = 1`
 /// and  ``variation = 2``. Any other value of `variation` invokes their average.
-pub fn gamma_nsm(c: &mut Cache, nf: u8, variation: u8) -> Complex<f64> {
+pub(crate) fn gamma_nsm(c: &mut Cache, nf: u8, variation: u8) -> Complex<f64> {
     let n = c.n();
     let S1 = c.get(K::S1);
     let S2 = c.get(K::S2);
diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as4/gnsp.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as4/gnsp.rs
index fafe956ef..ef06aebab 100644
--- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as4/gnsp.rs
+++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as4/gnsp.rs
@@ -10,7 +10,7 @@ use crate::harmonics::log_functions::{lm11, lm11m1, lm12m1, lm13m1};
 /// Compute the singlet-like non-singlet anomalous dimension.
 ///
 /// See [gamma_nsm][super::gamma_nsm] for implementation details.
-pub fn gamma_nsp(c: &mut Cache, nf: u8, variation: u8) -> Complex<f64> {
+pub(crate) fn gamma_nsp(c: &mut Cache, nf: u8, variation: u8) -> Complex<f64> {
     let n = c.n();
     let S1 = c.get(K::S1);
     let S2 = c.get(K::S2);
diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as4/gnsv.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as4/gnsv.rs
index 07bb2bce5..59439f9c9 100644
--- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as4/gnsv.rs
+++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as4/gnsv.rs
@@ -70,7 +70,7 @@ fn gamma_nss(c: &mut Cache, nf: u8, variation: u8) -> Complex<f64> {
 /// Compute the valence non-singlet anomalous dimension.
 ///
 /// See [gamma_nsm][super::gamma_nsm] for implementation details.
-pub fn gamma_nsv(c: &mut Cache, nf: u8, variation: u8) -> Complex<f64> {
+pub(crate) fn gamma_nsv(c: &mut Cache, nf: u8, variation: u8) -> Complex<f64> {
     gnsm::gamma_nsm(c, nf, variation) + gamma_nss(c, nf, variation)
 }
 
diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as4/gps.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as4/gps.rs
index 70ae4d42b..a335e6f9c 100644
--- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as4/gps.rs
+++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as4/gps.rs
@@ -13,7 +13,7 @@ use crate::harmonics::log_functions::{lm11m1, lm12m1, lm12m2, lm13m1, lm13m2, lm
 /// first ten even moments together with small-x/large-x constraints.
 /// The two sets spanning the error estimate are called via `variation = 1`
 /// and `variation = 2`.  Any other value of `variation` invokes their average.
-pub fn gamma_ps(c: &mut Cache, nf: u8, variation: u8) -> Complex<f64> {
+pub(super) fn gamma_ps(c: &mut Cache, nf: u8, variation: u8) -> Complex<f64> {
     let n = c.n();
     let nf_ = nf as f64;
     let nf2 = nf_.pow(2);
diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as4/gqg.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as4/gqg.rs
index 14fb7a2ca..3a3009d61 100644
--- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as4/gqg.rs
+++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as4/gqg.rs
@@ -14,7 +14,7 @@ use crate::harmonics::log_functions::{lm11, lm12, lm13, lm14, lm14m1, lm15, lm15
 /// first 10 even moments together with small-x/large-x constraints.
 /// The two sets indicating the error estimate are called via `variation = 1`
 /// and `variation = 2`.  Any other value of `variation` invokes their average.
-pub fn gamma_qg(c: &mut Cache, nf: u8, variation: u8) -> Complex<f64> {
+pub(super) fn gamma_qg(c: &mut Cache, nf: u8, variation: u8) -> Complex<f64> {
     let n = c.n();
     let S1 = c.get(K::S1);
     let S2 = c.get(K::S2);
diff --git a/crates/ekore/src/constants.rs b/crates/ekore/src/constants.rs
index a8e1ad923..81d8ae360 100644
--- a/crates/ekore/src/constants.rs
+++ b/crates/ekore/src/constants.rs
@@ -1,121 +1,18 @@
 //! Global constants.
 
-/// The number of colors.
-///
-/// Defaults to $N_C = 3$.
-pub const NC: u8 = 3;
+#[path = "math_constants.rs"]
+pub(super) mod math_constants;
+#[path = "pid.rs"]
+pub mod pid;
+#[path = "qcd_constants.rs"]
+pub(super) mod qcd_constants;
 
-/// The normalization of fundamental generators.
-///
-/// Defaults to $T_R = 1/2$.
-pub const TR: f64 = 1.0 / 2.0;
-
-/// Second Casimir constant in the adjoint representation.
-///
-/// Defaults to $C_A = N_C = 3$.
-pub const CA: f64 = NC as f64;
-
-/// Second Casimir constant in the fundamental representation.
-///
-/// 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);
-
-/// Up quark electric charge square.
-///
-/// Defaults to $e_u^2 = 4./9$.
-pub const EU2: f64 = 4. / 9.;
-
-/// Down quark electric charge square.
-///
-/// Defaults to $e_d^2 = 1./9$.
-pub const ED2: f64 = 1. / 9.;
-
-/// 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;
-
-/// Riemann zeta function at z = 5.
-pub const ZETA5: f64 = 1.03692775514337;
-
-/// 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;
-
-/// singlet-like non-singlet up-sector |PID|.
-pub const PID_NSP_U: u16 = 10102;
-
-/// singlet-like non-singlet down-sector |PID|.
-pub const PID_NSP_D: u16 = 10103;
-
-/// valence-like non-singlet up-sector |PID|.
-pub const PID_NSM_U: u16 = 10202;
-
-/// valence-like non-singlet down-sector |PID|.
-pub const PID_NSM_D: u16 = 10203;
+pub(super) use math_constants::*;
+pub use pid::*;
+pub(super) use qcd_constants::*;
 
 /// Maximum QCD coupling power implemented.
 pub const MAX_ORDER_QCD: usize = 4;
 
 /// Maximum QED coupling power implemented.
 pub const MAX_ORDER_QED: usize = 2;
-
-/// |QED| electric charge combinations.
-pub struct ChargeCombinations {
-    pub nf: u8,
-}
-
-impl ChargeCombinations {
-    /// Number of up-like flavors.
-    pub fn nu(&self) -> u8 {
-        self.nf / 2
-    }
-
-    /// Number of down-like flavors.
-    pub fn nd(&self) -> u8 {
-        self.nf - self.nu()
-    }
-
-    /// Charge average.
-    pub fn e2avg(&self) -> f64 {
-        (self.nu() as f64 * EU2 + self.nd() as f64 * ED2) / (self.nf as f64)
-    }
-
-    /// Relative up contribution.
-    pub fn vu(&self) -> f64 {
-        self.nu() as f64 / (self.nf as f64)
-    }
-
-    /// Relative down contribution.
-    pub fn vd(&self) -> f64 {
-        self.nd() as f64 / (self.nf as f64)
-    }
-
-    /// Relative up contribution to charge difference.
-    pub fn vue2m(&self) -> f64 {
-        self.vu() * (EU2 - ED2)
-    }
-
-    /// Relative down contribution to charge difference.
-    pub fn vde2m(&self) -> f64 {
-        self.vd() * (EU2 - ED2)
-    }
-
-    /// Asymmetric charge combination.
-    pub fn e2delta(&self) -> f64 {
-        self.vde2m() - self.vue2m() + self.e2avg()
-    }
-}
diff --git a/crates/ekore/src/harmonics.rs b/crates/ekore/src/harmonics.rs
index 75dcfd41b..54850550e 100644
--- a/crates/ekore/src/harmonics.rs
+++ b/crates/ekore/src/harmonics.rs
@@ -1,11 +1,11 @@
 //! Tools to compute harmonic sums and related special functions.
 
 pub mod cache;
-pub mod g_functions;
-pub mod log_functions;
-pub mod polygamma;
-pub mod w1;
-pub mod w2;
-pub mod w3;
-pub mod w4;
-pub mod w5;
+mod g_functions;
+pub(super) mod log_functions;
+mod polygamma;
+mod w1;
+mod w2;
+mod w3;
+mod w4;
+mod w5;
diff --git a/crates/ekore/src/harmonics/cache.rs b/crates/ekore/src/harmonics/cache.rs
index 36bdb1ff1..7d6366aff 100644
--- a/crates/ekore/src/harmonics/cache.rs
+++ b/crates/ekore/src/harmonics/cache.rs
@@ -6,7 +6,7 @@ use crate::harmonics::{g_functions, w1, w2, w3, w4, w5};
 
 /// List of available elements.
 #[derive(Debug, PartialEq, Eq)]
-pub enum K {
+pub(crate) enum K {
     /// $S_1(N)$
     S1,
     /// $S_2(N)$
@@ -101,7 +101,7 @@ impl Cache {
     }
 
     /// Retrieve an element.
-    pub fn get(&mut self, k: K) -> Complex<f64> {
+    pub(crate) fn get(&mut self, k: K) -> Complex<f64> {
         let idx = k.idx();
         if let Some(val) = self.m[idx] {
             return val;
@@ -138,7 +138,7 @@ impl Cache {
 /// 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(
+pub(crate) fn recursive_harmonic_sum(
     base_value: Complex<f64>,
     n: Complex<f64>,
     iterations: usize,
diff --git a/crates/ekore/src/harmonics/g_functions.rs b/crates/ekore/src/harmonics/g_functions.rs
index 2cec4b7f3..55eb245a2 100644
--- a/crates/ekore/src/harmonics/g_functions.rs
+++ b/crates/ekore/src/harmonics/g_functions.rs
@@ -12,7 +12,7 @@ use num::{Zero, complex::Complex};
 ///
 /// 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(c: &mut Cache) -> Complex<f64> {
+pub(super) fn g3(c: &mut Cache) -> Complex<f64> {
     let N = c.n();
     let S1 = c.get(K::S1);
     const CS: [f64; 7] = [
diff --git a/crates/ekore/src/harmonics/log_functions.rs b/crates/ekore/src/harmonics/log_functions.rs
index 1db2bcf9f..05d9b0811 100644
--- a/crates/ekore/src/harmonics/log_functions.rs
+++ b/crates/ekore/src/harmonics/log_functions.rs
@@ -10,14 +10,14 @@ use num::complex::Complex;
 use num::pow;
 
 /// Mellin transform of $(1-x)\ln(1-x)$.
-pub fn lm11m1(c: &mut Cache) -> Complex<f64> {
+pub(crate) fn lm11m1(c: &mut Cache) -> Complex<f64> {
     let n = c.n();
     let S1 = c.get(K::S1);
     1. / (1. + n).powu(2) - S1 / (1. + n).powu(2) - S1 / (n * (1. + n).powu(2))
 }
 
 /// Mellin transform of $(1-x)\ln^2(1-x)$.
-pub fn lm12m1(c: &mut Cache) -> Complex<f64> {
+pub(crate) fn lm12m1(c: &mut Cache) -> Complex<f64> {
     let n = c.n();
     let S1 = c.get(K::S1);
     let S2 = c.get(K::S2);
@@ -27,7 +27,7 @@ pub fn lm12m1(c: &mut Cache) -> Complex<f64> {
 }
 
 /// Mellin transform of $(1-x)\ln^3(1-x)$.
-pub fn lm13m1(c: &mut Cache) -> Complex<f64> {
+pub(crate) fn lm13m1(c: &mut Cache) -> Complex<f64> {
     let n = c.n();
     let S1 = c.get(K::S1);
     let S2 = c.get(K::S2);
@@ -40,7 +40,7 @@ pub fn lm13m1(c: &mut Cache) -> Complex<f64> {
 }
 
 /// Mellin transform of $(1-x)\ln^4(1-x)$.
-pub fn lm14m1(c: &mut Cache) -> Complex<f64> {
+pub(crate) fn lm14m1(c: &mut Cache) -> Complex<f64> {
     let n = c.n();
     let S1 = c.get(K::S1);
     let S2 = c.get(K::S2);
@@ -66,7 +66,7 @@ pub fn lm14m1(c: &mut Cache) -> Complex<f64> {
 }
 
 /// Mellin transform of $(1-x)\ln^5(1-x)$.
-pub fn lm15m1(c: &mut Cache) -> Complex<f64> {
+pub(crate) fn lm15m1(c: &mut Cache) -> Complex<f64> {
     let n = c.n();
     let S1 = c.get(K::S1);
     let S2 = c.get(K::S2);
@@ -93,14 +93,14 @@ pub fn lm15m1(c: &mut Cache) -> Complex<f64> {
 }
 
 /// Mellin transform of $\ln(1-x)$.
-pub fn lm11(c: &mut Cache) -> Complex<f64> {
+pub(crate) fn lm11(c: &mut Cache) -> Complex<f64> {
     let n = c.n();
     let S1 = c.get(K::S1);
     -S1 / n
 }
 
 /// Mellin transform of $\ln^2(1-x)$.
-pub fn lm12(c: &mut Cache) -> Complex<f64> {
+pub(crate) fn lm12(c: &mut Cache) -> Complex<f64> {
     let n = c.n();
     let S1 = c.get(K::S1);
     let S2 = c.get(K::S2);
@@ -108,7 +108,7 @@ pub fn lm12(c: &mut Cache) -> Complex<f64> {
 }
 
 /// Mellin transform of $\ln^3(1-x)$.
-pub fn lm13(c: &mut Cache) -> Complex<f64> {
+pub(crate) fn lm13(c: &mut Cache) -> Complex<f64> {
     let n = c.n();
     let S1 = c.get(K::S1);
     let S2 = c.get(K::S2);
@@ -117,7 +117,7 @@ pub fn lm13(c: &mut Cache) -> Complex<f64> {
 }
 
 /// Mellin transform of $\ln^4(1-x)$.
-pub fn lm14(c: &mut Cache) -> Complex<f64> {
+pub(crate) fn lm14(c: &mut Cache) -> Complex<f64> {
     let n = c.n();
     let S1 = c.get(K::S1);
     let S2 = c.get(K::S2);
@@ -127,7 +127,7 @@ pub fn lm14(c: &mut Cache) -> Complex<f64> {
 }
 
 /// Mellin transform of $\ln^5(1-x)$.
-pub fn lm15(c: &mut Cache) -> Complex<f64> {
+pub(crate) fn lm15(c: &mut Cache) -> Complex<f64> {
     let n = c.n();
     let S1 = c.get(K::S1);
     let S2 = c.get(K::S2);
@@ -143,14 +143,14 @@ pub fn lm15(c: &mut Cache) -> Complex<f64> {
 }
 
 /// Mellin transform of $(1-x)^2\ln(1-x)$.
-pub fn lm11m2(c: &mut Cache) -> Complex<f64> {
+pub(crate) fn lm11m2(c: &mut Cache) -> Complex<f64> {
     let n = c.n();
     let S1 = c.get(K::S1);
     (5. + 3. * n - (2. * (1. + n) * (2. + n) * S1) / n) / ((1. + n).powu(2) * (2. + n).powu(2))
 }
 
 /// Mellin transform of $(1-x)^2\ln^2(1-x)$.
-pub fn lm12m2(c: &mut Cache) -> Complex<f64> {
+pub(crate) fn lm12m2(c: &mut Cache) -> Complex<f64> {
     let n = c.n();
     let S1 = c.get(K::S1);
     let S2 = c.get(K::S2);
@@ -162,7 +162,7 @@ pub fn lm12m2(c: &mut Cache) -> Complex<f64> {
 }
 
 /// Mellin transform of $(1-x)^2\ln^3(1-x)$.
-pub fn lm13m2(c: &mut Cache) -> Complex<f64> {
+pub(crate) fn lm13m2(c: &mut Cache) -> Complex<f64> {
     let n = c.n();
     let S1 = c.get(K::S1);
     let S2 = c.get(K::S2);
@@ -179,7 +179,7 @@ pub fn lm13m2(c: &mut Cache) -> Complex<f64> {
 }
 
 /// Mellin transform of $(1-x)^2\ln^4(1-x)$.
-pub fn lm14m2(c: &mut Cache) -> Complex<f64> {
+pub(crate) fn lm14m2(c: &mut Cache) -> Complex<f64> {
     let n = c.n();
     let S1 = c.get(K::S1);
     let S2 = c.get(K::S2);
diff --git a/crates/ekore/src/harmonics/polygamma.rs b/crates/ekore/src/harmonics/polygamma.rs
index 9263bcccf..77bd41d4a 100644
--- a/crates/ekore/src/harmonics/polygamma.rs
+++ b/crates/ekore/src/harmonics/polygamma.rs
@@ -10,7 +10,7 @@ use std::f64::consts::PI;
 /// given by [\[KOLBIG1972221\]](crate::bib::KOLBIG1972221).
 ///
 /// TODO: introduce back errors
-pub fn cern_polygamma(Z: Complex<f64>, K: usize) -> Complex<f64> {
+pub(super) fn cern_polygamma(Z: Complex<f64>, K: usize) -> Complex<f64> {
     // const DELTA: f64 = 5e-13;
     const R1: f64 = 1.0;
     const HF: f64 = R1 / 2.0;
diff --git a/crates/ekore/src/harmonics/w1.rs b/crates/ekore/src/harmonics/w1.rs
index 2421d12d8..092b1163f 100644
--- a/crates/ekore/src/harmonics/w1.rs
+++ b/crates/ekore/src/harmonics/w1.rs
@@ -8,19 +8,19 @@ use crate::harmonics::polygamma::cern_polygamma;
 ///
 /// $$S_1(N) = \sum\limits_{j=1}^N \frac 1 j = \psi_0(N+1)+\gamma_E$$
 /// with $\psi_0(N)$ the digamma function and $\gamma_E$ the Euler-Mascheroni constant.
-pub fn S1(N: Complex<f64>) -> Complex<f64> {
+pub(super) fn S1(N: Complex<f64>) -> Complex<f64> {
     cern_polygamma(N + 1.0, 0) + 0.577_215_664_901_532_9
 }
 
 /// Analytic continuation of harmonic sum $S_{-1}(N)$ for even moments.
-pub fn Sm1e(c: &mut Cache) -> Complex<f64> {
+pub(super) fn Sm1e(c: &mut Cache) -> Complex<f64> {
     let hS1 = c.get(K::S1);
     let hS1h = c.get(K::S1h);
     hS1h - hS1
 }
 
 /// Analytic continuation of harmonic sum $S_{-1}(N)$ for odd moments.
-pub fn Sm1o(c: &mut Cache) -> Complex<f64> {
+pub(super) fn Sm1o(c: &mut Cache) -> Complex<f64> {
     let hS1 = c.get(K::S1);
     let hS1mh = c.get(K::S1mh);
     hS1mh - hS1
diff --git a/crates/ekore/src/harmonics/w2.rs b/crates/ekore/src/harmonics/w2.rs
index 643ed19ed..4d3c0672f 100644
--- a/crates/ekore/src/harmonics/w2.rs
+++ b/crates/ekore/src/harmonics/w2.rs
@@ -9,19 +9,19 @@ use crate::harmonics::polygamma::cern_polygamma;
 ///
 /// $$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> {
+pub(super) fn S2(N: Complex<f64>) -> Complex<f64> {
     -cern_polygamma(N + 1.0, 1) + ZETA2
 }
 
 /// Analytic continuation of harmonic sum $S_{-2}(N)$ for even moments.
-pub fn Sm2e(c: &mut Cache) -> Complex<f64> {
+pub(super) fn Sm2e(c: &mut Cache) -> Complex<f64> {
     let hS2 = c.get(K::S2);
     let hS2h = c.get(K::S2h);
     1. / 2. * hS2h - hS2
 }
 
 /// Analytic continuation of harmonic sum $S_{-2}(N)$ for odd moments.
-pub fn Sm2o(c: &mut Cache) -> Complex<f64> {
+pub(super) fn Sm2o(c: &mut Cache) -> Complex<f64> {
     let hS2 = c.get(K::S2);
     let hS2mh = c.get(K::S2mh);
     1. / 2. * hS2mh - hS2
diff --git a/crates/ekore/src/harmonics/w3.rs b/crates/ekore/src/harmonics/w3.rs
index 81bc0317c..76680624a 100644
--- a/crates/ekore/src/harmonics/w3.rs
+++ b/crates/ekore/src/harmonics/w3.rs
@@ -12,33 +12,33 @@ use crate::harmonics::polygamma::cern_polygamma;
 ///
 /// $$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> {
+pub(super) fn S3(N: Complex<f64>) -> Complex<f64> {
     0.5 * cern_polygamma(N + 1.0, 2) + ZETA3
 }
 
 /// Analytic continuation of harmonic sum $S_{-3}(N)$ for even moments.
-pub fn Sm3e(c: &mut Cache) -> Complex<f64> {
+pub(super) fn Sm3e(c: &mut Cache) -> Complex<f64> {
     let hS3 = c.get(K::S3);
     let hS3h = c.get(K::S3h);
     1. / (2.).pow(2) * hS3h - hS3
 }
 
 /// Analytic continuation of harmonic sum $S_{-3}(N)$ for odd moments.
-pub fn Sm3o(c: &mut Cache) -> Complex<f64> {
+pub(super) fn Sm3o(c: &mut Cache) -> Complex<f64> {
     let hS3 = c.get(K::S3);
     let hS3mh = c.get(K::S3mh);
     1. / (2.).pow(2) * hS3mh - hS3
 }
 
 /// Analytic continuation of harmonic sum $S_{-2,1}(N)$ for even moments.
-pub fn Sm21e(c: &mut Cache) -> Complex<f64> {
+pub(super) fn Sm21e(c: &mut Cache) -> Complex<f64> {
     let hSm1 = c.get(K::Sm1e);
     let mut cp1 = Cache::new(c.n() + 1.);
     ZETA2 * hSm1 - 5. / 8. * ZETA3 + ZETA2 * LN_2 - g3(&mut cp1)
 }
 
 /// Analytic continuation of harmonic sum $S_{-2,1}(N)$ for odd moments.
-pub fn Sm21o(c: &mut Cache) -> Complex<f64> {
+pub(super) fn Sm21o(c: &mut Cache) -> Complex<f64> {
     let hSm1 = c.get(K::Sm1o);
     let mut cp1 = Cache::new(c.n() + 1.);
     ZETA2 * hSm1 - 5. / 8. * ZETA3 + ZETA2 * LN_2 + g3(&mut cp1)
diff --git a/crates/ekore/src/harmonics/w4.rs b/crates/ekore/src/harmonics/w4.rs
index 66835ac5f..1bcb351b6 100644
--- a/crates/ekore/src/harmonics/w4.rs
+++ b/crates/ekore/src/harmonics/w4.rs
@@ -8,6 +8,6 @@ use crate::harmonics::polygamma::cern_polygamma;
 ///
 /// $$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> {
+pub(super) fn S4(N: Complex<f64>) -> Complex<f64> {
     ZETA4 - 1.0 / 6.0 * cern_polygamma(N + 1.0, 3)
 }
diff --git a/crates/ekore/src/harmonics/w5.rs b/crates/ekore/src/harmonics/w5.rs
index e77730a5b..6978dff3b 100644
--- a/crates/ekore/src/harmonics/w5.rs
+++ b/crates/ekore/src/harmonics/w5.rs
@@ -8,6 +8,6 @@ use crate::harmonics::polygamma::cern_polygamma;
 ///
 /// $$S_5(N) = \sum\limits_{j=1}^N \frac 1 {j^5} = - \frac 1 24 \psi_4(N+1)+\zeta(5)$$
 /// with $\psi_4(N)$ the 4rd polygamma function and $\zeta$ the Riemann zeta function.
-pub fn S5(N: Complex<f64>) -> Complex<f64> {
+pub(super) fn S5(N: Complex<f64>) -> Complex<f64> {
     ZETA5 + 1.0 / 24.0 * cern_polygamma(N + 1.0, 4)
 }
diff --git a/crates/ekore/src/lib.rs b/crates/ekore/src/lib.rs
index 6505e6957..3ca2ed118 100644
--- a/crates/ekore/src/lib.rs
+++ b/crates/ekore/src/lib.rs
@@ -8,4 +8,4 @@ pub mod bib;
 pub mod constants;
 pub mod harmonics;
 pub mod operator_matrix_elements;
-pub mod util;
+mod util;
diff --git a/crates/ekore/src/math_constants.rs b/crates/ekore/src/math_constants.rs
new file mode 100644
index 000000000..ec42c6ea0
--- /dev/null
+++ b/crates/ekore/src/math_constants.rs
@@ -0,0 +1,17 @@
+//! Math Constants.
+
+/// Riemann zeta function at z = 2.
+///
+/// $\zeta(2) = \pi^2 / 6$.
+pub(crate) const ZETA2: f64 = 1.6449340668482264;
+
+/// Riemann zeta function at z = 3.
+pub(crate) const ZETA3: f64 = 1.2020569031595942;
+
+/// Riemann zeta function at z = 4.
+///
+/// $\zeta(4) = \pi^4 / 90$.
+pub(crate) const ZETA4: f64 = 1.082323233711138;
+
+/// Riemann zeta function at z = 5.
+pub(crate) const ZETA5: f64 = 1.03692775514337;
diff --git a/crates/ekore/src/operator_matrix_elements/unpolarized/spacelike.rs b/crates/ekore/src/operator_matrix_elements/unpolarized/spacelike.rs
index b20125ddd..7dee10e7d 100644
--- a/crates/ekore/src/operator_matrix_elements/unpolarized/spacelike.rs
+++ b/crates/ekore/src/operator_matrix_elements/unpolarized/spacelike.rs
@@ -3,8 +3,8 @@ use crate::constants::MAX_ORDER_QCD;
 use crate::harmonics::cache::Cache;
 use num::Zero;
 use num::complex::Complex;
-pub mod as1;
-pub mod as2;
+mod as1;
+mod as2;
 
 /// Compute the tower of the singlet |OME|.
 ///
diff --git a/crates/ekore/src/operator_matrix_elements/unpolarized/spacelike/as1.rs b/crates/ekore/src/operator_matrix_elements/unpolarized/spacelike/as1.rs
index e3b1d5a6b..b9f5331c4 100644
--- a/crates/ekore/src/operator_matrix_elements/unpolarized/spacelike/as1.rs
+++ b/crates/ekore/src/operator_matrix_elements/unpolarized/spacelike/as1.rs
@@ -10,7 +10,7 @@ use crate::harmonics::cache::{Cache, K};
 /// Compute heavy-heavy |OME|.
 ///
 /// Implements Eq. (20a) of [\[Ball:2015tna\]](crate::bib::Ball2015tna).
-pub fn A_hh(c: &mut Cache, _nf: u8, L: f64) -> Complex<f64> {
+pub(super) fn A_hh(c: &mut Cache, _nf: u8, L: f64) -> Complex<f64> {
     let N = c.n();
     let S1m = c.get(K::S1) - 1. / N;
     let S2m = c.get(K::S2) - 1. / N.powu(2);
@@ -29,7 +29,7 @@ pub fn A_hh(c: &mut Cache, _nf: u8, L: f64) -> Complex<f64> {
 /// Compute gluon-heavy |OME|.
 ///
 /// Implements Eq. (20b) of [\[Ball:2015tna\]](crate::bib::Ball2015tna).
-pub fn A_gh(c: &mut Cache, _nf: u8, L: f64) -> Complex<f64> {
+pub(super) fn A_gh(c: &mut Cache, _nf: u8, L: f64) -> Complex<f64> {
     let N = c.n();
     let agh_l1 = (2. + N + N.powu(2)) / (N * (N.powu(2) - 1.));
     let agh_l0 = (-4. + 2. * N + N.powu(2) * (15. + N * (3. + N - N.powu(2))))
@@ -40,7 +40,7 @@ pub fn A_gh(c: &mut Cache, _nf: u8, L: f64) -> Complex<f64> {
 /// Compute heavy-gluon |OME|.
 ///
 /// Implements Eq. (B.2) of [\[Buza:1996wv\]](crate::bib::Buza1996wv).
-pub fn A_hg(c: &mut Cache, _nf: u8, L: f64) -> Complex<f64> {
+pub(super) fn A_hg(c: &mut Cache, _nf: u8, L: f64) -> Complex<f64> {
     let N = c.n();
     let den = 1. / (N * (N + 1.) * (2. + N));
     let num = 2. * (2. + N + N.powu(2));
@@ -50,12 +50,12 @@ pub fn A_hg(c: &mut Cache, _nf: u8, L: f64) -> Complex<f64> {
 /// Compute gluon-gluon |OME|.
 ///
 /// Implements Eq. (B.6) of [\[Buza:1996wv\]](crate::bib::Buza1996wv).
-pub fn A_gg(_c: &mut Cache, _nf: u8, L: f64) -> Complex<f64> {
+pub(super) fn A_gg(_c: &mut Cache, _nf: u8, L: f64) -> Complex<f64> {
     cmplx!(-2.0 / 3.0 * L, 0.)
 }
 
 /// Compute the |NLO| singlet |OME|.
-pub fn A_singlet(c: &mut Cache, nf: u8, L: f64) -> [[Complex<f64>; 3]; 3] {
+pub(super) fn A_singlet(c: &mut Cache, nf: u8, L: f64) -> [[Complex<f64>; 3]; 3] {
     [
         [A_gg(c, nf, L), Complex::<f64>::zero(), A_gh(c, nf, L)],
         [
@@ -68,7 +68,7 @@ pub fn A_singlet(c: &mut Cache, nf: u8, L: f64) -> [[Complex<f64>; 3]; 3] {
 }
 
 /// Compute the |NLO| non-singlet |OME|.
-pub fn A_ns(c: &mut Cache, nf: u8, L: f64) -> [[Complex<f64>; 2]; 2] {
+pub(super) fn A_ns(c: &mut Cache, nf: u8, L: f64) -> [[Complex<f64>; 2]; 2] {
     [
         [Complex::<f64>::zero(), Complex::<f64>::zero()],
         [Complex::<f64>::zero(), A_hh(c, nf, L)],
diff --git a/crates/ekore/src/operator_matrix_elements/unpolarized/spacelike/as2.rs b/crates/ekore/src/operator_matrix_elements/unpolarized/spacelike/as2.rs
index d82fcd89c..3d1e8bc88 100644
--- a/crates/ekore/src/operator_matrix_elements/unpolarized/spacelike/as2.rs
+++ b/crates/ekore/src/operator_matrix_elements/unpolarized/spacelike/as2.rs
@@ -12,7 +12,7 @@ use crate::operator_matrix_elements::unpolarized::spacelike::as1;
 /// |NNLO| light-light non-singlet |OME|.
 ///
 /// Implements Eq. (B.4) of [\[Buza:1996wv\]](crate::bib::Buza1996wv).
-pub fn A_qq_ns(c: &mut Cache, _nf: u8, L: f64) -> Complex<f64> {
+pub(super) fn A_qq_ns(c: &mut Cache, _nf: u8, L: f64) -> Complex<f64> {
     let N = c.n();
     let S1 = c.get(K::S1);
     let S2 = c.get(K::S2);
@@ -44,7 +44,7 @@ pub fn A_qq_ns(c: &mut Cache, _nf: u8, L: f64) -> Complex<f64> {
 /// |NNLO| heavy-light pure-singlet |OME|.
 ///
 /// Implements Eq. (B.1) of [\[Buza:1996wv\]](crate::bib::Buza1996wv).
-pub fn A_hq_ps(c: &mut Cache, _nf: u8, L: f64) -> Complex<f64> {
+pub(super) fn A_hq_ps(c: &mut Cache, _nf: u8, L: f64) -> Complex<f64> {
     let N = c.n();
     let S2 = c.get(K::S2);
     let F1M = 1.0 / (N - 1.0) * (ZETA2 - (S2 - 1.0 / N.powu(2)));
@@ -85,7 +85,7 @@ pub fn A_hq_ps(c: &mut Cache, _nf: u8, L: f64) -> Complex<f64> {
 ///
 /// Implements Eq. (B.3) of [\[Buza:1996wv\]](crate::bib::Buza1996wv).
 /// The expession for ``A_Hg_l0`` comes form [\[Bierenbaum:2009zt\]](crate::bib::Bierenbaum2009zt).
-pub fn A_hg(c: &mut Cache, _nf: u8, L: f64) -> Complex<f64> {
+pub(super) fn A_hg(c: &mut Cache, _nf: u8, L: f64) -> Complex<f64> {
     let N = c.n();
     let S1 = c.get(K::S1);
     let S2 = c.get(K::S2);
@@ -185,7 +185,7 @@ pub fn A_hg(c: &mut Cache, _nf: u8, L: f64) -> Complex<f64> {
 /// |NNLO| gluon-quark |OME|.
 ///
 /// Implements Eq. (B.5) of [\[Buza:1996wv\]](crate::bib::Buza1996wv).
-pub fn A_gq(c: &mut Cache, _nf: u8, L: f64) -> Complex<f64> {
+pub(super) fn A_gq(c: &mut Cache, _nf: u8, L: f64) -> Complex<f64> {
     let N = c.n();
     let S1 = c.get(K::S1);
     let S2 = c.get(K::S2);
@@ -212,7 +212,7 @@ pub fn A_gq(c: &mut Cache, _nf: u8, L: f64) -> Complex<f64> {
 /// |NNLO| gluon-gluon |OME|.
 ///
 /// Implements Eq. (B.7) of [\[Buza:1996wv\]](crate::bib::Buza1996wv).
-pub fn A_gg(c: &mut Cache, _nf: u8, L: f64) -> Complex<f64> {
+pub(super) fn A_gg(c: &mut Cache, _nf: u8, L: f64) -> Complex<f64> {
     let N = c.n();
     let S1 = c.get(K::S1);
     let S1m = S1 - 1. / N;
@@ -260,7 +260,12 @@ pub fn A_gg(c: &mut Cache, _nf: u8, L: f64) -> Complex<f64> {
 }
 
 /// |NNLO| singlet |OME|.
-pub fn A_singlet(c: &mut Cache, nf: u8, L: f64, is_msbar_mass: bool) -> [[Complex<f64>; 3]; 3] {
+pub(super) fn A_singlet(
+    c: &mut Cache,
+    nf: u8,
+    L: f64,
+    is_msbar_mass: bool,
+) -> [[Complex<f64>; 3]; 3] {
     let A_hq_2 = A_hq_ps(c, nf, L);
     let A_qq_2 = A_qq_ns(c, nf, L);
     let mut A_hg_2 = A_hg(c, nf, L);
@@ -280,7 +285,7 @@ pub fn A_singlet(c: &mut Cache, nf: u8, L: f64, is_msbar_mass: bool) -> [[Comple
 }
 
 /// |NNLO| non-singlet |OME|.
-pub fn A_ns(c: &mut Cache, nf: u8, L: f64) -> [[Complex<f64>; 2]; 2] {
+pub(super) fn A_ns(c: &mut Cache, nf: u8, L: f64) -> [[Complex<f64>; 2]; 2] {
     [
         [A_qq_ns(c, nf, L), Complex::<f64>::zero()],
         [Complex::<f64>::zero(), Complex::<f64>::zero()],
diff --git a/crates/ekore/src/pid.rs b/crates/ekore/src/pid.rs
new file mode 100644
index 000000000..dfed0e2b9
--- /dev/null
+++ b/crates/ekore/src/pid.rs
@@ -0,0 +1,79 @@
+//! PID Constants
+
+/// Up quark electric charge square.
+///
+/// Defaults to $e_u^2 = 4./9$.
+pub(crate) const EU2: f64 = 4. / 9.;
+
+/// Down quark electric charge square.
+///
+/// Defaults to $e_d^2 = 1./9$.
+pub(crate) const ED2: f64 = 1. / 9.;
+
+/// 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;
+
+/// singlet-like non-singlet up-sector |PID|.
+pub const PID_NSP_U: u16 = 10102;
+
+/// singlet-like non-singlet down-sector |PID|.
+pub const PID_NSP_D: u16 = 10103;
+
+/// valence-like non-singlet up-sector |PID|.
+pub const PID_NSM_U: u16 = 10202;
+
+/// valence-like non-singlet down-sector |PID|.
+pub const PID_NSM_D: u16 = 10203;
+
+/// |QED| electric charge combinations.
+pub(crate) struct ChargeCombinations {
+    pub nf: u8,
+}
+
+impl ChargeCombinations {
+    /// Number of up-like flavors.
+    pub(crate) fn nu(&self) -> u8 {
+        self.nf / 2
+    }
+
+    /// Number of down-like flavors.
+    pub(crate) fn nd(&self) -> u8 {
+        self.nf - self.nu()
+    }
+
+    /// Charge average.
+    pub(crate) fn e2avg(&self) -> f64 {
+        (self.nu() as f64 * EU2 + self.nd() as f64 * ED2) / (self.nf as f64)
+    }
+
+    /// Relative up contribution.
+    pub(crate) fn vu(&self) -> f64 {
+        self.nu() as f64 / (self.nf as f64)
+    }
+
+    /// Relative down contribution.
+    pub(crate) fn vd(&self) -> f64 {
+        self.nd() as f64 / (self.nf as f64)
+    }
+
+    /// Relative up contribution to charge difference.
+    pub(crate) fn vue2m(&self) -> f64 {
+        self.vu() * (EU2 - ED2)
+    }
+
+    /// Relative down contribution to charge difference.
+    pub(crate) fn vde2m(&self) -> f64 {
+        self.vd() * (EU2 - ED2)
+    }
+
+    /// Asymmetric charge combination.
+    pub(crate) fn e2delta(&self) -> f64 {
+        self.vde2m() - self.vue2m() + self.e2avg()
+    }
+}
diff --git a/crates/ekore/src/qcd_constants.rs b/crates/ekore/src/qcd_constants.rs
new file mode 100644
index 000000000..8c757dc65
--- /dev/null
+++ b/crates/ekore/src/qcd_constants.rs
@@ -0,0 +1,21 @@
+//! QCD Constants.
+
+/// The number of colors.
+///
+/// Defaults to $N_C = 3$.
+pub(crate) const NC: u8 = 3;
+
+/// The normalization of fundamental generators.
+///
+/// Defaults to $T_R = 1/2$.
+pub(crate) const TR: f64 = 1.0 / 2.0;
+
+/// Second Casimir constant in the adjoint representation.
+///
+/// Defaults to $C_A = N_C = 3$.
+pub(crate) const CA: f64 = NC as f64;
+
+/// Second Casimir constant in the fundamental representation.
+///
+/// Defaults to $C_F = \frac{N_C^2-1}{2N_C} = 4/3$.
+pub(crate) const CF: f64 = ((NC * NC - 1) as f64) / ((2 * NC) as f64);