diff --git a/benchmarks/lha_paper_bench.py b/benchmarks/lha_paper_bench.py
index 806456ea6..65bc64fb8 100644
--- a/benchmarks/lha_paper_bench.py
+++ b/benchmarks/lha_paper_bench.py
@@ -1,6 +1,6 @@
 # -*- coding: utf-8 -*-
 """
-Benchmark to :cite:`Giele:2002hx` (LO + NLO) and :cite:`Dittmar:2005ed` (NNLO)
+Benchmark to :cite:`Giele:2002hx` (LO + NLO) and :cite:`Dittmar:2005ed` (NNLO).
 """
 import numpy as np
 from banana import register
@@ -35,9 +35,7 @@
 
 
 class LHABenchmark(Runner):
-    """
-    Globally set the external program to LHA
-    """
+    """Globally set the external program to LHA."""
 
     def __init__(self):
         super().__init__()
@@ -46,36 +44,34 @@ def __init__(self):
         self.rotate_to_evolution_basis = True
 
     def plain_theory(self, pto):
-        """
-        Plain theories at given PTO.
+        """Generate plain theories at given PTO.
 
         Parameters
         ----------
-            pto : int
-                perturbation order
+        pto : int
+            perturbation order
 
         Returns
         -------
-            list(dict)
-                theory updates
+        list(dict)
+            theory updates
         """
         th = self.theory.copy()
         th.update({"PTO": pto})
         return [th]
 
     def sv_theories(self, pto):
-        """
-        Scale variation theories.
+        """Generate scale variation theories.
 
         Parameters
         ----------
-            pto : int
-                perturbation order
+        pto : int
+            perturbation order
 
         Returns
         -------
-            list(dict)
-                theory updates
+        list(dict)
+            theory updates
         """
         low = self.theory.copy()
         low["PTO"] = pto
@@ -89,29 +85,27 @@ def sv_theories(self, pto):
 
     @staticmethod
     def skip_pdfs(_theory):
-        """
-        Adjust skip_pdf by the used theory
+        """Adjust skip_pdf by the used theory.
 
         Parameters
         ----------
-            theory_updates : list(dict)
-                theory updates
+        theory : dict
+            theory update
 
         Returns
         -------
-            list :
-                current skip_pdf
+        list
+            allowed PDFs in LHA
         """
         return default_skip_pdfs
 
     def run_lha(self, theory_updates):
-        """
-        Enforce operators and PDF
+        """Enforce operator grid and PDF.
 
         Parameters
         ----------
-            theory_updates : list(dict)
-                theory updates
+        theory_updates : list(dict)
+            theory updates
         """
         self.run(
             theory_updates,
@@ -126,16 +120,16 @@ def run_lha(self, theory_updates):
         )
 
     def benchmark_plain(self, pto):
-        """Plain configuration"""
+        """Run plain configuration."""
         self.run_lha(self.plain_theory(pto))
 
     def benchmark_sv(self, pto):
-        """Scale variations"""
+        """Run scale variations."""
         self.run_lha(self.sv_theories(pto))
 
 
 class BenchmarkVFNS(LHABenchmark):
-    """Variable Flavor Number Scheme"""
+    """Provide |VFNS| settings."""
 
     def __init__(self):
         super().__init__()
@@ -153,7 +147,7 @@ def __init__(self):
 
 
 class BenchmarkFFNS(LHABenchmark):
-    """Fixed Flavor Number Scheme"""
+    """Provide |FFNS| settings."""
 
     def __init__(self):
         super().__init__()
@@ -171,19 +165,29 @@ def __init__(self):
 
     @staticmethod
     def skip_pdfs(theory):
+        """Adjust skip_pdf by the used theory.
+
+        Parameters
+        ----------
+        theory : dict
+            theory update
+
+        Returns
+        -------
+        list
+            allowed PDFs in FFNS LHA
+        """
         ffns_skip_pdfs = default_skip_pdfs.copy()
         # remove bottom
         ffns_skip_pdfs.extend([-5, 5, "T24"])
-        # in NNLO V8 becomes available
+        # in NNLO also V8 gets removed
         if theory["PTO"] >= 2:
             ffns_skip_pdfs.remove("V8")
         return ffns_skip_pdfs
 
 
 class BenchmarkRunner(BenchmarkVFNS):
-    """
-    Generic benchmark runner using the LHA VFNS settings
-    """
+    """Generic benchmark runner using the LHA |VFNS| settings."""
 
     def __init__(self, external):
         super().__init__()
@@ -191,18 +195,12 @@ def __init__(self, external):
         self.sandbox = True
 
     def benchmark_sv(self, pto):
-        """
-        Scale variations
+        """Run scale variations.
 
         Parameters
         ----------
