Limit Q2 range for expanded as evolution.

src/eko/couplings.py

@@ -57,6 +57,14 @@ def couplings_mod_ev(mod_ev: EvolutionMethod) -> CouplingEvolutionMethod:
     raise ValueError(f"Unknown evolution mode {mod_ev}")
 
 
+@nb.njit(cache=True)
+def check_expanded_low_q2_limit(values):
+    """Check if the expanded evolution is feasible."""
+    for v in values:
+        if v < 0.0:
+            raise ValueError("Too low Q2, can't evolve alphas with truncated solution.")
+
+
 @nb.njit(cache=True)
 def exact_lo(ref, beta0, lmu):
     r"""Compute expanded solution at |LO|.
@@ -103,6 +111,7 @@ def expanded_nlo(ref, beta0, b1, lmu):
         coupling at target scale :math:`a(\mu_R^2)`
     """
     den = 1.0 + beta0 * ref * lmu
+    check_expanded_low_q2_limit((den,))
     a_LO = exact_lo(ref, beta0, lmu)
     as_NLO = a_LO * (1 - b1 * a_LO * np.log(den))
     return as_NLO
@@ -318,16 +327,11 @@ def couplings_expanded_alphaem_running(
     # order[0] is always >=1
     if not decoupled_running:
         if order[1] >= 1:
-            res_as += (
-                -couplings_ref[0] ** 2
-                * b_qcd((2, 1), nf)
-                * np.log(1 + beta0_qcd * couplings_ref[1] * lmu)
-            )
-            res_aem += (
-                -couplings_ref[1] ** 2
-                * b_qed((1, 2), nf, nl)
-                * np.log(1 + beta0_qed * couplings_ref[0] * lmu)
-            )
+            den_s = 1 + beta0_qcd * couplings_ref[1] * lmu
+            den_em = 1 + beta0_qed * couplings_ref[0] * lmu
+            check_expanded_low_q2_limit((den_s, den_em))
+            res_as += -couplings_ref[0] ** 2 * b_qcd((2, 1), nf) * np.log(den_s)
+            res_aem += -couplings_ref[1] ** 2 * b_qed((1, 2), nf, nl) * np.log(den_em)
     return np.array([res_as, res_aem])
 
 

tests/eko/test_couplings.py

@@ -13,6 +13,7 @@
 from eko.io.types import EvolutionMethod
 from eko.quantities.couplings import CouplingEvolutionMethod, CouplingsInfo
 from eko.quantities.heavy_quarks import MatchingScales, QuarkMassScheme
+from ekobox.cards import example
 
 masses = [m**2 for m in (2.0, 4.5, 175.0)]
 
@@ -327,3 +328,21 @@ def benchmark_expanded_n3lo(self):
         np.testing.assert_allclose(
             mathematica_val, as_N3LO.a(Q2)[0] - as_NNLO.a(Q2)[0], rtol=5e-7
         )
+
+
+def test_expanded_low_bound():
+    # test low Q2 expanded solution limit
+    tcard = example.theory()
+    evmod = CouplingEvolutionMethod.EXPANDED
+    for order in [(2, 0), (3, 0), (4, 0), (2, 1), (3, 1), (4, 1)]:
+        sc = Couplings(
+            tcard.couplings,
+            order,
+            evmod,
+            masses=[(x.value) ** 2 for x in tcard.heavy.masses],
+            hqm_scheme=tcard.heavy.masses_scheme,
+            thresholds_ratios=[1, 1, 1],
+        )
+        assert not np.isnan(sc.a_s(0.04))
+        with pytest.raises(ValueError):
+            sc.a_s(0.03)