diff --git a/conda-recipe/meta.yaml b/conda-recipe/meta.yaml
index 23b2098e1e..60cd770671 100644
--- a/conda-recipe/meta.yaml
+++ b/conda-recipe/meta.yaml
@@ -52,7 +52,7 @@ requirements:
         - sphinxcontrib-bibtex
         - curio >=1.0
         - pineappl >=0.6.2
-        - eko >=0.14.0
+        - eko >=0.14.1
         - fiatlux
 
 test:
diff --git a/nnpdfcpp/data/theory.db b/nnpdfcpp/data/theory.db
index 1462a62b8a..e8169b528a 100644
Binary files a/nnpdfcpp/data/theory.db and b/nnpdfcpp/data/theory.db differ
diff --git a/validphys2/src/validphys/photon/alpha.py b/validphys2/src/validphys/photon/alpha.py
new file mode 100644
index 0000000000..574f12d67c
--- /dev/null
+++ b/validphys2/src/validphys/photon/alpha.py
@@ -0,0 +1,254 @@
+"""Module that implements the alpha running"""
+import numpy as np
+from scipy.integrate import solve_ivp
+
+from eko import beta
+from n3fit.io.writer import XGRID
+
+from .constants import ME, MMU, MQL, MTAU
+
+
+class Alpha:
+    """Class that solves the RGE for alpha_qed"""
+    def __init__(self, theory, q2max):
+        self.theory = theory
+        self.alpha_em_ref = theory["alphaqed"]
+        self.qref = theory["Qref"]
+        self.qed_order = theory['QED']
+
+        # compute and store thresholds
+        self.thresh_c = self.theory["kcThr"] * self.theory["mc"]
+        self.thresh_b = self.theory["kbThr"] * self.theory["mb"]
+        self.thresh_t = self.theory["ktThr"] * self.theory["mt"]
+        if self.theory["MaxNfAs"] <= 5:
+            self.thresh_t = np.inf
+        if self.theory["MaxNfAs"] <= 4:
+            self.thresh_b = np.inf
+        if self.theory["MaxNfAs"] <= 3:
+            self.thresh_c = np.inf
+
+        if self.theory["ModEv"] == "TRN":
+            self.alphaem_fixed_flavor = self.alphaem_fixed_flavor_trn
+            self.thresh, self.alphaem_thresh = self.compute_alphaem_at_thresholds()
+        elif self.theory["ModEv"] == "EXA":
+            self.alphaem_fixed_flavor = self.alphaem_fixed_flavor_exa
+            self.thresh, self.alphaem_thresh = self.compute_alphaem_at_thresholds()
+
+            xmin = XGRID[0]
+            qmin = xmin * theory["MP"] / np.sqrt(1 - xmin)
+            # use a lot of interpolation points since it is a long path 1e-9 -> 1e4
+            self.q = np.geomspace(qmin, np.sqrt(q2max), 500, endpoint=True)
+
+            # add threshold points in the q list since alpha is not smooth there
+            self.q = np.append(
+                self.q, [ME, MMU, MQL, self.thresh_c, MTAU, self.thresh_b, self.qref, self.thresh_t]
+            )
+            self.q = self.q[np.isfinite(self.q)]
+            self.q.sort()
+
+            self.alpha_vec = np.array([self.alpha_em(q_) for q_ in self.q])
+            self.alpha_em = self.interpolate_alphaem
+        else:
+            raise ValueError(f"Evolution mode not recognized: {self.theory['ModEv']}")
+
+    def interpolate_alphaem(self, q):
+        r"""
+        Interpolate precomputed values of alpha_em.
+
+        Parameters
+        ----------
+        q: float
+            value in which alpha_em is computed
+
+        Returns
+        -------
+        alpha_em: float
+            electromagnetic coupling
+        """
+        return np.interp(q, self.q, self.alpha_vec)
+
+    def alpha_em(self, q):
+        r"""
+        Compute the value of the running alphaem.
+
+        Parameters
+        ----------
+        q: float
+            value in which alphaem is computed
+
+        Returns
+        -------
+        alpha_em: numpy.ndarray
+            electromagnetic coupling
+        """
+        nf, nl = self.find_region(q)
+        return self.alphaem_fixed_flavor(
+            q, self.alphaem_thresh[(nf, nl)], self.thresh[(nf, nl)], nf, nl
+        )
+
+    def alphaem_fixed_flavor_trn(self, q, alphaem_ref, qref, nf, nl):
+        """
+        Compute the running alphaem for nf fixed at, at least, NLO, using truncated method.
+        In this case the RGE for alpha_em is solved decoupling it from the RGE for alpha_s
+        (so the mixed terms are removed). alpha_s will just be unused.
+
+        Parameters
+        ----------
+        q : float
+            target scale
+        alph_aem_ref : float
+            reference value of alpha_em
+        qref: float
+            reference scale
+        nf: int
+            number of flavors
+        nl: int
+            number of leptons
+
+        Returns
+        -------
+        alpha_em at NLO : float
+            target value of a
+        """
+        if (nf, nl) == (0, 0):
+            return alphaem_ref
+        lmu = 2 * np.log(q / qref)
+        den = 1.0 + self.betas_qed[(nf, nl)][0] * alphaem_ref * lmu
+        alpha_LO = alphaem_ref / den
+        alpha_NLO = alpha_LO * (
+            1 - self.betas_qed[(nf, nl)][1] / self.betas_qed[(nf, nl)][0] * alpha_LO * np.log(den)
+        )
+        return alpha_NLO
+
+    def alphaem_fixed_flavor_exa(self, q, alphaem_ref, qref, nf, nl):
+        """
+        Compute numerically the running alphaem for nf fixed.
+
+        Parameters
+        ----------
+        q : float
+            target scale
+        alph_aem_ref : float
+            reference value of alpha_em
+        qref: float
+            reference scale
+        nf: int
+            number of flavors
+        nl: int
+            number of leptons
+
+        Returns
+        -------
+        alpha_em: float
+            target value of a
+        """
+        # at LO in QED the TRN solution is exact
+        if self.qed_order == 1:
+            return self.alphaem_fixed_flavor_trn(q, alphaem_ref, qref, nf, nl)
+        
+        u = 2 * np.log(q / qref)
+
+        # solve RGE
+        res = solve_ivp(
+            _rge, (0, u), (alphaem_ref,), args=[self.betas_qed[(nf, nl)]], method="Radau", rtol=1e-6
+        )
+        # taking fist (and only) element of y since it is a 1-D differential equation.
+        #  then the last element of the array which corresponds to alpha(q)
+        return res.y[0][-1]
+
+    def compute_alphaem_at_thresholds(self):
+        """
+        Compute and store alphaem at thresholds to speed up the calling
+        to alpha_em inside fiatlux:
+        when q is in a certain range (e.g. thresh_c < q < thresh_b) and qref in a different one
+        (e.g. thresh_b < q < thresh_t) we need to evolve from qref to thresh_b with nf=5 and then
+        from thresh_b to q with nf=4. Given that the value of alpha at thresh_b is always the same
+        we can avoid computing the first step storing the values of alpha in the threshold points.
+        It is done for qref in a generic range (not necessarly qref=91.2).
+
+        """
+        # determine nfref
+        nfref, nlref = self.find_region(self.qref)
+
+        thresh_list = [ME, MMU, MQL, self.thresh_c, MTAU, self.thresh_b, self.thresh_t]
+        thresh_list.append(self.qref)
+        thresh_list.sort()
+
+        thresh = {}
+
+        for mq in thresh_list:
+            eps = -1e-6 if mq < self.qref else 1e-6
+            if np.isfinite(mq):
+                nf_, nl_ = self.find_region(mq + eps)
+                thresh[(nf_, nl_)] = mq
+
+        regions = list(thresh.keys())
+
+        self.betas_qed = self.compute_betas(regions)
+
+        start = regions.index((nfref, nlref))
+
+        alphaem_thresh = {(nfref, nlref): self.alpha_em_ref}
+
+        for i in range(start + 1, len(regions)):
+            alphaem_thresh[regions[i]] = self.alphaem_fixed_flavor(
+                thresh[regions[i]],
+                alphaem_thresh[regions[i - 1]],
+                thresh[regions[i - 1]],
+                *regions[i - 1],
+            )
+
+        for i in reversed(range(0, start)):
+            alphaem_thresh[regions[i]] = self.alphaem_fixed_flavor(
+                thresh[regions[i]],
+                alphaem_thresh[regions[i + 1]],
+                thresh[regions[i + 1]],
+                *regions[i + 1],
+            )
+
+        return thresh, alphaem_thresh
+
+    def compute_betas(self, regions):
+        """Set values of betaQCD and betaQED."""
+        betas_qed = {}
+        for nf, nl in regions:
+            betas_qed[(nf, nl)] = [
+                beta.beta_qed_aem2(nf, nl) / (4 * np.pi),
+                beta.beta_qed_aem3(nf, nl) / (4 * np.pi) ** 2 if self.theory['QED'] == 2 else 0.,
+            ]
+        return betas_qed
+
+    def find_region(self, q):
+        if q < ME:
+            nf = 0
+            nl = 0
+        elif q < MMU:
+            nf = 0
+            nl = 1
+        elif q < MQL:
+            nf = 0
+            nl = 2
+        elif q < self.thresh_c:
+            nf = 3
+            nl = 2
+        elif q < MTAU:
+            nf = 4
+            nl = 2
+        elif q < self.thresh_b:
+            nf = 4
+            nl = 3
+        elif q < self.thresh_t:
+            nf = 5
+            nl = 3
+        else:
+            nf = 6
+            nl = 3
+        return nf, nl
+
+
+def _rge(_t, alpha_t, beta_qed_vec):
+    """RGEs for the running of alphaem"""
+    rge_qed = -(alpha_t**2) * (
+        beta_qed_vec[0] + np.sum([alpha_t ** (k + 1) * beta for k, beta in enumerate(beta_qed_vec[1:])])
+    )
+    return rge_qed
diff --git a/validphys2/src/validphys/photon/compute.py b/validphys2/src/validphys/photon/compute.py
index 56872540be..453f45c8ba 100644
--- a/validphys2/src/validphys/photon/compute.py
+++ b/validphys2/src/validphys/photon/compute.py
@@ -1,19 +1,19 @@
-"""Script that calls fiatlux to add the photon PDF."""
+"""Module that calls fiatlux to add the photon PDF."""
 import logging
 import tempfile
 
 import fiatlux
 import numpy as np
