diff --git a/nnpdfcpp/data/theory.db b/nnpdfcpp/data/theory.db
index 26df4c3830..bf1cae6622 100644
Binary files a/nnpdfcpp/data/theory.db and b/nnpdfcpp/data/theory.db differ
diff --git a/validphys2/src/validphys/config.py b/validphys2/src/validphys/config.py
index e94a0f8cd4..e54981d5c4 100644
--- a/validphys2/src/validphys/config.py
+++ b/validphys2/src/validphys/config.py
@@ -1404,7 +1404,7 @@ def produce_defaults(
         default_filter_settings_recorded_spec_=None,
     ):
         """Produce default values for filters taking into account the
-        values of ``q2min``, ` `w2min`` and ``maxTau`` defined at namespace
+        values of ``q2min``, ``w2min`` and ``maxTau`` defined at namespace
         level and those inside a ``filter_defaults`` mapping.
         """
         from validphys.filters import default_filter_settings_input
@@ -1416,7 +1416,7 @@ def produce_defaults(
             raise ConfigError("q2min defined multiple times with different values")
         if w2min is not None and "w2min" in filter_defaults and w2min != filter_defaults["w2min"]:
             raise ConfigError("w2min defined multiple times with different values")
-        
+
         if (
             maxTau is not None
             and "maxTau" in filter_defaults
@@ -1442,7 +1442,7 @@ def produce_defaults(
         if w2min is not None and defaults_loaded:
             log.warning("Using w2min from runcard")
             filter_defaults["w2min"] = w2min
-            
+
         if maxTau is not None and defaults_loaded:
             log.warning("Using maxTau from runcard")
             filter_defaults["maxTau"] = maxTau
diff --git a/validphys2/src/validphys/scalevariations/pointprescriptions.yaml b/validphys2/src/validphys/scalevariations/pointprescriptions.yaml
index bffeb0e35f..cfe7a83f8d 100644
--- a/validphys2/src/validphys/scalevariations/pointprescriptions.yaml
+++ b/validphys2/src/validphys/scalevariations/pointprescriptions.yaml
@@ -9,3 +9,15 @@
 '7 point': ['(1, 1)', '(2, 1)', '(0.5, 1)', '(1, 2)', '(1, 0.5)', '(2, 2)', '(0.5, 0.5)']
 '7original point': ['(1, 1)', '(2, 1)', '(0.5, 1)', '(1, 2)', '(1, 0.5)', '(2, 2)', '(0.5, 0.5)']
 '9 point': ['(1, 1)', '(2, 1)', '(0.5, 1)', '(1, 2)', '(1, 0.5)', '(2, 2)', '(0.5, 0.5)', '(2, 0.5)', '(0.5, 2)']
+# N3LO IHOU Anomalous dimension theory covmat
+'n3lo ad covmat': ['(0, 0, 0, 0)','(1, 0, 0, 0)','(2, 0, 0, 0)','(3, 0, 0, 0)','(4, 0, 0, 0)','(5, 0, 0, 0)','(6, 0, 0, 0)','(7, 0, 0, 0)','(8, 0, 0, 0)','(9, 0, 0, 0)','(10, 0, 0, 0)','(11, 0, 0, 0)','(12, 0, 0, 0)','(13, 0, 0, 0)','(14, 0, 0, 0)','(15, 0, 0, 0)','(16, 0, 0, 0)','(17, 0, 0, 0)','(18, 0, 0, 0)','(19, 0, 0, 0)','(0, 1, 0, 0)','(0, 2, 0, 0)','(0, 3, 0, 0)','(0, 4, 0, 0)','(0, 5, 0, 0)','(0, 6, 0, 0)','(0, 7, 0, 0)','(0, 8, 0, 0)','(0, 9, 0, 0)','(0, 10, 0, 0)','(0, 11, 0, 0)','(0, 12, 0, 0)','(0, 13, 0, 0)','(0, 14, 0, 0)','(0, 15, 0, 0)','(0, 16, 0, 0)','(0, 17, 0, 0)','(0, 18, 0, 0)','(0, 19, 0, 0)','(0, 20, 0, 0)','(0, 21, 0, 0)','(0, 0, 1, 0)','(0, 0, 2, 0)','(0, 0, 3, 0)','(0, 0, 4, 0)','(0, 0, 5, 0)','(0, 0, 6, 0)','(0, 0, 7, 0)','(0, 0, 8, 0)','(0, 0, 9, 0)','(0, 0, 10, 0)','(0, 0, 11, 0)','(0, 0, 12, 0)','(0, 0, 13, 0)','(0, 0, 14, 0)','(0, 0, 15, 0)','(0, 0, 0, 1)','(0, 0, 0, 2)','(0, 0, 0, 3)','(0, 0, 0, 4)','(0, 0, 0, 5)','(0, 0, 0, 6)']
+# N3LO full IHOU: Anomalous dimension theory covmat + DIS massive coefficient functions
+'n3lo ihou': ['(0, 0, 0, 0)','(1, 0, 0, 0)','(2, 0, 0, 0)','(3, 0, 0, 0)','(4, 0, 0, 0)','(5, 0, 0, 0)','(6, 0, 0, 0)','(7, 0, 0, 0)','(8, 0, 0, 0)','(9, 0, 0, 0)','(10, 0, 0, 0)','(11, 0, 0, 0)','(12, 0, 0, 0)','(13, 0, 0, 0)','(14, 0, 0, 0)','(15, 0, 0, 0)','(16, 0, 0, 0)','(17, 0, 0, 0)','(18, 0, 0, 0)','(19, 0, 0, 0)','(0, 1, 0, 0)','(0, 2, 0, 0)','(0, 3, 0, 0)','(0, 4, 0, 0)','(0, 5, 0, 0)','(0, 6, 0, 0)','(0, 7, 0, 0)','(0, 8, 0, 0)','(0, 9, 0, 0)','(0, 10, 0, 0)','(0, 11, 0, 0)','(0, 12, 0, 0)','(0, 13, 0, 0)','(0, 14, 0, 0)','(0, 15, 0, 0)','(0, 16, 0, 0)','(0, 17, 0, 0)','(0, 18, 0, 0)','(0, 19, 0, 0)','(0, 20, 0, 0)','(0, 21, 0, 0)','(0, 0, 1, 0)','(0, 0, 2, 0)','(0, 0, 3, 0)','(0, 0, 4, 0)','(0, 0, 5, 0)','(0, 0, 6, 0)','(0, 0, 7, 0)','(0, 0, 8, 0)','(0, 0, 9, 0)','(0, 0, 10, 0)','(0, 0, 11, 0)','(0, 0, 12, 0)','(0, 0, 13, 0)','(0, 0, 14, 0)','(0, 0, 15, 0)','(0, 0, 0, 1)','(0, 0, 0, 2)','(0, 0, 0, 3)','(0, 0, 0, 4)','(0, 0, 0, 5)','(0, 0, 0, 6)','(-1, -1)','(1, 1)']
+# N3LO full IHOU + 7 point scale variations 
+'n3lo full thcovmat': ['(0, 0, 0, 0)','(1, 0, 0, 0)','(2, 0, 0, 0)','(3, 0, 0, 0)','(4, 0, 0, 0)','(5, 0, 0, 0)','(6, 0, 0, 0)','(7, 0, 0, 0)','(8, 0, 0, 0)','(9, 0, 0, 0)','(10, 0, 0, 0)','(11, 0, 0, 0)','(12, 0, 0, 0)','(13, 0, 0, 0)','(14, 0, 0, 0)','(15, 0, 0, 0)','(16, 0, 0, 0)','(17, 0, 0, 0)','(18, 0, 0, 0)','(19, 0, 0, 0)','(0, 1, 0, 0)','(0, 2, 0, 0)','(0, 3, 0, 0)','(0, 4, 0, 0)','(0, 5, 0, 0)','(0, 6, 0, 0)','(0, 7, 0, 0)','(0, 8, 0, 0)','(0, 9, 0, 0)','(0, 10, 0, 0)','(0, 11, 0, 0)','(0, 12, 0, 0)','(0, 13, 0, 0)','(0, 14, 0, 0)','(0, 15, 0, 0)','(0, 16, 0, 0)','(0, 17, 0, 0)','(0, 18, 0, 0)','(0, 19, 0, 0)','(0, 20, 0, 0)','(0, 21, 0, 0)','(0, 0, 1, 0)','(0, 0, 2, 0)','(0, 0, 3, 0)','(0, 0, 4, 0)','(0, 0, 5, 0)','(0, 0, 6, 0)','(0, 0, 7, 0)','(0, 0, 8, 0)','(0, 0, 9, 0)','(0, 0, 10, 0)','(0, 0, 11, 0)','(0, 0, 12, 0)','(0, 0, 13, 0)','(0, 0, 14, 0)','(0, 0, 15, 0)','(0, 0, 0, 1)','(0, 0, 0, 2)','(0, 0, 0, 3)','(0, 0, 0, 4)','(0, 0, 0, 5)','(0, 0, 0, 6)','(2, 1)', '(0.5, 1)', '(1, 2)', '(1, 0.5)', '(2, 2)', '(0.5, 0.5)','(-1, -1)','(1, 1)']
+# N3LO full IHOU + 3 point scale variations for datasets with no N3LO correcttions 
+'n3lo 3pt missing': ['(0, 0, 0, 0)','(1, 0, 0, 0)','(2, 0, 0, 0)','(3, 0, 0, 0)','(4, 0, 0, 0)','(5, 0, 0, 0)','(6, 0, 0, 0)','(7, 0, 0, 0)','(8, 0, 0, 0)','(9, 0, 0, 0)','(10, 0, 0, 0)','(11, 0, 0, 0)','(12, 0, 0, 0)','(13, 0, 0, 0)','(14, 0, 0, 0)','(15, 0, 0, 0)','(16, 0, 0, 0)','(17, 0, 0, 0)','(18, 0, 0, 0)','(19, 0, 0, 0)','(0, 1, 0, 0)','(0, 2, 0, 0)','(0, 3, 0, 0)','(0, 4, 0, 0)','(0, 5, 0, 0)','(0, 6, 0, 0)','(0, 7, 0, 0)','(0, 8, 0, 0)','(0, 9, 0, 0)','(0, 10, 0, 0)','(0, 11, 0, 0)','(0, 12, 0, 0)','(0, 13, 0, 0)','(0, 14, 0, 0)','(0, 15, 0, 0)','(0, 16, 0, 0)','(0, 17, 0, 0)','(0, 18, 0, 0)','(0, 19, 0, 0)','(0, 20, 0, 0)','(0, 21, 0, 0)','(0, 0, 1, 0)','(0, 0, 2, 0)','(0, 0, 3, 0)','(0, 0, 4, 0)','(0, 0, 5, 0)','(0, 0, 6, 0)','(0, 0, 7, 0)','(0, 0, 8, 0)','(0, 0, 9, 0)','(0, 0, 10, 0)','(0, 0, 11, 0)','(0, 0, 12, 0)','(0, 0, 13, 0)','(0, 0, 14, 0)','(0, 0, 15, 0)','(0, 0, 0, 1)','(0, 0, 0, 2)','(0, 0, 0, 3)','(0, 0, 0, 4)','(0, 0, 0, 5)','(0, 0, 0, 6)', '(1, 0.5 missing)', '(1, 2 missing)','(-1, -1)','(1, 1)']
+# N3LO full IHOU + 3 point scale variations for hadronic dasasets
+'n3lo 3pt hadronic': ['(0, 0, 0, 0)','(1, 0, 0, 0)','(2, 0, 0, 0)','(3, 0, 0, 0)','(4, 0, 0, 0)','(5, 0, 0, 0)','(6, 0, 0, 0)','(7, 0, 0, 0)','(8, 0, 0, 0)','(9, 0, 0, 0)','(10, 0, 0, 0)','(11, 0, 0, 0)','(12, 0, 0, 0)','(13, 0, 0, 0)','(14, 0, 0, 0)','(15, 0, 0, 0)','(16, 0, 0, 0)','(17, 0, 0, 0)','(18, 0, 0, 0)','(19, 0, 0, 0)','(0, 1, 0, 0)','(0, 2, 0, 0)','(0, 3, 0, 0)','(0, 4, 0, 0)','(0, 5, 0, 0)','(0, 6, 0, 0)','(0, 7, 0, 0)','(0, 8, 0, 0)','(0, 9, 0, 0)','(0, 10, 0, 0)','(0, 11, 0, 0)','(0, 12, 0, 0)','(0, 13, 0, 0)','(0, 14, 0, 0)','(0, 15, 0, 0)','(0, 16, 0, 0)','(0, 17, 0, 0)','(0, 18, 0, 0)','(0, 19, 0, 0)','(0, 20, 0, 0)','(0, 21, 0, 0)','(0, 0, 1, 0)','(0, 0, 2, 0)','(0, 0, 3, 0)','(0, 0, 4, 0)','(0, 0, 5, 0)','(0, 0, 6, 0)','(0, 0, 7, 0)','(0, 0, 8, 0)','(0, 0, 9, 0)','(0, 0, 10, 0)','(0, 0, 11, 0)','(0, 0, 12, 0)','(0, 0, 13, 0)','(0, 0, 14, 0)','(0, 0, 15, 0)','(0, 0, 0, 1)','(0, 0, 0, 2)','(0, 0, 0, 3)','(0, 0, 0, 4)','(0, 0, 0, 5)','(0, 0, 0, 6)', '(1, 0.5 hadronic)', '(1, 2 hadronic)','(-1, -1)','(1, 1)']
+# N3LO 7 point scale variations
+'n3lo 7 point': ['(0, 0, 0, 0)', '(2, 1)', '(0.5, 1)', '(1, 2)', '(1, 0.5)', '(2, 2)', '(0.5, 0.5)']
diff --git a/validphys2/src/validphys/scalevariations/scalevariationtheoryids.yaml b/validphys2/src/validphys/scalevariations/scalevariationtheoryids.yaml
index 34e3ab0b71..6d1af1cff6 100644
--- a/validphys2/src/validphys/scalevariations/scalevariationtheoryids.yaml
+++ b/validphys2/src/validphys/scalevariations/scalevariationtheoryids.yaml
@@ -137,6 +137,87 @@ scale_variations_for:
       (2, 1): 437
       (2, 2): 438
 
+  # N3LO: MHOU + IHOU
+  - theoryid: 722
+    variations:
+      (0, 0, 0, 0): 722 # central th
+      (1, 0, 0, 0): 723 # anomalous dimensions variations
+      (2, 0, 0, 0): 724
+      (3, 0, 0, 0): 725
+      (4, 0, 0, 0): 726
+      (5, 0, 0, 0): 727
+      (6, 0, 0, 0): 728
+      (7, 0, 0, 0): 729
+      (8, 0, 0, 0): 730
+      (9, 0, 0, 0): 731
+      (10, 0, 0, 0): 732
+      (11, 0, 0, 0): 733
+      (12, 0, 0, 0): 734
+      (13, 0, 0, 0): 735
+      (14, 0, 0, 0): 736
+      (15, 0, 0, 0): 737
+      (16, 0, 0, 0): 738
+      (17, 0, 0, 0): 739
+      (18, 0, 0, 0): 740
+      (19, 0, 0, 0): 741
+      (0, 1, 0, 0): 742
+      (0, 2, 0, 0): 743
+      (0, 3, 0, 0): 744
+      (0, 4, 0, 0): 745
+      (0, 5, 0, 0): 746
+      (0, 6, 0, 0): 747
+      (0, 7, 0, 0): 748
+      (0, 8, 0, 0): 749
+      (0, 9, 0, 0): 750
+      (0, 10, 0, 0): 751
+      (0, 11, 0, 0): 752
+      (0, 12, 0, 0): 753
+      (0, 13, 0, 0): 754
+      (0, 14, 0, 0): 755
+      (0, 15, 0, 0): 756
+      (0, 16, 0, 0): 757
+      (0, 17, 0, 0): 758
+      (0, 18, 0, 0): 759
+      (0, 19, 0, 0): 760
+      (0, 20, 0, 0): 761
+      (0, 21, 0, 0): 762
+      (0, 0, 1, 0): 763
+      (0, 0, 2, 0): 764
+      (0, 0, 3, 0): 765
+      (0, 0, 4, 0): 766
+      (0, 0, 5, 0): 767
+      (0, 0, 6, 0): 768
+      (0, 0, 7, 0): 769
+      (0, 0, 8, 0): 770
+      (0, 0, 9, 0): 771
+      (0, 0, 10, 0): 772
+      (0, 0, 11, 0): 773
+      (0, 0, 12, 0): 774
+      (0, 0, 13, 0): 775
+      (0, 0, 14, 0): 776
+      (0, 0, 15, 0): 777
+      (0, 0, 0, 1): 778
+      (0, 0, 0, 2): 779
+      (0, 0, 0, 3): 780
+      (0, 0, 0, 4): 781
+      (0, 0, 0, 5): 782
+      (0, 0, 0, 6): 783
+      (0.5, 1): 784 # scale variations
+      (2, 1): 785
+      (0.5, 0.5): 786
+      (1, 0.5): 787
+      (2, 0.5): 788
+      (0.5, 2): 789
+      (1, 2): 790
+      (2, 2): 791
+      (-1, -1): 792 # DIS coeff functions fact variations
+      (1, 1): 793
+      # IHOU + MHOU missing prescription: see https://github.com/NNPDF/papers/blob/e2ac1832cf4a36dab83a696564eaa75a4e55f5d2/minutes/minutes-2023-08-18.txt#L148-L157
+      (1, 0.5 missing): 794  # As 787 but DIS and N3LO DY from 722
+      (1, 2 missing): 795 # As 790 but DIS and N3LO DY from 722
+      (1, 0.5 hadronic): 796  # As 787 but DIS from 722
+      (1, 2 hadronic): 797 # As 790 but DIS from 722
+
   - theoryid: 600
     variations:
       (0.5, 0.5): 601
diff --git a/validphys2/src/validphys/theorycovariance/construction.py b/validphys2/src/validphys/theorycovariance/construction.py
index 5feb6835df..db1a370d26 100644
--- a/validphys2/src/validphys/theorycovariance/construction.py
+++ b/validphys2/src/validphys/theorycovariance/construction.py
@@ -333,6 +333,45 @@ def covmat_9pt(name1, name2, deltas1, deltas2):
         ) + (1 / 8) * (np.outer((deltas1[2] + deltas1[3]), (deltas2[2] + deltas2[3])))
     return s
 
+def covmat_n3lo_singlet(name1, name2, deltas1, deltas2):
+    """Returns theory covariance sub-matrix for all the
+    singlet splitting function variations.
+    """
+    n3lo_vars_dict = {
+        "gg": 19,
+        "gq": 21,
+        "qg": 15,
+        "qq": 6,
+    }
+    s_singlet_ad = 0
+    cnt = 0
+    for n_var in n3lo_vars_dict.values():
+        s_singlet_ad += covmat_n3lo_ad(name1, name2, deltas1[cnt:cnt+n_var], deltas2[cnt:cnt+n_var])
+        cnt += n_var
+    return s_singlet_ad
+
+
+def covmat_n3lo_ad(name1, name2, deltas1, deltas2):
+    """Returns theory covariance sub-matrix for each of the
+    singlet splitting function variations.
+
+    Normalization is given by:
+        (n_pt - 1)
+    
+    where:
+        * n_pt = number of point presctiption
+    """
+    norm = len(deltas1)
+    if name1 == name2:
+        s = sum(np.outer(d, d) for d in deltas1)
+    else:
+        s = 0
+        for i, d1 in enumerate(deltas1):
+            for j, d2 in enumerate(deltas2):
+                if i == j:
+                    s += np.outer(d1, d2)
+    return 1 / norm * s
+
 
 @check_correct_theory_combination
 def covs_pt_prescrip(
@@ -384,6 +423,31 @@ def covs_pt_prescrip(
                     s = covmat_7pt(name1, name2, deltas1, deltas2)
             elif l == 9:
                 s = covmat_9pt(name1, name2, deltas1, deltas2)
+            # n3lo ad variation prescriprion
+            elif l == 62:
+                s = covmat_n3lo_singlet(name1, name2, deltas1, deltas2)
+            # n3lo ihou prescriprion
+            elif l == 64:
+                s_ad = covmat_n3lo_singlet(name1, name2, deltas1[:-2], deltas2[:-2])
+                s_cf = covmat_3pt(name1, name2, deltas1[-2:], deltas2[-2:])
+                s = s_ad + s_cf
+            # n3lo 3 pt MHOU see also 
+            # see https://github.com/NNPDF/papers/blob/e2ac1832cf4a36dab83a696564eaa75a4e55f5d2/minutes/minutes-2023-08-18.txt#L148-L157
+            elif l == 66:
+                s_ad = covmat_n3lo_singlet(name1, name2, deltas1[:-4], deltas2[:-4])
+                s_mhou = covmat_3pt(name1, name2, deltas1[-4:-2], deltas2[-4:-2])
+                s_cf = covmat_3pt(name1, name2, deltas1[-2:], deltas2[-2:])
+                s = s_ad + s_cf + s_mhou
+            # n3lo full covmat prescriprion
+            elif l == 70:
+                # spit deltas and compose thcovmat
+                # splittin functions variatons
+                s_ad = covmat_n3lo_singlet(name1, name2, deltas1[:-8], deltas2[:-8])
+                # scale variations
+                s_mhou = covmat_7pt(name1, name2, deltas1[-8:-2], deltas2[-8:-2])
+                # massive coefficient function variations
+                s_cf = covmat_3pt(name1, name2, deltas1[-2:], deltas2[-2:])
+                s = s_ad + s_cf + s_mhou
             start_locs = (start_proc[name1], start_proc[name2])
             covmats[start_locs] = s
     return covmats
diff --git a/validphys2/src/validphys/theorycovariance/output.py b/validphys2/src/validphys/theorycovariance/output.py
index dd959727e6..0487d81215 100644
--- a/validphys2/src/validphys/theorycovariance/output.py
+++ b/validphys2/src/validphys/theorycovariance/output.py
@@ -79,8 +79,8 @@ def matrix_plot_labels(df):
 def plot_covmat_heatmap(covmat, title):
     """Matrix plot of a covariance matrix."""
     df = covmat
-    df.sort_index(0, inplace=True)
-    df.sort_index(1, inplace=True)
+    df.sort_index(axis=0, inplace=True)
+    df.sort_index(axis=1, inplace=True)
     oldindex = df.index.tolist()
     newindex = sorted(oldindex, key=_get_key)
     # reindex index
@@ -181,8 +181,8 @@ def _get_key(element):
 def plot_corrmat_heatmap(corrmat, title):
     """Matrix plot of a correlation matrix"""
     df = corrmat
-    df.sort_index(0, inplace=True)
-    df.sort_index(1, inplace=True)
+    df.sort_index(axis=0, inplace=True)
+    df.sort_index(axis=1, inplace=True)
     oldindex = df.index.tolist()
     newindex = sorted(oldindex, key=_get_key)
     # reindex index
@@ -350,18 +350,18 @@ def plot_diag_cov_comparison(
     plot_index = theory_covmat_custom.index
     sqrtdiags_th = np.sqrt(np.diag(theory_covmat_custom)) / data
     sqrtdiags_th = pd.DataFrame(sqrtdiags_th.values, index=plot_index)
-    sqrtdiags_th.sort_index(0, inplace=True)
+    sqrtdiags_th.sort_index(axis=0, inplace=True)
     oldindex = sqrtdiags_th.index.tolist()
     newindex = sorted(oldindex, key=_get_key)
     sqrtdiags_th = sqrtdiags_th.reindex(newindex)
     sqrtdiags_exp = np.sqrt(np.diag(procs_covmat)) / data
     sqrtdiags_exp = pd.DataFrame(sqrtdiags_exp.values, index=plot_index)
-    sqrtdiags_exp.sort_index(0, inplace=True)
+    sqrtdiags_exp.sort_index(axis=0, inplace=True)
     sqrtdiags_exp = sqrtdiags_exp.reindex(newindex)
     df_total = theory_covmat_custom + procs_covmat
     sqrtdiags_tot = np.sqrt(np.diag(df_total)) / data
     sqrtdiags_tot = pd.DataFrame(sqrtdiags_tot.values, index=plot_index)
-    sqrtdiags_tot.sort_index(0, inplace=True)
+    sqrtdiags_tot.sort_index(axis=0, inplace=True)
     sqrtdiags_tot = sqrtdiags_tot.reindex(newindex)
     fig, ax = plotutils.subplots(figsize=(20, 10))
     ax.plot(sqrtdiags_exp.values, ".", label="Experiment", color="orange")
@@ -401,12 +401,12 @@ def plot_diag_cov_impact(
     inv_tot = (np.diag(la.inv(matrix_theory + matrix_experiment))) ** (-0.5) / procs_data_values
     plot_index = theory_covmat_custom.index
     df_inv_exp = pd.DataFrame(inv_exp, index=plot_index)
-    df_inv_exp.sort_index(0, inplace=True)
+    df_inv_exp.sort_index(axis=0, inplace=True)
     oldindex = df_inv_exp.index.tolist()
     newindex = sorted(oldindex, key=_get_key)
     df_inv_exp = df_inv_exp.reindex(newindex)
     df_inv_tot = pd.DataFrame(inv_tot, index=plot_index)
-    df_inv_tot.sort_index(0, inplace=True)
+    df_inv_tot.sort_index(axis=0, inplace=True)
     df_inv_tot = df_inv_tot.reindex(newindex)
     fig, ax = plotutils.subplots()
     ax.plot(df_inv_exp.values, ".", label="Experiment", color="orange")
diff --git a/validphys2/src/validphys/theorycovariance/tests.py b/validphys2/src/validphys/theorycovariance/tests.py
index fff466d676..88f5460927 100644
--- a/validphys2/src/validphys/theorycovariance/tests.py
+++ b/validphys2/src/validphys/theorycovariance/tests.py
@@ -683,7 +683,7 @@ def eigenvector_plot(evals_nonzero_basis, shx_vector):
         [processnames, dsnames, ids], names=("process", "dataset", "id")
     )
     f = pd.DataFrame(f.values, index=tripleindex)
-    f.sort_index(0, inplace=True)
+    f.sort_index(axis=0, inplace=True)
     oldindex = f.index.tolist()
     newindex = sorted(oldindex, key=_get_key)
     f = f.reindex(newindex)
@@ -737,12 +737,12 @@ def deltamiss_plot(theory_shift_test, allthx_vector, evals_nonzero_basis, shx_ve
         [processnames, dsnames, ids], names=("process", "dataset", "id")
     )
     f = pd.DataFrame(f.values, index=tripleindex)
-    f.sort_index(0, inplace=True)
+    f.sort_index(axis=0, inplace=True)
     oldindex = f.index.tolist()
     newindex = sorted(oldindex, key=_get_key)
     f = f.reindex(newindex)
     fmiss = pd.DataFrame(fmiss, index=tripleindex)
-    fmiss.sort_index(0, inplace=True)
+    fmiss.sort_index(axis=0, inplace=True)
     fmiss = fmiss.reindex(newindex)
     # Plotting
     fig, ax = plotutils.subplots(figsize=(20, 10))
@@ -786,15 +786,15 @@ def shift_diag_cov_comparison(allthx_vector, shx_vector, thx_covmat, thx_vector)
         [processnames, dsnames, ids], names=("process", "dataset", "id")
     )
     matrix = pd.DataFrame(matrix.values, index=tripleindex, columns=tripleindex)
-    matrix.sort_index(0, inplace=True)
-    matrix.sort_index(1, inplace=True)
+    matrix.sort_index(axis=0, inplace=True)
+    matrix.sort_index(axis=1, inplace=True)
     oldindex = matrix.index.tolist()
     newindex = sorted(oldindex, key=_get_key)
     matrix = matrix.reindex(newindex)
     matrix = (matrix.T.reindex(newindex)).T
     sqrtdiags = np.sqrt(np.diag(matrix))
     fnorm = pd.DataFrame(fnorm.values, index=tripleindex)
-    fnorm.sort_index(0, inplace=True)
+    fnorm.sort_index(axis=0, inplace=True)
     fnorm = fnorm.reindex(newindex)
     # Plotting
     fig, ax = plotutils.subplots(figsize=(20, 10))
diff --git a/validphys2/src/validphys/theorycovariance/theorycovarianceutils.py b/validphys2/src/validphys/theorycovariance/theorycovarianceutils.py
index 369d0b258c..0e04dccc4f 100644
--- a/validphys2/src/validphys/theorycovariance/theorycovarianceutils.py
+++ b/validphys2/src/validphys/theorycovariance/theorycovarianceutils.py
@@ -18,7 +18,7 @@ def check_correct_theory_combination_internal(
     """Checks that a valid theory combination corresponding to an existing
     prescription has been inputted"""
     l = len(theoryids)
-    check(l in {3, 5, 7, 9}, "Expecting exactly 3, 5, 7 or 9 theories, but got {l}.")
+    check(l in {3, 5, 7, 9, 62, 64, 66, 70}, f"Expecting exactly 3, 5, 7, 9, 62, 64, 66 or 70 theories, but got {l}.")
     opts = {"bar", "nobar"}
     xifs = [theoryid.get_description()["XIF"] for theoryid in theoryids]
     xirs = [theoryid.get_description()["XIR"] for theoryid in theoryids]
@@ -53,6 +53,62 @@ def check_correct_theory_combination_internal(
     elif l == 7:
         correct_xifs = [1.0, 2.0, 0.5, 1.0, 1.0, 2.0, 0.5]
         correct_xirs = [1.0, 1.0, 1.0, 2.0, 0.5, 2.0, 0.5]
+    elif l in [62, 64, 66, 70]:
+        # check Anomalous dimensions variations
+        n3lo_vars_dict = {
+            "gg": 19,
+            "gq": 21,
+            "qg": 15,
+            "qq": 6,
+        }
+        # TODO: for the moment fish the n3lo_ad_variation from the comments
+        n3lo_vars_list = []
+        id_max = None
+        if l == 70:
+            id_max = -8
+        elif l == 66:
+            id_max = -4
+        elif l == 64:
+            id_max = -2
+        for theoryid in theoryids[:id_max]:
+            n3lo_vars_list.append([int(val) for val in theoryid.get_description()["Comments"][28:-1].split(",")])
+        full_var_list = [[0,0,0,0]]
+        for entry, max_var in enumerate(n3lo_vars_dict.values()):
+            for idx in range(1, max_var + 1):
+                base_var = [0,0,0,0]
+                base_var[entry] = idx
+                full_var_list.append(base_var)
+        check(
+            n3lo_vars_list == full_var_list,
+            f"Theories do not include the full list of N3LO variation but {n3lo_vars_list}"
+        )
+        if l == 70:
+            # check Scale variations
+            varied_xifs = [xifs[0]]
+            varied_xirs = [xirs[0]]
+            varied_xifs.extend(xifs[-8:-2])
+            varied_xirs.extend(xirs[-8:-2])
+            correct_xifs = [1.0, 2.0, 0.5, 1.0, 1.0, 2.0, 0.5]
+            correct_xirs = [1.0, 1.0, 1.0, 2.0, 0.5, 2.0, 0.5]
+            check(
+                varied_xifs == correct_xifs and varied_xirs == correct_xirs,
+                "Choice of input theories does not correspond to a valid "
+                "prescription for theory covariance matrix calculation",
+            )
+        elif l == 66:
+            # check Scale variations
+            varied_xifs = [xifs[0]]
+            varied_xirs = [xirs[0]]
+            varied_xifs.extend(xifs[-4:-2])
+            varied_xirs.extend(xirs[-4:-2])
+            correct_xifs = [1.0, 1.0, 1.0]
+            correct_xirs = [1.0, 0.5, 2.0]
+            check(
+                varied_xifs == correct_xifs and varied_xirs == correct_xirs,
+                "Choice of input theories does not correspond to a valid "
+                "prescription for theory covariance matrix calculation",
+            )
+        return
     else:
         correct_xifs = [1.0, 2.0, 0.5, 1.0, 1.0, 2.0, 0.5, 2.0, 0.5]
         correct_xirs = [1.0, 1.0, 1.0, 2.0, 0.5, 2.0, 0.5, 0.5, 2.0]