-            pto : int
-                perturbation order
-
-        Returns
-        -------
-            list(dict)
-                theory updates
+        pto : int
+            perturbation order
         """
         high, low = self.sv_theories(pto)
 
@@ -223,15 +221,15 @@ def benchmark_sv(self, pto):
 if __name__ == "__main__":
 
     # Benchmark to LHA
-    # obj = BenchmarkVFNS()
-    obj = BenchmarkFFNS()
+    obj = BenchmarkVFNS()
+    # obj = BenchmarkFFNS()
 
     # obj.benchmark_plain(0)
-    # obj.benchmark_sv(2)
+    obj.benchmark_sv(2)
 
     # # VFNS benchmarks with LHA settings
-    programs = ["LHA", "pegasus", "apfel"]
-    for p in programs:
-        obj = BenchmarkRunner(p)
-        # obj.benchmark_plain(2)
-        obj.benchmark_sv(2)
+    # programs = ["LHA", "pegasus", "apfel"]
+    # for p in programs:
+    #     obj = BenchmarkRunner(p)
+    #     # obj.benchmark_plain(2)
+    #     obj.benchmark_sv(2)
diff --git a/doc/source/theory/pQCD.rst b/doc/source/theory/pQCD.rst
index 3f50a7a4c..ad8422671 100644
--- a/doc/source/theory/pQCD.rst
+++ b/doc/source/theory/pQCD.rst
@@ -156,11 +156,12 @@ This method provides many advantages:
     * it naturally incorporates renormalization group invariance,
       as the perturbative order increases, estimates of |MHOU| decrease;
     * the same procedure can be used for any perturbative process,
-      since the scale dependence of the strong coupling :math:`a_s(\mu^2)` and of PDFs is universal;
+      since the scale dependence of the strong coupling :math:`a_s(\mu^2)` and of |PDF| are universal;
 
 However, there is no unique prescription to determine the specific range of the scale variation,
-the most used prescription specify to vary the factor :math:`\mu_F/\mu_R` in the range:
-:math:`1/2 \le \mu_F/\mu_R \le 2`.
+the most common prescription specify to vary the factor :math:`\mu_F/\mu_R` in the range:
+:math:`1/2 \le \mu_F/\mu_R \le 2`. In the following we express this additional dependency as a function
+of :math:`k = \ln(\mu_F^2/\mu_R^2)`
 
 This variation can be performed at least at two different levels during the |PDF|
 evolution, always evaluating the strong coupling at :math:`\mu_R^2`.
@@ -169,30 +170,28 @@ evolution, always evaluating the strong coupling at :math:`\mu_R^2`.
       and the anomalous dimension are then modified using :cite:`Vogt:2004ns`:
 
         .. math ::
-            & \gamma^{(1)}(N) \to \gamma^{(1)}(N) - \beta_0 k \gamma^{(0)} \\
-            & \gamma^{(2)}(N) \to \gamma^{(2)}(N) - 2 \beta_0 k \gamma^{(1)} - ( \beta_1 k - \beta_0^2 k^2) \gamma^{(0)} \\
-            & \gamma^{(3)}(N) \to \gamma^{(3)}(N) - 3 \beta_0 k \gamma^{(2)} - ( 2 \beta_1 k - 3 \beta_0^2 k^2) \gamma^{(1)} - (\beta_2 k - \frac{5}{2} \beta_1 \beta_0 k^2 + \beta_0^3 k^3) \gamma^{(0)} \\
-            & k = \ln(\mu_F^2/\mu_R^2)
+            \gamma^{(1)}(N) &\to \gamma^{(1)}(N) - \beta_0 k \gamma^{(0)} \\
+            \gamma^{(2)}(N) &\to \gamma^{(2)}(N) - 2 \beta_0 k \gamma^{(1)} - ( \beta_1 k - \beta_0^2 k^2) \gamma^{(0)} \\
+            \gamma^{(3)}(N) &\to \gamma^{(3)}(N) - 3 \beta_0 k \gamma^{(2)} - ( 2 \beta_1 k - 3 \beta_0^2 k^2) \gamma^{(1)} - (\beta_2 k - \frac{5}{2} \beta_1 \beta_0 k^2 + \beta_0^3 k^3) \gamma^{(0)}
 
       This procedure corresponds to Eq. (3.32) of :cite:`AbdulKhalek:2019ihb`, and we recommend to use it along with
       ``ModEv='iterate-exact'`` in order to be in agreement with the treatment of the evolution integral expansion.
 
 
-    * In ``ModSV='expanded'`` the |EKO| is multiplied by an additional kernel, such that
-      the scale variation is applied to the whole evolution operator:
+    * In ``ModSV='expanded'`` the full |EKO| is multiplied by an additional kernel:
 
         .. math ::
             \tilde{\mathbf{E}}(a_s \leftarrow a_s^0) & = \tilde{\mathbf{K}}(a_s) \tilde{\mathbf{E}}(a_s \leftarrow a_s^0) \\
             \tilde{\mathbf{K}}(a_s) & = 1 - k \gamma + \frac{1}{2} k^2 \left ( \gamma^{2} - \beta \frac{\partial \gamma}{\partial a_s} \right ) \\
-            & + \frac{1}{6} k^3 \left [ - \beta \frac{\partial}{\partial a_s} \left( \beta \frac{\partial \gamma}{\partial a_s} \right) + 3 \beta \frac{\partial \gamma}{\partial a_s} \gamma - \gamma^3 \right ] + \mathcal{O}(k^4)
+            & \hspace{10pt} + \frac{1}{6} k^3 \left [ - \beta \frac{\partial}{\partial a_s} \left( \beta \frac{\partial \gamma}{\partial a_s} \right) + 3 \beta \frac{\partial \gamma}{\partial a_s} \gamma - \gamma^3 \right ] + \mathcal{O}(k^4)
 
-      where scale variation kernel is expanded consistently order by order in :math:`a_s`,
+      where the scale variation kernel :math:`\tilde{\mathbf{K}}` is expanded consistently order by order in :math:`a_s`,
       leading to:
 
         .. math ::
-            \tilde{\mathbf{K}}(a_s) \approx & 1 - a_s k \gamma^{(0)} + a_s^2 \left [ - k \gamma^{(1)} + \frac{1}{2} k^2 \gamma^{(0)} (\beta_0 + \gamma^{(0)}) \right ] \\
-            & + a_s^3 \left [ -k \gamma^{(2)} + \frac{1}{2} k^2 \left(\beta_1 \gamma^{(0)} + 2 \gamma^{(1)} (\beta_0 + \gamma^{(0)} ) \right) \right. \\
-            & \left. - \frac{1}{6} k^3 \gamma^{(0)} \left(2 \beta_0^2 + 3 \beta_0 \gamma^{(0)}+\gamma^{(0),2} \right) \right] + \mathcal{O}(a^4)
+            \tilde{\mathbf{K}}(a_s) &\approx 1 - a_s k \gamma^{(0)} + a_s^2 \left [ - k \gamma^{(1)} + \frac{1}{2} k^2 \gamma^{(0)} (\beta_0 + \gamma^{(0)}) \right ] \\
+            & \hspace{10pt} + a_s^3 \left [ -k \gamma^{(2)} + \frac{1}{2} k^2 \left(\beta_1 \gamma^{(0)} + 2 \beta_0\gamma^{(1)}  + \gamma^{(1)}\gamma^{(0)} + \gamma^{(0)}\gamma^{(1)} \right) \right. \\
+            & \hspace{35pt} \left. - \frac{1}{6} k^3 \gamma^{(0)} \left(2 \beta_0^2 + 3 \beta_0 \gamma^{(0)} + \left(\gamma^{(0)}\right)^2 \right) \right] + \mathcal{O}(a_s^4)
 
 
       In this way the dependence of the |EKO| on :math:`k` is factorized outside the unvaried evolution kernel.
@@ -200,12 +199,11 @@ evolution, always evaluating the strong coupling at :math:`\mu_R^2`.
       It corresponds to Eq. (3.35) of :cite:`AbdulKhalek:2019ihb`, and we recommend to use it along with
       ``ModEv='truncated'`` in order to keep consistency with the evolution integral expansion.
 
-
-By construction, the corrections of the order :math:`\mathcal{O}(k^n)` will appear
-at the order :math:`n` in the expansion :math:`a_s`.
-This happens because :math:`\beta \approx \mathcal{O}(a_s^2)`, :math:`\gamma \approx \mathcal{O}(a_s)`
-and the contribution proportional to :math:`\mathcal{O}(k^n)` is originated
-by the `n-th` derivative in :math:`\gamma` :cite:`AbdulKhalek:2019ihb`.
+      By construction, the corrections of the order :math:`\mathcal{O}(k^n)` will appear
+      at the order :math:`n` in the expansion :math:`a_s`.
+      This happens because :math:`\beta \approx \mathcal{O}(a_s^2)`, :math:`\gamma \approx \mathcal{O}(a_s)`
+      and the contribution proportional to :math:`\mathcal{O}(k^n)` is originated
+      by the `n-th` derivative in :math:`\gamma` :cite:`AbdulKhalek:2019ihb`.
 
 Furthermore the distance between the varied |EKO| and the unvaried one will decrease while
 keeping higher order terms in :math:`a_s`
diff --git a/src/eko/anomalous_dimensions/as4/ggg.py b/src/eko/anomalous_dimensions/as4/ggg.py
index 8e83a1d71..939cda60d 100644
--- a/src/eko/anomalous_dimensions/as4/ggg.py
+++ b/src/eko/anomalous_dimensions/as4/ggg.py
@@ -137,17 +137,17 @@ def gamma_gg_nf3(n, sx):
 
 
 @nb.njit(cache=True)
-def gamma_gg_nf2(n, sx):
+def gamma_gg_nf2(_n, _sx):
     return 0
 
 
 @nb.njit(cache=True)
-def gamma_gg_nf1(n, sx):
+def gamma_gg_nf1(_n, _sx):
     return 0
 
 
 @nb.njit(cache=True)
-def gamma_gg_nf0(n, sx):
+def gamma_gg_nf0(_n, _sx):
     return 0
 
 
diff --git a/src/eko/anomalous_dimensions/as4/gps.py b/src/eko/anomalous_dimensions/as4/gps.py
index 31b34f9e4..20878f063 100644
--- a/src/eko/anomalous_dimensions/as4/gps.py
+++ b/src/eko/anomalous_dimensions/as4/gps.py
@@ -79,12 +79,12 @@ def gamma_ps_nf3(n, sx):
 
 
 @nb.njit(cache=True)
-def gamma_ps_nf2(n, sx):
+def gamma_ps_nf2(_n, _sx):
     return 0
 
 
 @nb.njit(cache=True)
-def gamma_ps_nf1(n, sx):
+def gamma_ps_nf1(_n, _sx):
     return 0
 
 
diff --git a/src/eko/anomalous_dimensions/as4/gqg.py b/src/eko/anomalous_dimensions/as4/gqg.py
index 028875186..76047a249 100644
--- a/src/eko/anomalous_dimensions/as4/gqg.py
+++ b/src/eko/anomalous_dimensions/as4/gqg.py
@@ -329,17 +329,17 @@ def gamma_qg_nf3(n, sx):
 
 
 @nb.njit(cache=True)
-def gamma_qg_nf2(n, sx):
+def gamma_qg_nf2(_n, _sx):
     return 0
 
 
 @nb.njit(cache=True)
-def gamma_qg_nf1(n, sx):
+def gamma_qg_nf1(_n, _sx):
     return 0
 
 
 @nb.njit(cache=True)
-def gamma_qg_nf0(n, sx):
+def gamma_qg_nf0(_n, _sx):
     return 0
 
 
diff --git a/src/eko/evolution_operator/__init__.py b/src/eko/evolution_operator/__init__.py
index 9ab8fb5d7..5fa8def51 100644
--- a/src/eko/evolution_operator/__init__.py
+++ b/src/eko/evolution_operator/__init__.py
@@ -133,47 +133,49 @@ def quad_ker(
     ev_op_iterations,
     ev_op_max_order,
     sv_mode,
+    is_threshold,
 ):
-    """
-    Raw evolution kernel inside quad.
+    """Raw evolution kernel inside quad.
 
     Parameters
     ----------
-        u : float
-            quad argument
-        order : tuple(int,int)
-            perturbation orders
-        method : str
-            method
-        mode0: int
-            pid for first sector element
-        mode1 : int
-            pid for second sector element
-        is_log : boolean
-            is a logarithmic interpolation
-        logx : float
-            Mellin inversion point
-        areas : tuple
-            basis function configuration
-        as1 : float
-            target coupling value
-        as0 : float
-            initial coupling value
-        nf : int
-            number of active flavors
-        L : float
-            logarithm of the squared ratio of factorization and renormalization scale
-        ev_op_iterations : int
-            number of evolution steps
-        ev_op_max_order : int
-            perturbative expansion order of U
-        sv_mode: int, `enum.IntEnum`
-            scale variation mode, see `eko.scale_variations.Modes`
+    u : float
+        quad argument
+    order : int
+        perturbation order
+    mode0: int
+        pid for first sector element
+    mode1 : int
+        pid for second sector element
+    method : str
+        method
+    is_log : boolean
+        is a logarithmic interpolation
+    logx : float
+        Mellin inversion point
+    areas : tuple
+        basis function configuration
+    a1 : float
+        target coupling value
+    a0 : float
+        initial coupling value
+    nf : int
+        number of active flavors
+    L : float
+        logarithm of the squared ratio of factorization and renormalization scale
+    ev_op_iterations : int
+        number of evolution steps
+    ev_op_max_order : int
+        perturbative expansion order of U
+    sv_mode: int, `enum.IntEnum`
+        scale variation mode, see `eko.scale_variations.Modes`
+    is_threshold : boolean
+        is this an itermediate threshold operator?
 
     Returns
     -------
-        ker : float
-            evaluated integration kernel
+    float
+        evaluated integration kernel
     """
     ker_base = QuadKerBase(u, is_log, logx, mode0)
     integrand = ker_base.integrand(areas)