-from scipy.integrate import solve_ivp, trapezoid
+from scipy.integrate import trapezoid
 from scipy.interpolate import interp1d
 import yaml
 
-from eko import beta
 from eko.io import EKO
 from n3fit.io.writer import XGRID
 from validphys.n3fit_data import replica_luxseed
 
 from . import structure_functions as sf
+from .alpha import Alpha
 
 log = logging.getLogger(__name__)
 
@@ -225,202 +225,3 @@ def generate_errors(self, replica):
         u, s, _ = np.linalg.svd(self.error_matrix, full_matrices=False)
         errors = u @ (s * rng.normal(size=7))
         return errors
-
-
-class Alpha:
-    def __init__(self, theory, q2max):
-        self.theory = theory
-        self.alpha_em_ref = theory["alphaqed"]
-        self.qref = self.theory["Qref"]
-        self.betas_qed = self.compute_betas()
-
-        # compute and store thresholds
-        self.thresh_c = self.theory["kcThr"] * self.theory["mc"]
-        self.thresh_b = self.theory["kbThr"] * self.theory["mb"]
-        self.thresh_t = self.theory["ktThr"] * self.theory["mt"]
-        if self.theory["MaxNfAs"] <= 5:
-            self.thresh_t = np.inf
-        if self.theory["MaxNfAs"] <= 4:
-            self.thresh_b = np.inf
-        if self.theory["MaxNfAs"] <= 3:
-            self.thresh_c = np.inf
-
-        if self.theory["ModEv"] == "TRN":
-            self.alphaem_fixed_flavor = self.alphaem_fixed_flavor_trn
-            self.thresh, self.alphaem_thresh = self.compute_alphaem_at_thresholds()
-        elif self.theory["ModEv"] == "EXA":
-            self.alphaem_fixed_flavor = self.alphaem_fixed_flavor_exa
-            self.thresh, self.alphaem_thresh = self.compute_alphaem_at_thresholds()
-
-            xmin = XGRID[0]
-            qmin = xmin * theory["MP"] / np.sqrt(1 - xmin)
-            # use a lot of interpolation points since it is a long path 1e-9 -> 1e4
-            self.q = np.geomspace(qmin, np.sqrt(q2max), 500, endpoint=True)
-
-            # add threshold points in the q list since alpha is not smooth there
-            self.q = np.append(self.q, [self.thresh_c, self.thresh_b, self.thresh_t])
-            self.q = self.q[np.isfinite(self.q)]
-            self.q.sort()
-
-            self.alpha_vec = np.array([self.alpha_em(q_) for q_ in self.q])
-            self.alpha_em = self.interpolate_alphaem
-        else:
-            raise ValueError(f"Evolution mode not recognized: {self.theory['ModEv']}")
-
-    def interpolate_alphaem(self, q):
-        r"""
-        Interpolate precomputed values of alpha_em.
-
-        Parameters
-        ----------
-        q: float
-            value in which alpha_em is computed
-
-        Returns
-        -------
-        alpha_em: float
-            electromagnetic coupling
-        """
-        return np.interp(q, self.q, self.alpha_vec)
-
-    def alpha_em(self, q):
-        r"""
-        Compute the value of the running alphaem.
-
-        Parameters
-        ----------
-        q: float
-            value in which alphaem is computed
-
-        Returns
-        -------
-        alpha_em: numpy.ndarray
-            electromagnetic coupling
-        """
-        if q < self.thresh_c:
-            nf = 3
-        elif q < self.thresh_b:
-            nf = 4
-        elif q < self.thresh_t:
-            nf = 5
-        else:
-            nf = 6
-        return self.alphaem_fixed_flavor(q, self.alphaem_thresh[nf], self.thresh[nf], nf)
-
-    def alphaem_fixed_flavor_trn(self, q, alphaem_ref, qref, nf):
-        """
-        Compute the running alphaem for nf fixed at NLO, using truncated method.
-        In this case the RGE for alpha_em is solved decoupling it from the RGE for alpha_s
-        (so the mixed terms are removed). alpha_s will just be unused.
-
-        Parameters
-        ----------
-        q : float
-            target scale
-        alph_aem_ref : float
-            reference value of alpha_em
-        qref: float
-            reference scale
-        nf: int
-            number of flavors
-
-        Returns
-        -------
-        alpha_em at NLO : float
-            target value of a
-        """
-        alpha_ref = alphaem_ref
-        lmu = 2 * np.log(q / qref)
-        den = 1.0 + self.betas_qed[nf][0] * alpha_ref * lmu
-        alpha_LO = alpha_ref / den
-        alpha_NLO = alpha_LO * (1 - self.betas_qed[nf][1] * alpha_LO * np.log(den))
-        return alpha_NLO
-
-    def alphaem_fixed_flavor_exa(self, q, alphaem_ref, qref, nf):
-        """
-        Compute numerically the running alphaem for nf fixed.
-
-        Parameters
-        ----------
-        q : float
-            target scale
-        alph_aem_ref : float
-            reference value of alpha_em
-        qref: float
-            reference scale
-        nf: int
-            number of flavors
-
-        Returns
-        -------
-        alpha_em: float
-            target value of a
-        """
-        u = 2 * np.log(q / qref)
-
-        # solve RGE
-        res = solve_ivp(
-            rge, (0, u), (alphaem_ref,), args=[self.betas_qed[nf]], method="Radau", rtol=1e-6
-        )
-        return res.y[0][-1]
-
-    def compute_alphaem_at_thresholds(self):
-        """
-        Compute and store alphaem at thresholds to speed up the calling
-        to alpha_em inside fiatlux:
-        when q is in a certain range (e.g. thresh_c < q < thresh_b) and qref in a different one
-        (e.g. thresh_b < q < thresh_t) we need to evolve from qref to thresh_b with nf=5 and then
-        from thresh_b to q with nf=4. Given that the value of alpha at thresh_b is always the same
-        we can avoid computing the first step storing the values of alpha in the threshold points.
-        It is done for qref in a generic range (not necessarly qref=91.2).
-
-        """
-        # determine nfref
-        if self.qref < self.thresh_c:
-            nfref = 3
-        elif self.qref < self.thresh_b:
-            nfref = 4
-        elif self.qref < self.thresh_t:
-            nfref = 5
-        else:
-            nfref = 6
-
-        thresh_list = [self.thresh_c, self.thresh_b, self.thresh_t]
-        thresh_list.insert(nfref - 3, self.qref)
-
-        thresh = {nf: thresh_list[nf - 3] for nf in range(3, self.theory["MaxNfAs"] + 1)}
-
-        alphaem_thresh = {nfref: self.alpha_em_ref}
-
-        # determine the values of alphaem in the threshold points, depending on the value of qref
-        for nf in range(nfref + 1, self.theory["MaxNfAs"] + 1):
-            alphaem_thresh[nf] = self.alphaem_fixed_flavor(
-                thresh[nf], alphaem_thresh[nf - 1], thresh[nf - 1], nf - 1
-            )
-
-        for nf in reversed(range(3, nfref)):
-            alphaem_thresh[nf] = self.alphaem_fixed_flavor(
-                thresh[nf], alphaem_thresh[nf + 1], thresh[nf + 1], nf + 1
-            )
-
-        return thresh, alphaem_thresh
-
-    def compute_betas(self):
-        """Set values of betaQCD and betaQED."""
-        betas_qed = {}
-        for nf in range(3, 6 + 1):
-            vec_qed = [beta.beta_qed_aem2(nf) / (4 * np.pi)]
-            for ord in range(1, self.theory['QED']):
-                vec_qed.append(beta.b_qed((0, ord + 2), nf) / (4 * np.pi) ** ord)
-            betas_qed[nf] = vec_qed
-        return betas_qed
-
-
-def rge(_t, alpha, beta_qed_vec):
-    """RGEs for the running of alphaem"""
-    rge_qed = (
-        -(alpha**2)
-        * beta_qed_vec[0]
-        * (1 + np.sum([alpha ** (k + 1) * b for k, b in enumerate(beta_qed_vec[1:])]))
-    )
-    return rge_qed
diff --git a/validphys2/src/validphys/photon/constants.py b/validphys2/src/validphys/photon/constants.py
index d820fb12a9..3c536736cc 100644
--- a/validphys2/src/validphys/photon/constants.py
+++ b/validphys2/src/validphys/photon/constants.py
@@ -1,4 +1,8 @@
 EU2 = 4.0 / 9
 ED2 = 1.0 / 9
