Document QED solutions

[alecandido]

4x4 singlet exact

• even the exact solution with 4 flavors requires a numerical "upgrade"
  • for pure QCD the 2x2 vector is solved with projectors $\textbf{e}_{\pm}$
  • this does not scale to 4x4
  • but the solution in components its simple and extremely similar NxN-pegasus-truncated.pdf
    

Inconsistent $a a_s \cdot f_2$ solution

• we can implement the inconsiste $a a_s$ solution, neglecting the commutators in the first picture, if this corresponds to neglect $a a_s \cdot f_2$ terms in the second picture
  • this has to be proven
  • and then it has to be checked numerically not to be too drastic (coupling-wise is a correction of order NNLO QCD, since $a a_s \sim a_s^3$
  

$\varepsilon$-trick

• we can implement the $\varepsilon$ trick, but we won't
  • nevertheless it's good to document that we can
    
    the expression in the previous picture is easily proven to be equivalent to Eq. (11) of https://github.com/scarrazza/apfel/blob/master/doc/pdfs/TruncatedSolution.pdf
  • this can be implemented in the following replacing the a_half with a_half * epsilon after pure QED and LO QCD (using contracted gammas), and running with:
    • $\varepsilon = 0$ for LO exact
    • $\varepsilon = 1$ for $N^kLO$ exact
    • $\varepsilon = \pm 10^{-5}$ to implement the $\varepsilon$-trick (and then recombine before the Mellin inversion)

eko/src/eko/kernels/singlet.py

Lines 331 to 351 in 13197aa

	a_steps = utils.geomspace(a0, a1, 1 + ev_op_iterations)
	beta_vec = [beta.beta_qcd((2, 0), nf), beta.beta_qcd((3, 0), nf)]
	if order[0] >= 3:
	beta_vec.append(beta.beta_qcd((4, 0), nf))
	if order[0] >= 4:
	beta_vec.append(beta.beta_qcd((5, 0), nf))
	e = np.identity(2, np.complex_)
	al = a_steps[0]
	for ah in a_steps[1:]:
	a_half = (ah + al) / 2.0
	delta_a = ah - al
	gamma_summed = np.zeros((2, 2), dtype=np.complex_)
	beta_summed = 0
	for i in range(order[0]):
	gamma_summed += gamma_singlet[i] * a_half**i
	beta_summed += beta_vec[i] * a_half ** (i + 1)
	ln = gamma_summed / beta_summed * delta_a
	ek = np.ascontiguousarray(ad.exp_singlet(ln)[0])
	e = ek @ e
	al = ah
	return e