@@ -199,10 +201,10 @@ def quad_ker(
             ev_op_max_order,
         )
         # scale var expanded is applied on the kernel
-        if sv_mode == sv.Modes.expanded:
-            ker = np.ascontiguousarray(ker) @ np.ascontiguousarray(
+        if sv_mode == sv.Modes.expanded and not is_threshold:
+            ker = np.ascontiguousarray(
                 sv.expanded.singlet_variation(gamma_singlet, as1, order, nf, L)
-            )
+            ) @ np.ascontiguousarray(ker)
         ker = select_singlet_element(ker, mode0, mode1)
     else:
         gamma_ns = ad.gamma_ns(order, mode0, ker_base.n, nf)
@@ -217,40 +219,43 @@ def quad_ker(
             nf,
             ev_op_iterations,
         )
-        if sv_mode == sv.Modes.expanded:
-            ker = ker * sv.expanded.non_singlet_variation(gamma_ns, as1, order, nf, L)
+        if sv_mode == sv.Modes.expanded and not is_threshold:
+            ker = sv.expanded.non_singlet_variation(gamma_ns, as1, order, nf, L) * ker
 
     # recombine everthing
     return np.real(ker * integrand)
 
 
-class Operator:
-    """
-    Internal representation of a single EKO.
+class Operator(sv.ModeMixin):
+    """Internal representation of a single EKO.
 
     The actual matrices are computed upon calling :meth:`compute`.
 
     Parameters
     ----------