-NL = 3
 NC = 3
+
+MTAU = 1.777
+MMU = 0.1057
+ME = 0.0005110
+MQL = 1.0
diff --git a/validphys2/src/validphys/tests/photon/test_alpha.py b/validphys2/src/validphys/tests/photon/test_alpha.py
new file mode 100644
index 0000000000..7db54055b6
--- /dev/null
+++ b/validphys2/src/validphys/tests/photon/test_alpha.py
@@ -0,0 +1,243 @@
+import numpy as np
+
+from eko import beta
+from eko.couplings import Couplings, expanded_qed
+from eko.quantities.couplings import CouplingEvolutionMethod, CouplingsInfo
+from eko.quantities.heavy_quarks import QuarkMassScheme
+from validphys.api import API
+from validphys.photon import constants
+from validphys.photon.alpha import Alpha
+
+from ..conftest import THEORY_QED
+
+
+def test_masses_init():
+    """test thresholds values in Alpha class"""
+    test_theory = API.theoryid(theoryid=THEORY_QED)
+    theory = test_theory.get_description()
+    alpha = Alpha(theory, 1e8)
+    np.testing.assert_equal(alpha.thresh_t, np.inf)
+    np.testing.assert_almost_equal(alpha.thresh_b, theory["mb"])
+    np.testing.assert_almost_equal(alpha.thresh_c, theory["mc"])
+
+
+def test_set_thresholds_alpha_em():
+    """test value of alpha_em at threshold values"""
+    test_theory = API.theoryid(theoryid=THEORY_QED)
+    theory = test_theory.get_description()
+
+    for modev in ['EXA', 'TRN']:
+        theory['ModEv'] = modev
+
+        alpha = Alpha(theory, 1e8)
+        alpha_ref = theory["alphaqed"]
+
+        # test all the thresholds
+        np.testing.assert_almost_equal(alpha.alpha_em_ref, theory["alphaqed"])
+        np.testing.assert_almost_equal(alpha.thresh[(5, 3)], theory["Qedref"])
+        np.testing.assert_almost_equal(alpha.thresh[(4, 3)], theory["mb"])
+        np.testing.assert_almost_equal(alpha.thresh[(4, 2)], constants.MTAU)
+        np.testing.assert_almost_equal(alpha.thresh[(3, 2)], theory["mc"])
+        np.testing.assert_almost_equal(alpha.thresh[(0, 2)], constants.MQL)
+        np.testing.assert_almost_equal(alpha.thresh[(0, 1)], constants.MMU)
+        np.testing.assert_almost_equal(alpha.thresh[(0, 0)], constants.ME)
+
+        # test some values of alpha at the threshold points
+        np.testing.assert_almost_equal(alpha.alphaem_thresh[(5, 3)], theory["alphaqed"])
+        np.testing.assert_almost_equal(
+            alpha.alphaem_thresh[(4, 3)],
+            alpha.alphaem_fixed_flavor(theory["mb"], alpha_ref, theory["Qedref"], 5, 3),
+        )
+        np.testing.assert_almost_equal(
+            alpha.alphaem_thresh[(4, 2)],
+            alpha.alphaem_fixed_flavor(
+                constants.MTAU, alpha.alphaem_thresh[(4, 3)], theory["mb"], 4, 3
+            ),
+        )
+        np.testing.assert_almost_equal(
+            alpha.alphaem_thresh[(3, 2)],
+            alpha.alphaem_fixed_flavor(
+                theory["mc"], alpha.alphaem_thresh[(4, 2)], constants.MTAU, 4, 2
+            ),
+        )
+        np.testing.assert_almost_equal(
+            alpha.alphaem_thresh[(0, 2)],
+            alpha.alphaem_fixed_flavor(
+                constants.MQL, alpha.alphaem_thresh[(3, 2)], theory["mc"], 3, 2
+            ),
+        )
+        np.testing.assert_almost_equal(
+            alpha.alphaem_thresh[(0, 1)],
+            alpha.alphaem_fixed_flavor(
+                constants.MMU, alpha.alphaem_thresh[(0, 2)], constants.MQL, 0, 2
+            ),
+        )
+        np.testing.assert_almost_equal(
+            alpha.alphaem_thresh[(0, 0)],
+            alpha.alphaem_fixed_flavor(
+                constants.ME, alpha.alphaem_thresh[(0, 1)], constants.MMU, 0, 1
+            ),
+        )
+
+        np.testing.assert_equal(len(alpha.alphaem_thresh), 7)
+        np.testing.assert_equal(len(alpha.thresh), 7)
+
+
+def test_couplings_exa():
+    """
+    Benchmark the exact running of alpha:
+    Testing both the FFS and the VFS with the results from EKO.
+    Test also the values of the couplings at the threshold values.
+    """
+    test_theory = API.theoryid(theoryid=THEORY_QED)
+    theory = test_theory.get_description()
+    for k in [1]:
+        theory["mb"] *= k
+        mass_list = [theory["mc"], theory["mb"], theory["mt"]]
+
+        alpha = Alpha(theory, 1e8)
+        couplings = CouplingsInfo.from_dict(
+            dict(
+                alphas=theory["alphas"],
+                alphaem=theory["alphaqed"],
+                scale=theory["Qref"],
+                num_flavs_ref=None,
+                max_num_flavs=theory["MaxNfAs"],
+                em_running=True,
+            )
+        )
+        eko_alpha = Couplings(
+            couplings,
+            (theory["PTO"] + 1, theory["QED"]),
+            method=CouplingEvolutionMethod.EXACT,
+            masses=[m**2 for m in mass_list],
+            hqm_scheme=QuarkMassScheme.POLE,
+            thresholds_ratios=[1.0, 1.0, 1.0],
+        )
+        eko_alpha.decoupled_running = True
+        alpha_ref = theory["alphaqed"]
+        for q in [5, 10, 20, 50, 80, 100, 200]:
+            np.testing.assert_allclose(
+                alpha.alphaem_fixed_flavor(q, alpha_ref, theory["Qref"], 5, 3),
+                alpha.alpha_em(q),
+                rtol=5e-6,
+            )
+            np.testing.assert_allclose(
+                alpha.alphaem_fixed_flavor(q, alpha_ref, theory["Qref"], 5, 3),
+                eko_alpha.compute_exact_alphaem_running(
+                    np.array([0.118, alpha_ref]) / (4 * np.pi), 5, 3, theory["Qref"] ** 2, q**2
+                )[1]
+                * 4
+                * np.pi,
+                rtol=1e-7,
+            )
+
+        # TODO: these tests have to be switched on once the varying nl is implemented in eko
+
+        for q in [1, 2, 3, 4]:
+            np.testing.assert_allclose(
+                alpha.alpha_em(q), eko_alpha.a_em(q**2) * 4 * np.pi, rtol=5e-6
+            )
+        for nf in range(3, theory["MaxNfAs"]):
+            np.testing.assert_allclose(
+                alpha.alphaem_thresh[(nf, 2 if nf == 3 else 3)],
+                eko_alpha.a_em(mass_list[nf - 3] ** 2, nf) * 4 * np.pi,
+                rtol=3e-7,
+            )
+
+
+def test_exa_interpolation():
+    """test the accuracy of the alphaem interpolation"""
+    test_theory = API.theoryid(theoryid=THEORY_QED)
+    theory = test_theory.get_description()
+
+    alpha = Alpha(theory, 1e8)
+    for q in np.geomspace(5.0, 1e4, 1000, endpoint=True):
+        np.testing.assert_allclose(
+            alpha.alpha_em(q),
+            alpha.alphaem_fixed_flavor(q, theory["alphaqed"], theory["Qref"], 5, 3),
+            rtol=1e-5,
+        )
+    for q in np.geomspace(1.778, 4.9, 20):
+        np.testing.assert_allclose(
+            alpha.alpha_em(q),
+            alpha.alphaem_fixed_flavor(q, alpha.alphaem_thresh[(4, 3)], alpha.thresh_b, 4, 3),
+            rtol=1e-5,
+        )
+    for q in np.geomspace(0.2, 0.8, 20):
+        np.testing.assert_allclose(
+            alpha.alpha_em(q),
+            alpha.alphaem_fixed_flavor(q, alpha.alphaem_thresh[(0, 2)], constants.MQL, 0, 2),
+            rtol=1e-5,
+        )
+
+
+def test_couplings_trn():
+    """benchmark the truncated running of alpha with EKO"""
+    test_theory = API.theoryid(theoryid=THEORY_QED)
+    theory = test_theory.get_description()
+    theory["ModEv"] = 'TRN'
+    alpha = Alpha(theory, 1e8)
+    alpha_ref = theory["alphaqed"]
+
+    for q in [80, 10, 5]:
+        np.testing.assert_allclose(
+            alpha.alphaem_fixed_flavor(q, alpha_ref, theory["Qref"], 5, 3),
+            alpha.alpha_em(q),
+            rtol=1e-10,
+        )
+        np.testing.assert_allclose(
+            alpha.alphaem_fixed_flavor(q, alpha_ref, theory["Qref"], 5, 3),
+            expanded_qed(
+                alpha_ref / (4 * np.pi),
+                theory["QED"],
+                beta.beta_qed((0, 2), 5, 3),
+                [beta.b_qed((0, i + 2), 5, 3) for i in range(theory["QED"])],
+                2 * np.log(q / theory["Qref"]),
+            )
+            * 4
+            * np.pi,
+            rtol=1e-7,
+        )
+
+
+def test_trn_vs_exa():
+    """benchmark the exact running of alpha vs truncated"""
+    # does this test make sense?
+    test_theory = API.theoryid(theoryid=THEORY_QED)
+    theory_exa = test_theory.get_description()
+    theory_trn = theory_exa.copy()
+    theory_trn["ModEv"] = 'TRN'
+
+    alpha_exa = Alpha(theory_exa, 1e8)
+    alpha_trn = Alpha(theory_trn, 1e8)
+
+    for q in [1, 3, 10, 50, 80]:
+        np.testing.assert_allclose(alpha_exa.alpha_em(q), alpha_trn.alpha_em(q), rtol=2e-3)
+
+
+def test_betas():
+    """test betas for different nf"""
+    test_theory = API.theoryid(theoryid=THEORY_QED)
+    alpha = Alpha(test_theory.get_description(), 1e8)
+    vec_beta0 = [
+        -0.0,
+        -0.10610329539459688,
+        -0.21220659078919377,
+        -0.42441318157838753,
+        -0.5658842421045168,
+        -0.6719875374991137,
+        -0.7073553026306458,
+    ]
+    vec_beta1 = [
+        -0.0,
+        -0.025330295910584444,
+        -0.05066059182116889,
+        -0.06754745576155852,
+        -0.0825580014863493,
+        -0.10788829739693376,
+        -0.10882645650473316,
+    ]
+    for i, (nf, nl) in enumerate([(0, 0), (0, 1), (0, 2), (3, 2), (4, 2), (4, 3), (5, 3)]):
+        np.testing.assert_allclose(alpha.betas_qed[(nf, nl)][0], vec_beta0[i], rtol=1e-7)
+        np.testing.assert_allclose(alpha.betas_qed[(nf, nl)][1], vec_beta1[i], rtol=1e-7)
diff --git a/validphys2/src/validphys/tests/photon/test_compute.py b/validphys2/src/validphys/tests/photon/test_compute.py
index 0d4c7ea57d..0cd51c561c 100644
--- a/validphys2/src/validphys/tests/photon/test_compute.py
+++ b/validphys2/src/validphys/tests/photon/test_compute.py
@@ -4,17 +4,14 @@
 import numpy as np
 import yaml
 
