diff --git a/doc/source/theory/pQCD.rst b/doc/source/theory/pQCD.rst
index 2223b327c..27c5cb543 100644
--- a/doc/source/theory/pQCD.rst
+++ b/doc/source/theory/pQCD.rst
@@ -86,8 +86,9 @@ At |NLO|, the singlet entry of the quark-quark anomalous dimension can be decomp
     \gamma^{(1)}_{qq} =\gamma^{(1)}_{ps} + \gamma^{(1)}_{ns,+}
 
 The non-singlet sector in the polarized case swaps the plus and minus non-singlet relative to the unpolarized case.
-This is because the polarized non-singlet splitting functions are defined as the difference between the probability of the polarized parton splitting into daughter partons of the same flavour
-and daughters splitting into a different flavours and opposite helicity. The first moments of the anomalous dimensions are:
+This is because the polarized non-singlet splitting functions are defined as the difference between the probability of the polarized parton splitting 
+into daughter partons of the same flavour and same helicity and daughters of a different flavours and opposite helicity. 
+The first moments of the anomalous dimensions are:
 
 .. math ::
     \gamma^{(1)}_{ns,+} (N=1) &= 0 \\
@@ -102,6 +103,12 @@ At |NNLO| the non-singlet is further decomposed into the helicity difference qua
 where :math:`\gamma^{(2)}_{ns,-}` is the minus flavour asymmetry non-singlet and :math:`\gamma^{(2)}_{ns,s}` the sea-like polarized non-singlet.
 The singlet entry :math:`\gamma^{(2)}_{qq}` is defined as above in the |NLO| case.
 
+Finally the violation of the axial current conservation :math:`\bar{\psi} \gamma_\mu \gamma_5 \bar{\psi}` only through 
+loop corrections impose the following relations to the singlet splittings at all orders :cite:`Moch:2014sna` :
+
+.. math ::
+     \gamma^{(n)}_{gg} & = - \beta_n \\
+     \gamma^{(n)}_{ps} & = - 2 n_f \gamma^{(n-1)}_{gq}
 
 Unified Splitting Functions
 ---------------------------
diff --git a/tests/ekore/anomalous_dimensions/polarized/space_like/test_ad_as1.py b/tests/ekore/anomalous_dimensions/polarized/space_like/test_ad_as1.py
index 4b7b24d38..ad0a3363d 100755
--- a/tests/ekore/anomalous_dimensions/polarized/space_like/test_ad_as1.py
+++ b/tests/ekore/anomalous_dimensions/polarized/space_like/test_ad_as1.py
@@ -2,7 +2,7 @@
 import numpy as np
 
 import ekore.anomalous_dimensions.polarized.space_like.as1 as as1
-from eko import constants
+from eko import beta, constants
 from ekore import harmonics
 
 NF = 5
@@ -30,3 +30,10 @@ def test_gluon_momentum():
 def test_qg_helicity_conservation():
     N = complex(1.0, 0.0)
     np.testing.assert_almost_equal(as1.gamma_qg(N, NF), 0)
+
+
+def test_axial_anomaly():
+    # violation of the axial current conservation happens only through loops
+    N = complex(1.0, 0.0)
+    cache = harmonics.cache.reset()
+    np.testing.assert_allclose(as1.gamma_gg(N, cache, NF), -beta.beta_qcd_as2(NF))
diff --git a/tests/ekore/anomalous_dimensions/polarized/space_like/test_ad_as2.py b/tests/ekore/anomalous_dimensions/polarized/space_like/test_ad_as2.py
index b77fd7a83..06ed35541 100755
--- a/tests/ekore/anomalous_dimensions/polarized/space_like/test_ad_as2.py
+++ b/tests/ekore/anomalous_dimensions/polarized/space_like/test_ad_as2.py
@@ -1,7 +1,9 @@
 # Test NLO Polarized splitting functions
 import numpy as np
 
+import ekore.anomalous_dimensions.polarized.space_like.as1 as as1
 import ekore.anomalous_dimensions.polarized.space_like.as2 as as2
+from eko import beta
 from eko.constants import CA, CF, TR
 from ekore import harmonics
 
@@ -32,7 +34,10 @@ def test_qg_momentum():
     cache = harmonics.cache.reset()
     np.testing.assert_allclose(
         -as2.gamma_qg(N, nf, cache),
-        4 * nf * (0.574074 * CF - 2 * CA * (-7 / 18 + 1 / 6 * (5 - np.pi**2 / 3))) * TR,
+        4
+        * nf
+        * (0.574074 * CF - 2 * CA * (-7 / 18 + 1 / 6 * (5 - np.pi**2 / 3)))
+        * TR,
     )
 
 
@@ -55,5 +60,16 @@ def test_gg_momentum():
     cache = harmonics.cache.reset()
     np.testing.assert_almost_equal(
         -as2.gamma_gg(N, nf, cache),
-        4 * (-1.7537256813471833 * CA**2 + ((29 * CA) / 27 - (28 * CF) / 27) * nf * TR),
+        4
+        * (-1.7537256813471833 * CA**2 + ((29 * CA) / 27 - (28 * CF) / 27) * nf * TR),
+    )
+
+
+def test_axial_anomaly():
+    # violation of the axial current conservation happens only through loops
+    N = complex(1.0, 0.0)
+    cache = harmonics.cache.reset()
+    np.testing.assert_allclose(
+        as2.gamma_gg(N, nf, cache), -beta.beta_qcd_as3(nf), rtol=9e-7
     )
+    np.testing.assert_allclose(as2.gamma_ps(N, nf), -2 * nf * as1.gamma_gq(N))
diff --git a/tests/ekore/anomalous_dimensions/polarized/space_like/test_ad_as3.py b/tests/ekore/anomalous_dimensions/polarized/space_like/test_ad_as3.py
index 82bd48675..fc75fedbe 100755
--- a/tests/ekore/anomalous_dimensions/polarized/space_like/test_ad_as3.py
+++ b/tests/ekore/anomalous_dimensions/polarized/space_like/test_ad_as3.py
@@ -1,8 +1,10 @@
 # Test NNLO Polarized splitting functions
 import numpy as np
 
+import ekore.anomalous_dimensions.polarized.space_like.as2 as as2
 import ekore.anomalous_dimensions.polarized.space_like.as3 as as3
 import ekore.anomalous_dimensions.unpolarized.space_like.as3 as as3_unpol
+from eko import beta
 from eko.constants import zeta2, zeta3
 from ekore import harmonics
 
@@ -54,3 +56,15 @@ def test_ns():
     np.testing.assert_allclose(
         as3_unpol.gamma_nsp(N, nf, cache), as3.gamma_nsm(N, nf, cache)
     )
+
+
+def test_axial_anomaly():
+    # violation of the axial current conservation happens only through loops
+    N = complex(1.0, 0.0)
+    cache = harmonics.cache.reset()
+    np.testing.assert_allclose(
+        as3.gamma_gg(N, nf, cache), -beta.beta_qcd_as4(nf), rtol=2e-5
+    )
+    np.testing.assert_allclose(
+        as3.gamma_ps(N, nf, cache), -2 * nf * as2.gamma_gq(N, nf, cache), rtol=3e-6
+    )