-        config : dict
-            configuration
-        managers : dict
-            managers
-        nf : int
-            number of active flavors
-        q2_from : float
-            evolution source
-        q2_to : float
-            evolution target
-        mellin_cut : float
-            cut to the upper limit in the mellin inversion
+    config : dict
+        configuration
+    managers : dict
+        managers
+    nf : int
+        number of active flavors
+    q2_from : float
+        evolution source
+    q2_to : float
+        evolution target
+    mellin_cut : float
+        cut to the upper limit in the mellin inversion
+    is_threshold : bool
+        is this an itermediate threshold operator?
     """
 
     log_label = "Evolution"
     # complete list of possible evolution operators labels
     full_labels = br.full_labels
 
-    def __init__(self, config, managers, nf, q2_from, q2_to=None, mellin_cut=5e-2):
+    def __init__(
+        self, config, managers, nf, q2_from, q2_to, mellin_cut=5e-2, is_threshold=False
+    ):
         self.config = config
         self.managers = managers
         self.nf = nf
@@ -258,6 +263,7 @@ def __init__(self, config, managers, nf, q2_from, q2_to=None, mellin_cut=5e-2):
         self.q2_to = q2_to
         # TODO make 'cut' external parameter?
         self._mellin_cut = mellin_cut
+        self.is_threshold = is_threshold
         self.op_members = {}
         self.order = config["order"]
 
@@ -266,20 +272,14 @@ def n_pools(self):
         n_pools = self.config["n_integration_cores"]
         if n_pools > 0:
             return n_pools
-        return os.cpu_count() + n_pools
+        # so we subtract from the maximum number
+        return max(os.cpu_count() + n_pools, 1)
 
     @property
     def fact_to_ren(self):
         r"""Returns the factor :math:`(\mu_F/\mu_R)^2`"""
         return self.config["fact_to_ren"]
 
-    @property
-    def sv_mode(self):
-        """Returns the scale variation mode"""
-        if self.config["ModSV"] is not None:
-            return sv.Modes[self.config["ModSV"]]
-        return sv.Modes.unvaried
-
     @property
     def int_disp(self):
         """Returns the interpolation dispatcher"""
@@ -290,27 +290,41 @@ def grid_size(self):
         """Returns the grid size"""
         return self.int_disp.xgrid.size
 
+    def mur2_shift(self, q2):
+        """Computes shifted renormalization scale.
+
+        Parameters
+        ----------
+        q2 : float
+            factorization scale
+
+        Returns
+        -------
+        float
+            renormalization scale
+        """
+        if self.sv_mode == sv.Modes.exponentiated:
+            return q2 / self.fact_to_ren
+        return q2
+
     @property
     def a_s(self):
         """Returns the computed values for :math:`a_s`"""
         sc = self.managers["strong_coupling"]
-        as0 = sc.a_s(
-            self.q2_from / self.fact_to_ren, fact_scale=self.q2_from, nf_to=self.nf
-        )
-        as1 = sc.a_s(
-            self.q2_to / self.fact_to_ren, fact_scale=self.q2_to, nf_to=self.nf
+        a0 = sc.a_s(
+            self.mur2_shift(self.q2_from), fact_scale=self.q2_from, nf_to=self.nf
         )
-        return (as0, as1)
+        a1 = sc.a_s(self.mur2_shift(self.q2_to), fact_scale=self.q2_to, nf_to=self.nf)
+        return (a0, a1)
 
     @property
     def labels(self):
-        """
-        Compute necessary sector labels to compute.
+        """Compute necessary sector labels to compute.
 
         Returns
         -------
-            labels : list(str)
-                sector labels
+        list(str)
+            sector labels
         """
         labels = []
         # the NS sector is dynamic
@@ -331,22 +345,21 @@ def labels(self):
         return labels
 
     def quad_ker(self, label, logx, areas):
-        """
-        Partially initialized integrand function
+        """Partially initialized integrand function.
 
         Parameters
         ----------
-            label: tuple
-                operator element pids
-            logx: float
-                Mellin inversion point
-            areas : tuple
-                basis function configuration
+        label: tuple
+            operator element pids
+        logx: float
+            Mellin inversion point
+        areas : tuple
+            basis function configuration
 
         Returns
         -------
-            quad_ker : functools.partial
-                partially initialized integration kernel
+        functools.partial
+            partially initialized integration kernel
 
         """
         return functools.partial(
@@ -366,6 +379,7 @@ def quad_ker(self, label, logx, areas):
             ev_op_iterations=self.config["ev_op_iterations"],
             ev_op_max_order=self.config["ev_op_max_order"],
             sv_mode=self.sv_mode,
+            is_threshold=self.is_threshold,
         )
 
     def initialize_op_members(self):
@@ -388,18 +402,17 @@ def run_op_integration(
         self,
         log_grid,
     ):
-        """
-        Run the integration for each grid point
+        """Run the integration for each grid point
 
         Parameters
         ----------
-            log_grid : tuple(k, logx)
-                log grid point with relative index
+        log_grid : tuple(k, logx)
+            log grid point with relative index
 
         Returns
         -------
-            column : list
-                computed operators at the give grid point
+        list
+            computed operators at the give grid point
 
         """
         column = []