-from eko import beta
-from eko.couplings import Couplings, expanded_qed
 from eko.io import EKO
-from eko.quantities.couplings import CouplingEvolutionMethod, CouplingsInfo
-from eko.quantities.heavy_quarks import QuarkMassScheme
 from n3fit.io.writer import XGRID
 from validphys.api import API
 from validphys.core import PDF as PDFset
 from validphys.loader import FallbackLoader
 from validphys.photon import structure_functions as sf
-from validphys.photon.compute import FIATLUX_DEFAULT, Alpha, Photon
+from validphys.photon.compute import FIATLUX_DEFAULT, Photon
+from validphys.photon.alpha import Alpha
 
 from ..conftest import PDF, THEORY_QED
 
@@ -47,169 +44,6 @@ def test_parameters_init():
     np.testing.assert_equal(photon.q_in, 100.0)
 
 
-def test_masses_init():
-    """test thresholds values in Alpha class"""
-    test_theory = API.theoryid(theoryid=THEORY_QED)
-    theory = test_theory.get_description()
-    alpha = Alpha(theory, 1e8)
-    np.testing.assert_equal(alpha.thresh_t, np.inf)
-    np.testing.assert_almost_equal(alpha.thresh_b, theory["mb"])
-    np.testing.assert_almost_equal(alpha.thresh_c, theory["mc"])
-
-
-def test_set_thresholds_alpha_em():
-    """test value of alpha_em at threshold values"""
-    test_theory = API.theoryid(theoryid=THEORY_QED)
-    theory = test_theory.get_description()
-
-    for modev in ['EXA', 'TRN']:
-        theory['ModEv'] = modev
-
-        alpha = Alpha(theory, 1e8)
-        alpha_ref = theory["alphaqed"]
-
-        np.testing.assert_almost_equal(alpha.alpha_em_ref, theory["alphaqed"])
-        np.testing.assert_almost_equal(alpha.thresh[5], theory["Qedref"])
-        np.testing.assert_almost_equal(alpha.thresh[4], theory["mb"])
-        np.testing.assert_almost_equal(alpha.thresh[3], theory["mc"])
-        np.testing.assert_almost_equal(alpha.alphaem_thresh[5], theory["alphaqed"])
-        np.testing.assert_almost_equal(
-            alpha.alphaem_thresh[4],
-            alpha.alphaem_fixed_flavor(theory["mb"], alpha_ref, theory["Qedref"], 5),
-        )
-        np.testing.assert_almost_equal(
-            alpha.alphaem_thresh[3],
-            alpha.alphaem_fixed_flavor(theory["mc"], alpha.alphaem_thresh[4], theory["mb"], 4),
-        )
-        np.testing.assert_equal(len(alpha.alphaem_thresh), 3)
-        np.testing.assert_equal(len(alpha.thresh), 3)
-
-
-def test_couplings_exa():
-    """
-    Benchmark the exact running of alpha:
-    Testing both the FFS and the VFS with the results from EKO.
-    Test also the values of the couplings at the threshold values.
-    """
-    test_theory = API.theoryid(theoryid=THEORY_QED)
-    theory = test_theory.get_description()
-    for k in [0.5, 1, 2]:
-        theory["mb"] *= k
-        mass_list = [theory["mc"], theory["mb"], theory["mt"]]
-
-        alpha = Alpha(theory, 1e8)
-        couplings = CouplingsInfo.from_dict(
-            dict(
-                alphas=theory["alphas"],
-                alphaem=theory["alphaqed"],
-                scale=theory["Qref"],
-                num_flavs_ref=None,
-                max_num_flavs=theory["MaxNfAs"],
-                em_running=True,
-            )
-        )
-        eko_alpha = Couplings(
-            couplings,
-            (theory["PTO"] + 1, theory["QED"]),
-            method=CouplingEvolutionMethod.EXACT,
-            masses=[m**2 for m in mass_list],
-            hqm_scheme=QuarkMassScheme.POLE,
-            thresholds_ratios=[1.0, 1.0, 1.0],
-        )
-        eko_alpha.decoupled_running = True
-        alpha_ref = theory["alphaqed"]
-        for q in [5, 10, 20, 50, 80, 100, 200]:
-            np.testing.assert_allclose(
-                alpha.alphaem_fixed_flavor(q, alpha_ref, theory["Qref"], 5),
-                alpha.alpha_em(q),
-                rtol=5e-6,
-            )
-            np.testing.assert_allclose(
-                alpha.alphaem_fixed_flavor(q, alpha_ref, theory["Qref"], 5),
-                eko_alpha.compute_exact_alphaem_running(
-                    np.array([0.118, alpha_ref]) / (4 * np.pi), 5, theory["Qref"] ** 2, q**2
-                )[1]
-                * 4
-                * np.pi,
-                rtol=1e-7,
-            )
-        for q in [1, 2, 3, 4]:
-            np.testing.assert_allclose(
-                alpha.alpha_em(q), eko_alpha.a_em(q**2) * 4 * np.pi, rtol=5e-6
-            )
-        for nf in range(3, theory["MaxNfAs"]):
-            np.testing.assert_allclose(
-                alpha.alphaem_thresh[nf],
-                eko_alpha.a_em(mass_list[nf - 3] ** 2, nf) * 4 * np.pi,
-                rtol=3e-7,
-            )
-
-
-def test_exa_interpolation():
-    """test the accuracy of the alphaem interpolation"""
-    test_theory = API.theoryid(theoryid=THEORY_QED)
-    theory = test_theory.get_description()
-
-    alpha = Alpha(theory, 1e8)
-    for q in np.geomspace(1.0, 1e4, 1000, endpoint=True):
-        np.testing.assert_allclose(alpha.alpha_em(q), alpha.alpha_em(q), rtol=1e-5)
-
-
-def test_couplings_trn():
-    """benchmark the truncated running of alpha with EKO"""
-    test_theory = API.theoryid(theoryid=THEORY_QED)
-    theory = test_theory.get_description()
-    theory["ModEv"] = 'TRN'
-    alpha = Alpha(theory, 1e8)
-    alpha_ref = theory["alphaqed"]
-
-    for q in [80, 10, 5]:
-        np.testing.assert_allclose(
-            alpha.alphaem_fixed_flavor(q, alpha_ref, theory["Qref"], 5),
-            alpha.alpha_em(q),
-            rtol=1e-10,
-        )
-        np.testing.assert_allclose(
-            alpha.alphaem_fixed_flavor(q, alpha_ref, theory["Qref"], 5),
-            expanded_qed(
-                alpha_ref / (4 * np.pi),
-                theory["QED"],
-                beta.beta_qed((0, 2), 5),
-                [beta.b_qed((0, i + 2), 5) for i in range(theory["QED"])],
-                2 * np.log(q / theory["Qref"]),
-            )
-            * 4
-            * np.pi,
-            rtol=1e-7,
-        )
-
-
-def test_trn_vs_exa():
-    """benchmark the exact running of alpha vs truncated"""
-    # does this test make sense?
-    test_theory = API.theoryid(theoryid=THEORY_QED)
-    theory_exa = test_theory.get_description()
-    theory_trn = theory_exa.copy()
-    theory_trn["ModEv"] = 'TRN'
-
-    alpha_exa = Alpha(theory_exa, 1e8)
-    alpha_trn = Alpha(theory_trn, 1e8)
-
-    for q in [1, 3, 10, 50, 80]:
-        np.testing.assert_allclose(alpha_exa.alpha_em(q), alpha_trn.alpha_em(q), rtol=2e-3)
-
-
-def test_betas():
-    """test betas for different nf"""
-    test_theory = API.theoryid(theoryid=THEORY_QED)
-    alpha = Alpha(test_theory.get_description(), 1e8)
-    vec_beta0 = [-0.5305164769729844, -0.6719875374991137, -0.7073553026306458, -0.8488263631567751]
-    vec_b1 = [0.17507043740108488, 0.1605510390839295, 0.1538497783221655, 0.1458920311675707]
-    for nf in range(3, 6 + 1):
-        np.testing.assert_allclose(alpha.betas_qed[nf][0], vec_beta0[nf - 3], rtol=1e-7)
-        np.testing.assert_allclose(alpha.betas_qed[nf][1], vec_b1[nf - 3], rtol=1e-7)
-
-
 def test_photon():
     """
     Test that photon coming out of Photon interpolator matches the photon array