$\mu$-$e$ two-loop matrix element

This work was done in collaboration with Alessandro Broggio, William J. Torres Bobadilla, Andrea Ferroglia, Manoj Kumar Mandal, Pierpaolo Mastrolia, Jonathan Ronca, Max Zoller

from pymule import *

Populating the interactive namespace from numpy and matplotlib

We can plot the two-loop matrix element with $m=0$ \begin{align} \mathcal{M}_{n}^{(2)} = 2{\rm Re}\Big[\mathcal{A}_n^{(2)} \times \mathcal{A}_n^{(0)*}\Big]+ \Big|\mathcal{A}_n^{(1)}\Big|^2 \end{align} We set $\alpha$ and $M$ to their physical values and

\begin{align} s &= m^2 + M^2 + 2 m \times 160\,{\rm GeV} = 0.1746\,{\rm GeV}^2\\ t &\in [-0.153069, -0.00102148]\,{\rm GeV}^2 \end{align}

We can now plot the finite part of the FDH matrix element

SUBROUTINE TESTMUEPLOT
use mue_mat_el
integer i
integer, parameter :: n = 100
real(kind=prec) :: alpha_real = 1/137.035999084_prec
real(kind=prec) :: tt, tmin, tmax, E1
real(kind=prec) :: twoloop(-4:0, 5)

call blockstart("mu-e plot")
call initflavour('muone160+')
musq = mm**2

E1 = (scms + Me**2 - Mm**2)/2/sqrt(scms)
tmin = -4*(E1**2 - Me**2)
tmax = 2*Me**2 - 2*Me*1000.

do i = 0, n-1
tt = tmin + (tmax - tmin) / n * (0.5 + i)
call emem(scms, tt, Mm**2, twoloop=twoloop)
write(7,*) tt, alpha_real**4*twoloop(0, :)
enddo

END SUBROUTINE

dat = np.loadtxt('two-loop-plot.txt.bz2')

for i in range(1, len(dat[0])):
    plot(dat[:-1,0], 1e5*dat[:-1,i])
ylim(-2.2, 4.255)
legend([
    f"$q^{7-i}Q^{1+i}$"
    for i in range(1,6)
], loc=3, ncol=3)
ylabel(r"$10^5\mathcal{M}_{n}^{(2)}\Big|_\text{finite}$")
xlabel(r"$t\,/\,{\rm MeV}^2$")
gcf().savefig("plots/two-loop.pdf")