@@ -434,7 +447,7 @@ def run_op_integration(
         return column
 
     def compute(self):
-        """compute the actual operators (i.e. run the integrations)"""
+        """Compute the actual operators (i.e. run the integrations)."""
         self.initialize_op_members()
 
         # skip computation ?
@@ -455,8 +468,8 @@ def compute(self):
         logger.info(
             "%s: µ_R^2 distance: %e -> %e",
             self.log_label,
-            self.q2_from / self.fact_to_ren,
-            self.q2_to / self.fact_to_ren,
+            self.mur2_shift(self.q2_from),
+            self.mur2_shift(self.q2_to),
         )
         if self.sv_mode.name != "unvaried":
             logger.info(
diff --git a/src/eko/evolution_operator/grid.py b/src/eko/evolution_operator/grid.py
index dbdc05fd8..db235f3d1 100644
--- a/src/eko/evolution_operator/grid.py
+++ b/src/eko/evolution_operator/grid.py
@@ -1,6 +1,5 @@
 # -*- coding: utf-8 -*-
-"""
-This module contains the :class:`OperatorGrid` class.
+"""This module contains the :class:`OperatorGrid` class.
 
 The first is the driver class of eko as it is the one that collects all the
 previously instantiated information and does the actual computation of the Q2s.
@@ -13,6 +12,7 @@
 
 from .. import basis_rotation as br
 from .. import matching_conditions, member
+from .. import scale_variations as sv
 from ..matching_conditions.operator_matrix_element import OperatorMatrixElement
 from ..thresholds import flavor_shift, is_downward_path
 from . import Operator, flavors, physical
@@ -20,7 +20,7 @@
 logger = logging.getLogger(__name__)
 
 
-class OperatorGrid:
+class OperatorGrid(sv.ModeMixin):
     """
     The operator grid is the driver class of the evolution.
 
@@ -185,7 +185,12 @@ def get_threshold_operators(self, path):
                 # Compute the operator and store it
                 logger.info("Prepare threshold operator")
                 op_th = Operator(
-                    self.config, self.managers, seg.nf, seg.q2_from, seg.q2_to
+                    self.config,
+                    self.managers,
+                    seg.nf,
+                    seg.q2_from,
+                    seg.q2_to,
+                    is_threshold=True,
                 )
                 op_th.compute()
                 self._threshold_operators[new_op_key] = op_th
@@ -198,18 +203,17 @@ def get_threshold_operators(self, path):
         return thr_ops
 
     def compute(self, q2grid=None):
-        """
-        Computes all ekos for the q2grid.
+        """Computes all ekos for the `q2grid`.
 
         Parameters
         ----------
-            q2grid: list(float)
-                List of q^2
+        q2grid: list(float)
+            List of :math:`Q^2`
 
         Returns
         -------
-            grid_return: list(dict)
-                List of ekos for each value of q^2
+        list(dict)
+            List of ekos for each value of :math:`Q^2`
         """
         # use input?
         if q2grid is None:
@@ -220,22 +224,27 @@ def compute(self, q2grid=None):
         # And now return the grid
         grid_return = {}
         for q2 in q2grid:
-            grid_return[q2] = self.generate(q2)
+            # shift path for expanded scheme
+            q2_gen = (
+                q2 * self.config["fact_to_ren"]
+                if self.sv_mode == sv.Modes.expanded
+                else q2
+            )
+            grid_return[q2] = self.generate(q2_gen)
         return grid_return
 
     def generate(self, q2):
-        """
-        Computes a single EKO.
+        r"""Computes a single EKO.
 
         Parameters
         ----------
-            q2: float
-                Target value of q^2
+        q2: float
+            Target value of :math:`Q^2`
 
         Returns
         -------
-            final_op: dict
-                eko E(q^2 <- q_0^2) in flavor basis as numpy array
+        dict
+            eko :math:`\mathbf E(Q^2 \leftarrow Q_0^2)` in flavor basis as numpy array
         """
         # The lists of areas as produced by the thresholds
         path = self.managers["thresholds_config"].path(q2)
@@ -287,13 +296,7 @@ def generate(self, q2):
                 final_op = final_op @ rot @ matching @ phys_op
 
         values, errors = final_op.to_flavor_basis_tensor()
-        fact_to_ren = self.config["fact_to_ren"]
         return {
             "operators": values,
             "operator_errors": errors,
-            "alphas": self.managers["strong_coupling"].a(
-                q2 / fact_to_ren, fact_scale=q2, nf_to=path[-1].nf
-            )[0]
-            * 4.0
-            * np.pi,
         }
diff --git a/src/eko/harmonics/f_functions/f11.py b/src/eko/harmonics/f_functions/f11.py
index 5ec1984e4..b0e46365e 100644
--- a/src/eko/harmonics/f_functions/f11.py
+++ b/src/eko/harmonics/f_functions/f11.py
@@ -1,4 +1,5 @@
 # -*- coding: utf-8 -*-
+# pylint: skip-file
 """This module contains implemtation of F11.
 
 When using it, please cite :cite:`Blumlein:2009ta`.
diff --git a/src/eko/harmonics/f_functions/f13.py b/src/eko/harmonics/f_functions/f13.py
index f540a834e..38046f88e 100644
--- a/src/eko/harmonics/f_functions/f13.py
+++ b/src/eko/harmonics/f_functions/f13.py
@@ -1,4 +1,5 @@
 # -*- coding: utf-8 -*-
+# pylint: skip-file
 """This module contains implemtation of F13.
 
 When using it, please cite :cite:`Blumlein:2009ta`.
diff --git a/src/eko/harmonics/f_functions/f14_f12.py b/src/eko/harmonics/f_functions/f14_f12.py
index de0745bb7..48c8b1168 100644
--- a/src/eko/harmonics/f_functions/f14_f12.py
+++ b/src/eko/harmonics/f_functions/f14_f12.py
@@ -1,5 +1,5 @@
 # -*- coding: utf-8 -*-
-# pylint: disable=too-many-lines
+# pylint: skip-file
 """This module contains implemtation of F14 - F12.
 
 When using it, please cite :cite:`Blumlein:2009ta`.
diff --git a/src/eko/harmonics/f_functions/f16.py b/src/eko/harmonics/f_functions/f16.py
index db9ecd017..26caa9353 100644
--- a/src/eko/harmonics/f_functions/f16.py
+++ b/src/eko/harmonics/f_functions/f16.py
@@ -1,5 +1,5 @@
 # -*- coding: utf-8 -*-
-# pylint: disable=too-many-lines
+# pylint: skip-file
 """This module contains implemtation of F16.
 
 When using it, please cite :cite:`Blumlein:2009ta`.
diff --git a/src/eko/harmonics/f_functions/f17.py b/src/eko/harmonics/f_functions/f17.py
index cc0f0a4cd..be9709821 100644
--- a/src/eko/harmonics/f_functions/f17.py
+++ b/src/eko/harmonics/f_functions/f17.py
@@ -1,4 +1,5 @@
 # -*- coding: utf-8 -*-
+# pylint: skip-file
 """This module contains implemtation of F17.
 
 When using it, please cite :cite:`Blumlein:2009ta`.
diff --git a/src/eko/harmonics/f_functions/f18.py b/src/eko/harmonics/f_functions/f18.py
index 344a2fe84..0a5d55948 100644
--- a/src/eko/harmonics/f_functions/f18.py
+++ b/src/eko/harmonics/f_functions/f18.py
@@ -1,5 +1,5 @@
 # -*- coding: utf-8 -*-
-# pylint: disable=too-many-lines
+# pylint: skip-file
 """This module contains implemtation of F18.
 
 When using it, please cite :cite:`Blumlein:2009ta`.
diff --git a/src/eko/harmonics/f_functions/f19.py b/src/eko/harmonics/f_functions/f19.py
index f372b89d9..76a180194 100644
--- a/src/eko/harmonics/f_functions/f19.py
+++ b/src/eko/harmonics/f_functions/f19.py
@@ -1,4 +1,5 @@
 # -*- coding: utf-8 -*-
+# pylint: skip-file
 """This module contains implemtation of F19.
 
 When using it, please cite :cite:`Blumlein:2009ta`.
diff --git a/src/eko/harmonics/f_functions/f20.py b/src/eko/harmonics/f_functions/f20.py
index 2db510935..23149d287 100644
--- a/src/eko/harmonics/f_functions/f20.py
+++ b/src/eko/harmonics/f_functions/f20.py
@@ -1,4 +1,5 @@
 # -*- coding: utf-8 -*-
+# pylint: skip-file
 """This module contains implemtation of F20.
 
 When using it, please cite :cite:`Blumlein:2009ta`.
diff --git a/src/eko/harmonics/f_functions/f21.py b/src/eko/harmonics/f_functions/f21.py
index a2aba7b90..81895821b 100644
--- a/src/eko/harmonics/f_functions/f21.py
+++ b/src/eko/harmonics/f_functions/f21.py
@@ -1,4 +1,5 @@
 # -*- coding: utf-8 -*-
+# pylint: skip-file
 """This module contains implemtation of F21.
 
 When using it, please cite :cite:`Blumlein:2009ta`.
diff --git a/src/eko/harmonics/f_functions/f9.py b/src/eko/harmonics/f_functions/f9.py
index 94d41fb2e..3099e9679 100644
--- a/src/eko/harmonics/f_functions/f9.py
+++ b/src/eko/harmonics/f_functions/f9.py
@@ -1,4 +1,5 @@
 # -*- coding: utf-8 -*-
+# pylint: skip-file
 """This module contains implemtation of F9.
 
 When using it, please cite :cite:`Blumlein:2009ta`.
diff --git a/src/eko/matching_conditions/as3/aHg.py b/src/eko/matching_conditions/as3/aHg.py
index 7482accd3..a7ea6dbbb 100644
--- a/src/eko/matching_conditions/as3/aHg.py
+++ b/src/eko/matching_conditions/as3/aHg.py
@@ -1,5 +1,5 @@
 # -*- coding: utf-8 -*-
-# pylint: disable=too-many-lines
+# pylint: skip-file
 import numba as nb
 import numpy as np
 
diff --git a/src/eko/matching_conditions/as3/aHgstfac.py b/src/eko/matching_conditions/as3/aHgstfac.py
index 2ddd57d70..8651be75c 100644
--- a/src/eko/matching_conditions/as3/aHgstfac.py
+++ b/src/eko/matching_conditions/as3/aHgstfac.py
@@ -1,4 +1,5 @@
 # -*- coding: utf-8 -*-
+# pylint: skip-file
 import numba as nb
 import numpy as np
 
diff --git a/src/eko/matching_conditions/as3/agg.py b/src/eko/matching_conditions/as3/agg.py
index e0bf47284..e5fd62576 100644
--- a/src/eko/matching_conditions/as3/agg.py
+++ b/src/eko/matching_conditions/as3/agg.py
@@ -1,5 +1,5 @@
 # -*- coding: utf-8 -*-
-# pylint: disable=too-many-lines
+# pylint: skip-file
 import numba as nb
 import numpy as np
 
diff --git a/src/eko/matching_conditions/operator_matrix_element.py b/src/eko/matching_conditions/operator_matrix_element.py
index 4e6533e0a..117308146 100644
--- a/src/eko/matching_conditions/operator_matrix_element.py
+++ b/src/eko/matching_conditions/operator_matrix_element.py
@@ -281,27 +281,26 @@ def quad_ker(
 
 
 class OperatorMatrixElement(Operator):
-    """
-    Internal representation of a single |OME|.
+    """Internal representation of a single |OME|.
 
     The actual matrices are computed upon calling :meth:`compute`.
 
     Parameters
     ----------
-        config : dict
-            configuration
-        managers : dict
-            managers
-        is_backward: bool
-            True for backward evolution
-        q2: float
-            matching scale
-        nf: int
-            number of active flavor below threshold
-        L: float
-            log of K threshold squared
-        is_msbar: bool
-            add the |MSbar| contribution
+    config : dict
+        configuration
+    managers : dict
+        managers
+    nf: int
+        number of active flavor below threshold
+    q2: float
+        matching scale
+    is_backward: bool
+        True for backward evolution
+    L: float
+        log of K threshold squared
+    is_msbar: bool
+        add the |MSbar| contribution
     """
 
     log_label = "Matching"
@@ -320,7 +319,7 @@ class OperatorMatrixElement(Operator):
     ]
 
     def __init__(self, config, managers, nf, q2, is_backward, L, is_msbar):
-        super().__init__(config, managers, nf, q2)
+        super().__init__(config, managers, nf, q2, None)
         self.backward_method = config["backward_inversion"] if is_backward else ""
         if is_backward:
             self.is_intrinsic = True
@@ -333,13 +332,12 @@ def __init__(self, config, managers, nf, q2, is_backward, L, is_msbar):
 
     @property
     def labels(self):
-        """
-        Compute necessary sector labels to compute.
+        """Computes the necessary sector labels to compute.
 
         Returns
         -------
-            labels : list(str)
-                sector labels
+        list(str)
+            sector labels
         """
 
         labels = []
@@ -375,23 +373,21 @@ def labels(self):
         return labels
 
     def quad_ker(self, label, logx, areas):
-        """
-        Partially initialized integrand function
+        """Partially initialized integrand function.
 
         Parameters
         ----------
-            label: tuple
-                operator element pids
-            logx: float
-                Mellin inversion point
-            areas : tuple
-                basis function configuration
+        label: tuple
+            operator element pids
+        logx: float
+            Mellin inversion point
+        areas : tuple
+            basis function configuration
 
         Returns
         -------
-            quad_ker : functools.partial
-                partially initialized integration kernel
-
+        functools.partial
+            partially initialized integration kernel
         """
         return functools.partial(
             quad_ker,
@@ -410,17 +406,15 @@ def quad_ker(self, label, logx, areas):
 
     @property
     def a_s(self):
-        """
-        Returns the computed values for :math:`a_s`.
+        """Returns the computed values for :math:`a_s`.
+
         Note that here you need to use :math:`a_s^{n_f+1}`
         """
         sc = self.managers["strong_coupling"]
-        return sc.a(self.q2_from / self.fact_to_ren, self.q2_from, nf_to=self.nf + 1)[0]
+        return sc.a_s(self.mur2_shift(self.q2_from), self.q2_from, nf_to=self.nf + 1)
 
     def compute(self):
-        """
-        compute the actual operators (i.e. run the integrations)
-        """
+        """Compute the actual operators (i.e. run the integrations)"""
         self.initialize_op_members()
 
         # At LO you don't need anything else
diff --git a/src/eko/scale_variations/__init__.py b/src/eko/scale_variations/__init__.py
index 853d95ce7..b16fac3e4 100644
--- a/src/eko/scale_variations/__init__.py
+++ b/src/eko/scale_variations/__init__.py
@@ -1,7 +1,5 @@
 # -*- coding: utf-8 -*-
-"""
-This module contain the possible scale variations integrals.
-"""
+"""This module contain the possible scale variations schemes."""
 import enum
 
 from . import expanded, exponentiated
@@ -13,3 +11,30 @@ class Modes(enum.IntEnum):
     unvaried = enum.auto()
     exponentiated = enum.auto()
     expanded = enum.auto()
+
+
+def sv_mode(s):
+    """Returns the scale variation mode.
+
+    Parameters
+    ----------
+    s : str
+        string representation
+
+    Returns
+    -------
+    enum.IntEnum
+        enum representation
+    """
+    if s is not None:
+        return Modes[s]
+    return Modes.unvaried
+
+
+class ModeMixin:
+    """Mixin to cast scale variation mode."""
+
+    @property
+    def sv_mode(self):
+        """Returns the scale variation mode"""
+        return sv_mode(self.config["ModSV"])
diff --git a/src/eko/scale_variations/expanded.py b/src/eko/scale_variations/expanded.py
index 7cdd12379..aeea77b76 100644
--- a/src/eko/scale_variations/expanded.py
+++ b/src/eko/scale_variations/expanded.py
@@ -1,8 +1,11 @@
 # -*- coding: utf-8 -*-
-r"""
-This module contains the scale variation operator in ``ModSV=expanded``
-"""
+r"""This module contains the scale variation operator for the expanded scheme (``ModSV=expanded``).
 
+The expressions can be obtained using Eqs. (3.33) and (3.38) of :cite:`AbdulKhalek:2019ihb`.
+Be aware that corresponding the signs of the ingredients there are a number of differences.
+However, the ultimate sign can be obtained by comparing to the exponentiated scheme in the
+trunctated solution.
+"""
 
 import numba as nb
 import numpy as np
@@ -11,157 +14,160 @@
 
 
 @nb.njit(cache=True)
-def gamma_1_variation(gamma, L):
-    r"""
-    Computes the |NLO| anomalous dimension variation.
+def variation_as1(gamma, L):
+    r"""Computes the |NLO| anomalous dimension variation.
 
     Parameters
     ----------
-        gamma : numpy.ndarray
-            anomalous dimensions
-        L : float
-            logarithmic ratio of factorization and renormalization scale
+    gamma : numpy.ndarray
+        anomalous dimensions
+    L : float
+        logarithmic ratio of factorization and renormalization scale
 
     Returns
     -------
-        gamma_1 : complex
-            variation to :math:`\gamma^{(1)}`
+    complex
+        variation at |NLO|
     """
     return -L * gamma[0]
 
 
 @nb.njit(cache=True)
-def gamma_2_variation(gamma, L, beta0, g0e2):
-    r"""
-    Computes the |NNLO| anomalous dimension variation.
+def variation_as2(gamma, L, beta0, g0e2):
+    r"""Computes the |NNLO| anomalous dimension variation.
 
     Parameters
     ----------
-        gamma : numpy.ndarray
-            anomalous dimensions
-        L : float
-            logarithmic ratio of factorization and renormalization scale
-        beta0: float
-            :math:`\beta_0`
-        g0e2: complex
-            :math:`\gamma^{(0),2}`
+    gamma : numpy.ndarray
+        anomalous dimensions
+    L : float
+        logarithmic ratio of factorization and renormalization scale
+    beta0: float
+        :math:`\beta_0`
+    g0e2: complex or numpy.ndarray
+        :math:`\left(\gamma^{(0)}\right)^2`
 
     Returns
     -------
-        gamma_2 : complex
-            variation to :math:`\gamma^{(2)}`
+    complex
+        variation at |NNLO|
     """
-    return -gamma[1] * L + 1 / 2 * (beta0 * gamma[0] + g0e2) * L**2
+    return -gamma[1] * L + 1.0 / 2.0 * (beta0 * gamma[0] + g0e2) * L**2
 
 
 @nb.njit(cache=True)
-def gamma_3_variation(gamma, L, beta0, beta1, g0e2, g0e3, g1g0):
-    r"""
-    Computes the |N3LO| anomalous dimension variation.
+def variation_as3(gamma, L, beta0, beta1, g0e2, g0e3, g1g0, g0g1):
+    r"""Computes the |N3LO| anomalous dimension variation.
 
     Parameters
     ----------
-        gamma : numpy.ndarray
-            anomalous dimensions
-        L : float
-            logarithmic ratio of factorization and renormalization scale
-        beta0: float
-            :math:`\beta_0`
-        beta0: float
-            :math:`\beta_1`
-        g0e2: complex
-            :math:`\gamma^{(0),2}`
-        g0e3: complex
-            :math:`\gamma^{(0),3}`
-        g1g0: complex
-            :math:`\gamma^{(1)} \gamma^{(0)}`
+    gamma : numpy.ndarray
+        anomalous dimensions
+    L : float
+        logarithmic ratio of factorization and renormalization scale
+    beta0: float
+        :math:`\beta_0`
+    beta0: float
+        :math:`\beta_1`
+    g0e2: complex or numpy.ndarray
+        :math:`\left(\gamma^{(0)}\right)^2`
+    g0e3: complex or numpy.ndarray
+        :math:`\left(\gamma^{(0)}\right)^3`
+    g1g0: complex or numpy.ndarray
+        :math:`\gamma^{(1)} \gamma^{(0)}`
+    g0g1: complex or numpy.ndarray
+        :math:`\gamma^{(0)} \gamma^{(1)}`
 
     Returns
     -------
-        gamma_3 : complex
-            variation to :math:`\gamma^{(3)}`
+    complex
+        variation at |N3LO|
     """
     return (
         -gamma[2] * L
-        + (1 / 2) * (beta1 * gamma[0] + 2 * beta0 * gamma[1] + 2 * g1g0) * L**2
-        - (1 / 6) * (2 * beta0**2 * gamma[0] + 3 * beta0 * g0e2 + g0e3) * L**3
+        + (1.0 / 2.0)
+        * (beta1 * gamma[0] + 2.0 * beta0 * gamma[1] + g1g0 + g0g1)
+        * L**2
+        - (1.0 / 6.0)
+        * (2.0 * beta0**2 * gamma[0] + 3.0 * beta0 * g0e2 + g0e3)
+        * L**3
     )
 
 
 @nb.njit(cache=True)
 def non_singlet_variation(gamma, a_s, order, nf, L):
-    """
-    Scale Variation non-singlet dispatcher
+    """Non-singlet scale variation dispatcher.
 
     Parameters
     ----------
-        gamma : numpy.ndarray
-            anomalous dimensions
-        a_s :  float
-            target coupling value
-        order : tuple(int,int)
-            perturbation order
-        nf : int
-            number of active flavors
-        L : float
-            logarithmic ratio of factorization and renormalization scale
+    gamma : numpy.ndarray
+        anomalous dimensions
+    a_s :  float
+        target coupling value
+    order : int
+        perturbation order
+    nf : int
+        number of active flavors
+    L : float
+        logarithmic ratio of factorization and renormalization scale
 
     Returns
     -------
-        sv_ker : numpy.ndarray
-            scale varion kernel
+    complex
+        scale variation kernel
     """
     sv_ker = 1.0
     if order[0] >= 2:
-        sv_ker += a_s * gamma_1_variation(gamma, L)
+        sv_ker += a_s * variation_as1(gamma, L)
     if order[0] >= 3:
         beta0 = beta.beta_qcd_as2(nf)
-        sv_ker += a_s**2 * gamma_2_variation(gamma, L, beta0, gamma[0] ** 2)
+        sv_ker += a_s**2 * variation_as2(gamma, L, beta0, gamma[0] ** 2)
     if order[0] >= 4:
         beta1 = beta.beta_qcd((3, 0), nf)
-        sv_ker += a_s**3 * gamma_3_variation(
-            gamma, L, beta0, beta1, gamma[0] ** 2, gamma[0] ** 3, gamma[0] * gamma[1]
+        g0g1 = gamma[0] * gamma[1]
+        sv_ker += a_s**3 * variation_as3(
+            gamma, L, beta0, beta1, gamma[0] ** 2, gamma[0] ** 3, g0g1, g0g1
         )
     return sv_ker
 
 
 @nb.njit(cache=True)
 def singlet_variation(gamma, a_s, order, nf, L):
-    """
-    Scale Variation singlet dispatcher
+    """Singlet scale variation dispatcher.
 
     Parameters
     ----------
-        gamma : numpy.ndarray
-            anomalous dimensions
-        a_s :  float
-            target coupling value
-        order : tuple(int,int)
-            perturbation order
-        nf : int
-            number of active flavors
-        L : float
-            logarithmic ratio of factorization and renormalization scale
+    gamma : numpy.ndarray
+        anomalous dimensions
+    a_s :  float
+        target coupling value
+    order : int
+        perturbation order
+    nf : int
+        number of active flavors
+    L : float
+        logarithmic ratio of factorization and renormalization scale
 
     Returns
     -------
-        sv_ker : numpy.ndarray
-            scale varion kernel
+    numpy.ndarray
+        scale variation kernel
     """
     sv_ker = np.eye(2, dtype=np.complex_)
     gamma = np.ascontiguousarray(gamma)
     if order[0] >= 2:
-        sv_ker += a_s * gamma_1_variation(gamma, L)
+        sv_ker += a_s * variation_as1(gamma, L)
     if order[0] >= 3:
         beta0 = beta.beta_qcd_as2(nf)
         gamma0e2 = gamma[0] @ gamma[0]
-        sv_ker += a_s**2 * gamma_2_variation(gamma, L, beta0, gamma0e2)
+        sv_ker += a_s**2 * variation_as2(gamma, L, beta0, gamma0e2)
     if order[0] >= 4:
         beta1 = beta.beta_qcd((3, 0), nf)
         gamma0e3 = gamma0e2 @ gamma[0]
         # here the product is not commutative
         g1g0 = gamma[1] @ gamma[0]
-        sv_ker += a_s**3 * gamma_3_variation(
-            gamma, L, beta0, beta1, gamma0e2, gamma0e3, g1g0
+        g0g1 = gamma[0] @ gamma[1]
+        sv_ker += a_s**3 * variation_as3(
+            gamma, L, beta0, beta1, gamma0e2, gamma0e3, g1g0, g0g1
         )
     return sv_ker
diff --git a/src/ekobox/utils.py b/src/ekobox/utils.py
index 858525f0b..e228d7cc4 100644
--- a/src/ekobox/utils.py
+++ b/src/ekobox/utils.py
@@ -40,7 +40,6 @@ def ekos_product(eko_ini, eko_fin, in_place=True):
 
     final_op_dict = {}
     final_op_error_dict = {}
-    final_alphas_dict = {}
     final_dict = {}
     for q2, op2 in ope2_dict.items():
         final_op_dict[q2] = np.einsum("ajbk,bkcl -> ajcl", ope1, op2)
@@ -49,11 +48,9 @@ def ekos_product(eko_ini, eko_fin, in_place=True):
             "ajbk,bkcl -> ajcl", ope1, ope2_error_dict[q2]
         ) + np.einsum("ajbk,bkcl -> ajcl", ope1_error, op2)
 
-        final_alphas_dict[q2] = eko_fin["Q2grid"][q2]["alphas"]
         final_dict[q2] = {
             "operators": final_op_dict[q2],
             "operator_errors": final_op_error_dict[q2],
-            "alphas": final_alphas_dict[q2],
         }
 
     final_eko = None
diff --git a/tests/eko/test_couplings.py b/tests/eko/test_couplings.py
index 21452f154..ede6a190c 100644
--- a/tests/eko/test_couplings.py
+++ b/tests/eko/test_couplings.py
@@ -77,6 +77,7 @@ def test_init(self):
                         ktThr=1.0,
                         MaxNfAs=6,
                         HQ="POLE",
+                        ModSV=None,
                     )
                     sc2 = Couplings.from_dict(setup)
                     assert sc2.q2_ref == scale_ref
diff --git a/tests/eko/test_ev_op_grid.py b/tests/eko/test_ev_op_grid.py
index 197778bd1..0b6fdbbc0 100644
--- a/tests/eko/test_ev_op_grid.py
+++ b/tests/eko/test_ev_op_grid.py
@@ -109,16 +109,12 @@ def test_compute_q2grid(self):
         opg = opgrid.compute()
         assert len(opg) == 2
         assert all(
-            k in op
-            for k in ["operators", "operator_errors", "alphas"]
-            for op in opg.values()
+            k in op for k in ["operators", "operator_errors"] for op in opg.values()
         )
         opg = opgrid.compute(3)
         assert len(opg) == 1
         assert all(
-            k in op
-            for k in ["operators", "operator_errors", "alphas"]
-            for op in opg.values()
+            k in op for k in ["operators", "operator_errors"] for op in opg.values()
         )
 
     def test_grid_computation_VFNS(self):
@@ -128,16 +124,6 @@ def test_grid_computation_VFNS(self):
         operators = opgrid.compute(qgrid_check)
         assert len(operators) == len(qgrid_check)
 
-    def test_alphas(self):
-        opgrid = self._get_operator_grid()
-        # q2 has not be precomputed - but should work nevertheless
-        opg = opgrid.compute(3)
-        sv_opgrid = self._get_operator_grid(
-            theory_update={"fact_to_ren_scale_ratio": 2.0}
-        )
-        sv_opg = sv_opgrid.compute(3)
-        assert opg[3]["alphas"] < sv_opg[3]["alphas"]
-
     def test_mod_expanded(self):
         theory_update = {
             "order": (1, 0),
diff --git a/tests/eko/test_ev_operator.py b/tests/eko/test_ev_operator.py
index d2ba997e8..795712c6c 100644
--- a/tests/eko/test_ev_operator.py
+++ b/tests/eko/test_ev_operator.py
@@ -45,6 +45,7 @@ def test_quad_ker(monkeypatch):
             ev_op_iterations=0,
             ev_op_max_order=(0, 0),
             sv_mode=1,
+            is_threshold=False,
         )
         np.testing.assert_allclose(res_ns, 0.0)
         res_s = quad_ker(
@@ -63,6 +64,7 @@ def test_quad_ker(monkeypatch):
             ev_op_iterations=0,
             ev_op_max_order=(0, 0),
             sv_mode=1,
+            is_threshold=False,
         )
         np.testing.assert_allclose(res_s, 1.0)
         res_s = quad_ker(
@@ -81,6 +83,7 @@ def test_quad_ker(monkeypatch):
             ev_op_iterations=0,
             ev_op_max_order=(0, 0),
             sv_mode=1,
+            is_threshold=False,
         )
         np.testing.assert_allclose(res_s, 0.0)
     for label in [(br.non_singlet_pids_map["ns+"], 0), (100, 100)]:
@@ -101,6 +104,7 @@ def test_quad_ker(monkeypatch):
                 ev_op_iterations=0,
                 ev_op_max_order=(1, 0),
                 sv_mode=sv,
+                is_threshold=False,
             )
             np.testing.assert_allclose(res_sv, 1.0)
 
