From 16d69f24a33e27ad57e451d29d1d8ecf3d992c0c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= Date: Tue, 9 Jan 2024 19:36:44 +0100 Subject: [PATCH] Implement varying nl in the coupling running --- src/eko/constants.py | 22 ++++++++++++++++ src/eko/couplings.py | 36 ++++++++++++++++++-------- src/eko/evolution_operator/__init__.py | 28 ++------------------ 3 files changed, 49 insertions(+), 37 deletions(-) diff --git a/src/eko/constants.py b/src/eko/constants.py index fba786f46..1e00d0c62 100644 --- a/src/eko/constants.py +++ b/src/eko/constants.py @@ -25,6 +25,9 @@ Defaults to :math:`\frac{N_C^2-1}{2N_C} = 4/3`. """ +MTAU = 1.777 +"""Mass of the tau.""" + eu2 = 4.0 / 9 """Up quarks charge squared.""" @@ -114,3 +117,22 @@ def charge_combinations(nf): vde2m = nd / nf * (eu2 - ed2) e2delta = vde2m - vue2m + e2avg return e2avg, vue2m, vde2m, e2delta + + +@nb.njit(cache=True) +def lepton_number(q2): + """Compute the number of up flavors. + + at the moment we never go below 1 GeV so we don't need muon and electron + + Parameters + ---------- + q : float + energy + + Returns + ------- + nl : int + + """ + return 3 if q2 > MTAU**2 else 2 diff --git a/src/eko/couplings.py b/src/eko/couplings.py index b0d39a070..bd6e658db 100644 --- a/src/eko/couplings.py +++ b/src/eko/couplings.py @@ -277,7 +277,7 @@ def expanded_qed(ref, order, beta0, b_vec, lmu): @nb.njit(cache=True) def couplings_expanded_alphaem_running( - order, couplings_ref, nf, scale_from, scale_to, decoupled_running + order, couplings_ref, nf, nl, scale_from, scale_to, decoupled_running ): r"""Compute coupled expanded expression of the couplings for running alphaem. @@ -308,8 +308,8 @@ def couplings_expanded_alphaem_running( beta0_qcd = beta_qcd((2, 0), nf) b_vec_qcd = [b_qcd((i + 2, 0), nf) for i in range(order[0])] res_as = expanded_qcd(couplings_ref[0], order[0], beta0_qcd, b_vec_qcd, lmu) - beta0_qed = beta_qed((0, 2), nf) - b_vec_qed = [b_qed((0, i + 2), nf) for i in range(order[1])] + beta0_qed = beta_qed((0, 2), nf, nl) + b_vec_qed = [b_qed((0, i + 2), nf, nl) for i in range(order[1])] res_aem = expanded_qed(couplings_ref[1], order[1], beta0_qed, b_vec_qed, lmu) # if order[0] >= 1 and order[1] >= 1: # order[0] is always >=1 @@ -521,7 +521,7 @@ def rge(_t, a, b_vec): ) return res.y[0][-1] - def compute_exact_alphaem_running(self, a_ref, nf, scale_from, scale_to): + def compute_exact_alphaem_running(self, a_ref, nf, nl, scale_from, scale_to): """Compute couplings via |RGE| with running alphaem. Parameters @@ -576,12 +576,12 @@ def compute_exact_alphaem_running(self, a_ref, nf, scale_from, scale_to): # while aem is constant return np.array([rge_qcd, a_ref[1]]) if self.order[1] >= 1: - beta_qed_vec = [beta_qed((0, 2), nf)] + beta_qed_vec = [beta_qed((0, 2), nf, nl)] if not self.decoupled_running: beta_qcd_mix = beta_qcd((2, 1), nf) - beta_qed_mix = beta_qed((1, 2), nf) # order[0] is always at least 1 + beta_qed_mix = beta_qed((1, 2), nf, nl) # order[0] is always at least 1 if self.order[1] >= 2: - beta_qed_vec.append(beta_qed((0, 3), nf)) + beta_qed_vec.append(beta_qed((0, 3), nf, nl)) # integration kernel def rge(_t, a, beta_qcd_vec, beta_qcd_mix, beta_qed_vec, beta_qed_mix): @@ -656,7 +656,7 @@ def compute_exact_fixed_alphaem(self, a_ref, nf, scale_from, scale_to): ) return np.array([rge_qcd, a_ref[1]]) - def compute(self, a_ref, nf, scale_from, scale_to): + def compute(self, a_ref, nf, nl, scale_from, scale_to): """Compute actual couplings. This is a wrapper in order to pass the computation to the corresponding method @@ -678,7 +678,7 @@ def compute(self, a_ref, nf, scale_from, scale_to): numpy.ndarray couplings at target scale :math:`a(Q^2)` """ - key = (float(a_ref[0]), float(a_ref[1]), nf, scale_from, float(scale_to)) + key = (float(a_ref[0]), float(a_ref[1]), nf, nl, scale_from, float(scale_to)) try: return self.cache[key].copy() except KeyError: @@ -686,7 +686,7 @@ def compute(self, a_ref, nf, scale_from, scale_to): if self.method == "exact": if self.alphaem_running: a_new = self.compute_exact_alphaem_running( - a_ref.astype(float), nf, scale_from, scale_to + a_ref.astype(float), nf, nl, scale_from, scale_to ) else: a_new = self.compute_exact_fixed_alphaem( @@ -698,6 +698,7 @@ def compute(self, a_ref, nf, scale_from, scale_to): self.order, a_ref.astype(float), nf, + nl, scale_from, float(scale_to), self.decoupled_running, @@ -737,7 +738,20 @@ def a( for k, seg in enumerate(path): # skip a very short segment, but keep the matching if not np.isclose(seg.origin, seg.target): - new_a = self.compute(final_a, seg.nf, seg.origin, seg.target) + nli = constants.lepton_number(seg.origin) + nlf = constants.lepton_number(seg.target) + if nli != nlf: + # it means that MTAU is between origin and target: + # first we evolve from origin to MTAU with nli leptons + a_tmp = self.compute( + final_a, seg.nf, nli, seg.origin, constants.MTAU**2 + ) + # then from MTAU to target with nlf leptons + new_a = self.compute( + a_tmp, seg.nf, nlf, constants.MTAU**2, seg.target + ) + else: + new_a = self.compute(final_a, seg.nf, nli, seg.origin, seg.target) else: new_a = final_a # apply matching conditions: see hep-ph/9706430 diff --git a/src/eko/evolution_operator/__init__.py b/src/eko/evolution_operator/__init__.py index 29e67c19c..8092d31a0 100644 --- a/src/eko/evolution_operator/__init__.py +++ b/src/eko/evolution_operator/__init__.py @@ -704,38 +704,14 @@ def compute_aem_list(self): as_list = np.array([self.a_s[0], self.a_s[1]]) a_half = np.zeros((ev_op_iterations, 2)) else: - as0 = self.a_s[0] - as1 = self.a_s[1] - aem0 = self.a_em[0] - aem1 = self.a_em[1] - q2ref = self.managers["couplings"].mu2_ref - delta_from = abs(self.q2_from - q2ref) - delta_to = abs(self.q2_to - q2ref) - # I compute the values in aem_list starting from the mu2 - # that is closer to mu_ref. - if delta_from > delta_to: - a_start = np.array([as1, aem1]) - mu2_start = self.q2_to - else: - a_start = np.array([as0, aem0]) - mu2_start = self.q2_from couplings = self.managers["couplings"] mu2_steps = utils.geomspace(self.q2_from, self.q2_to, 1 + ev_op_iterations) mu2_l = mu2_steps[0] - as_list = np.array( - [ - couplings.compute( - a_ref=a_start, nf=self.nf, scale_from=mu2_start, scale_to=mu2 - )[0] - for mu2 in mu2_steps - ] - ) + as_list = np.array([couplings.a_s(scale_to=mu2)[0] for mu2 in mu2_steps]) a_half = np.zeros((ev_op_iterations, 2)) for step, mu2_h in enumerate(mu2_steps[1:]): mu2_half = (mu2_h + mu2_l) / 2.0 - a_s, aem = couplings.compute( - a_ref=a_start, nf=self.nf, scale_from=mu2_start, scale_to=mu2_half - ) + a_s, aem = couplings.a(scale_to=mu2_half) a_half[step] = [a_s, aem] mu2_l = mu2_h return as_list, a_half