@@ -121,6 +125,7 @@ def test_quad_ker(monkeypatch):
         ev_op_iterations=0,
         ev_op_max_order=(0, 0),
         sv_mode=1,
+        is_threshold=False,
     )
     np.testing.assert_allclose(res_ns, 0.0)
 
@@ -195,6 +200,9 @@ def test_labels(self):
 
     def test_n_pools(self):
         excluded_cores = 3
+        # make sure we actually have more the those cores (e.g. on github we don't)
+        if os.cpu_count() <= excluded_cores:
+            return
         o = Operator(
             dict(
                 order=(2, 0),
@@ -205,9 +213,28 @@ def test_n_pools(self):
             {},
             3,
             1,
+            10,
         )
         assert o.n_pools == os.cpu_count() - excluded_cores
 
+    def test_exponentiated(self):
+        tcard = copy.deepcopy(theory_card)
+        tcard["fact_to_ren_scale_ratio"] = 2.0
+        tcard["ModSV"] = "exponentiated"
+        ocard = copy.deepcopy(operators_card)
+        g = OperatorGrid.from_dict(
+            tcard,
+            ocard,
+            ThresholdsAtlas.from_dict(tcard),
+            Couplings.from_dict(tcard),
+            InterpolatorDispatcher.from_dict(ocard),
+        )
+        # setup objs
+        o = Operator(g.config, g.managers, 3, 2.0, 10.0)
+        np.testing.assert_allclose(o.mur2_shift(40.0), 10.0)
+        o.compute()
+        self.check_lo(o)
+
     def test_compute_parallel(self, monkeypatch):
         tcard = copy.deepcopy(theory_card)
         ocard = copy.deepcopy(operators_card)
@@ -337,6 +364,7 @@ def quad_ker_pegasus(
                     ev_op_iterations,
                     10,
                     0,
+                    False,
                 ),
                 epsabs=1e-12,
                 epsrel=1e-5,
diff --git a/tests/eko/test_msbar_masses.py b/tests/eko/test_msbar_masses.py
index aef3de6da..88279ad26 100644
--- a/tests/eko/test_msbar_masses.py
+++ b/tests/eko/test_msbar_masses.py
@@ -30,6 +30,7 @@
     "Qmt": 174.9,
     "order": (3, 0),
     "ModEv": "TRN",
+    "ModSV": None,
 }
 
 
diff --git a/tests/eko/test_sv_expanded.py b/tests/eko/test_sv_expanded.py
index c13107868..eae6e9b96 100644
--- a/tests/eko/test_sv_expanded.py
+++ b/tests/eko/test_sv_expanded.py
@@ -63,7 +63,7 @@ def scheme_diff(g, k, pto, is_singlet):
         """
         :math:`(ker_A - ker_B)/ker_{unv}` for truncated expansion
         Effects due to non commutativity are neglected thus,
-        he accuracy of singlet quantities is slightly worst.
+        the accuracy of singlet quantities is slightly worse.
         """
         if pto[0] >= 2:
             diff = g[0] * k * a0