Implement QED exact solution

benchmarks/CT18_bench.py

@@ -53,6 +53,29 @@ def benchmark_nnlo(self, Q0=1.295, Q2grid=(1e4,)):
         ]
         self.run([theory_card], [operator_card], ["CT18NNLO"])
 
+    def benchmark_nnlo_qed(self, Q0=1.295, Q2grid=(1e4,)):
+        theory_card = base_theory.copy()
+        theory_card.update(
+            {
+                "alphas": 0.118000,
+                "alphaqed": 0.007496,
+                "PTO": 2,
+                "QED": 1,
+                "Q0": Q0,
+                "MaxNfPdf": 5,
+                "MaxNfAs": 5,
+                "alphaem_running": True,
+            }
+        )
+        operator_card = {"Q2grid": list(Q2grid)}
+        self.skip_pdfs = lambda _theory: [
+            -6,
+            6,
+            "Tu8",
+            "Vu8",
+        ]
+        self.run([theory_card], [operator_card], ["CT18qed"])
+
     def benchmark_znnlo(self, Q0=1.3, Q2grid=(1e4,)):
         theory_card = base_theory.copy()
         theory_card.update(
@@ -81,4 +104,4 @@ def benchmark_znnlo(self, Q0=1.3, Q2grid=(1e4,)):
 
 if __name__ == "__main__":
     b = BenchmarkCT18()
-    b.benchmark_nnlo()
+    b.benchmark_nnlo_qed()

benchmarks/NNPDF_bench.py

@@ -62,6 +62,35 @@ def benchmark_nlo(self, Q0=1.65, Q2grid=(100,)):
         self.run([theory_card], [operator_card], ["NNPDF31_nlo_as_0118"])
 
 
+class BenchmarkNNPDF31_luxqed(BenchmarkNNPDF):
+    """Benchmark NNPDF3.1_luxqed"""
+
+    def benchmark_nnlo(self, Q0=1.65, Q2grid=(100,)):
+        theory_card = {
+            **base_theory,
+            "PTO": 2,
+            "QED": 2,
+            "Q0": Q0,
+        }
+        theory_card.update(
+            {"ModEv": "iterate-exact", "FNS": "VFNS", "alphaem_running": True}
+        )
+
+        self.skip_pdfs = lambda _theory: [
+            -6,
+            6,
+            "Tu8",
+            "Vu8",
+        ]
+
+        operator_card = {
+            **base_operator,
+            "Q2grid": list(Q2grid),
+            "ev_op_iterations": 10,
+        }
+        self.run([theory_card], [operator_card], ["NNPDF31_nnlo_as_0118_luxqed"])
+
+
 class BenchmarkNNPDF40(BenchmarkNNPDF):
     """Benchmark NNPDF4.0"""
 
@@ -98,6 +127,8 @@ def benchmark_nnlo(self, Q0=1.65, Q2grid=(100,)):
     # nn31.benchmark_nlo(Q0=np.sqrt(low2), Q2grid=[10])
     # # test backward
     # #nn31.benchmark_nlo(Q0=np.sqrt(high2), Q2grid=[low2])
-    nn40 = BenchmarkNNPDF40()
+    nn31qed = BenchmarkNNPDF31_luxqed()
+    nn31qed.benchmark_nnlo()
+    # nn40 = BenchmarkNNPDF40()
     # nn40.benchmark_nnlo(Q2grid=[100])
-    nn40.benchmark_nnlo(Q0=np.sqrt(high2), Q2grid=[low2])
+    # nn40.benchmark_nnlo(Q0=np.sqrt(high2), Q2grid=[low2])

benchmarks/apfel_bench.py

@@ -187,12 +187,137 @@ def benchmark_sv(self, pto, svmode):
         self.run(ts, operators.build(operators.apfel_config), ["ToyLH"])
 
 
+class BenchmarkFFNS_qed(ApfelBenchmark):
+    """Benckmark FFNS"""
+
+    ffns_theory = {
+        "Qref": 91.1870,
+        "mc": 1.3,
+        "mb": 4.75,
+        "mt": 172.0,
+        "FNS": "FFNS",
+        "ModEv": [
+            "EXA",
+            # "EXP",
+            # "TRN",
+        ],
+        "NfFF": 5,
+        "kcThr": 0.0,
+        "kbThr": 0.0,
+        "ktThr": np.inf,
+        "Q0": 5.0,
+        "alphas": 0.118000,
+        "alphaqed": 0.007496,
+        "alphaem_running": True,
+    }
+    ffns_theory = tolist(ffns_theory)
+
+    def benchmark_plain(self, pto, qed):
+        """Plain configuration"""
+
+        th = self.ffns_theory.copy()
+        th.update({"PTO": [pto], "QED": [qed]})
+        self.skip_pdfs = lambda _theory: [
+            -6,
+            6,
+            "Tu8",
+            "Vu8",
+        ]
+        self.run(
+            cartesian_product(th),
+            operators.build(operators.apfel_config),
+            ["NNPDF31_nnlo_as_0118_luxqed"],
+        )
+
+    def benchmark_sv(self, pto, qed, svmode):
+        """Scale Variation"""
+
+        ts = []
+        th = self.ffns_theory.copy()
+        th.update(
+            {
+                "PTO": [pto],
+                "QED": [qed],
+                "XIR": [np.sqrt(0.5)],
+                "fact_to_ren_scale_ratio": [np.sqrt(2.0)],
+                "ModSV": [svmode],
+                "EScaleVar": [0],
+            }
+        )
+        ts.extend(cartesian_product(th))
+        th = self.ffns_theory.copy()
+        th.update(
+            {
+                "PTO": [pto],
+                "QED": [qed],
+                "XIR": [np.sqrt(2.0)],
+                "fact_to_ren_scale_ratio": [np.sqrt(0.5)],
+                "ModSV": [svmode],
+                "EScaleVar": [0],
+            }
+        )
+        ts.extend(cartesian_product(th))
+        self.skip_pdfs = lambda _theory: [
+            -6,
+            6,
+            "Tu8",
+            "Vu8",
+        ]
+        self.run(ts, operators.build(operators.apfel_config), ["ToyLH"])
+
+
+class BenchmarkVFNS_qed(ApfelBenchmark):
+    """Benckmark FFNS"""
+
+    vfns_theory = {
+        "Qref": 91.1870,
+        "mc": 1.3,
+        "mb": 4.75,
+        "mt": 172.0,
+        "FNS": "VFNS",
+        "ModEv": [
+            "EXA",
+            # "EXP",
+            # "TRN",
+        ],
+        "kcThr": 1.0,
+        "kbThr": 1.0,
+        "ktThr": 1.0,
+        "Q0": 1.25,
+        "alphas": 0.118000,
+        "alphaqed": 0.007496,
+        "alphaem_running": True,
+    }
+    vfns_theory = tolist(vfns_theory)
+
+    def benchmark_plain(self, pto, qed):
+        """Plain configuration"""
+
+        th = self.vfns_theory.copy()
+        th.update({"PTO": [pto], "QED": [qed]})
+        self.skip_pdfs = lambda _theory: [
+            -6,
+            6,
+            "Tu8",
+            "Vu8",
+        ]
+        self.run(
+            cartesian_product(th),
+            operators.build(operators.apfel_config),
+            ["NNPDF31_nnlo_as_0118_luxqed"],
+        )
+
+
 if __name__ == "__main__":
 
-    obj = BenchmarkVFNS()
+    # obj = BenchmarkVFNS()
     # obj = BenchmarkFFNS()
 
     # obj.benchmark_plain(2)
-    obj.benchmark_sv(2, "exponentiated")
+    # obj.benchmark_sv(2, "exponentiated")
     # obj.benchmark_kthr(2)
     # obj.benchmark_msbar(2)
+
+    # obj = BenchmarkFFNS_qed()
+    obj = BenchmarkVFNS_qed()
+    obj.benchmark_plain(2, 2)

doc/source/theory/DGLAP.rst

@@ -172,7 +172,7 @@ Here the strategies are:
 - for ``method in ['iterate-exact', 'iterate-expanded']`` we use a discretized path-ordering :cite:`Bonvini:2012sh`:
 
 .. math::
-    \ESk{1}{a_s}{a_s^0} = \prod\limits_{k=n}^{0} \ESk{1}{a_s^{k+1}}{a_s^{k}}\quad \text{with} a_s^{n+1} = a_s
+    \ESk{1}{a_s}{a_s^0} = \prod\limits_{k=n}^{0} \ESk{1}{a_s^{k+1}}{a_s^{k}}\quad \text{with}\quad a_s^{n+1} = a_s
 
 where the order of the product is such that later |EKO| are to the left and
 

extras/qed_matching/Makefile

@@ -0,0 +1,9 @@
+all: matching_and_rotation.pdf
+
+%.pdf: %.tex
+	pdflatex $<
+
+clean:
+	rm -f *.aux
+	rm -f *.log
+	rm -f *.pdf

extras/qed_matching/matching_and_rotation.tex

@@ -0,0 +1,141 @@
+\documentclass[a4paper,oneside]{article}
+\usepackage[utf8]{inputenc}
+\usepackage{xcolor}
+\usepackage{amsmath}
+\usepackage{amssymb}
+\usepackage{amsfonts}
+
+%\usepackage[a4paper,top=3cm,bottom=3cm,left=3cm,right=3cm]{geometry}
+\usepackage[a4paper,top=1.5cm,bottom=1.5cm,left=1.5cm,right=1.5cm]{geometry}
+
+
+\title{Matching and Basis Rotation for the Intrinsic Unified Evolution Basis}
+\author{Niccolò Laurenti}
+
+\date{}
+
+\begin{document}
+
+\maketitle
+
+In this document we will explain how the matching and the basis rotation are performed in the Intrinsic Unified Evolution Basis.
+
+\section{Matching}
+Fist of all, we have to perform the matching between the basis vectors in the two different flavor schemes, i.e.\ the $n_f$ and the $n_f+1$ flavor schemes. The matching of the gluon, light quark and heavy quark PDFs are the followings:
+\begin{align*}
+g^{(n_f+1)}&=A^S_{gg}g^{(n_f)}+A^S_{gq}\Sigma^{(n_f)}_{(n_f)}+A^S_{gH}h^{(n_f)} \\
+l^{(n_f+1)} &= A^{ns}_{qq} l^{(n_f)} + {1\over2n_f} A^S_{qg}g^{(n_f)} \quad \text{and the same for $\bar{l}$}\\
+h^{(n_f+1)} &= {1\over2} A^{S}_{Hg} g^{(n_f)}_{(n_f)}+ {1\over2} A^{ps}_{Hq} \Sigma^{(n_f)}_{(n_f)} + A_{HH}h^{(n_f)} \quad \text{and the same for $\bar{h}$}
+\end{align*}
+From the second relation we get that
+\begin{align*}
+l^{(n_f+1)}_+ &= A^{ns}_{qq} l^{(n_f)}_+ + {1\over n_f} A^S_{qg}g^{(n_f)} \\
+l^{(n_f+1)}_- &= A^{ns}_{qq} l^{(n_f)}_- \\
+\Sigma^{(n_f+1)}_{(n_f)}&=A^S_{qg}g^{(n_f)}+A^{ns}_{qq}\Sigma^{(n_f)}_{(n_f)} \\
+V^{(n_f+1)}_{(n_f)}&=A^{ns}_{qq}V^{(n_f)}_{(n_f)}
+\end{align*}
+while from the third relation we get that
+\begin{align*}
+h^{(n_f+1)}_+ &= A^{S}_{Hg} g^{(n_f)}_{(n_f)}+A^{ps}_{Hq} \Sigma^{(n_f)}_{(n_f)} + A_{HH}h^{(n_f)}_+ \\
+h^{(n_f+1)}_- &= A_{HH}h^{(n_f)}_-
+\end{align*}
+ The matching of the components $\Sigma_{\Delta(n_f)}$, $V_{\Delta(n_f)}$, $T_i$, $V_i$ are diagonal. For the $V$s this is trivial since they are composed by $l_-$. For $\Sigma_{\Delta(n_f)}$ instead: being it defined as
+ \begin{align*}
+ \Sigma^{(n_f+1)}_{\Delta(n_f)} & = \frac{n_d(n_f)}{n_u(n_f)} \sum_{i=1}^{n_u} u_{+\,i}^{(n_f+1)} - \sum_{i=1}^{n_d} d_{+\,i}^{(n_f+1)}
+ \end{align*}
+ the gluon contribution cancels, giving the relation
+ \begin{equation*}
+ \Sigma^{(n_f+1)}_{\Delta(n_f)}=A^{ns}_{qq}\Sigma^{(n_f)}_{\Delta(n_f)}
+\end{equation*}
+The same holds for the $T_i$ components.
+
+Observe that this holds up to NNLO, since that at N$^3$LO we have to consider also the pure singlet components of the light quark matching (I think that we have to add $\frac{1}{2n_f}A^{ps}_{qq}\Sigma^{nf}_{(nf)}$ to the matching of $l^{(n_f+1)}$, but I'm not 100\% sure. In this way the matching of $l_-$ remains diagonal, as it should be, and the same hold for $\Sigma_{\Delta(n_f)}$ and $T_i$).
+
+Up to second order the perturbative expansion of the matching terms is given by
+\begin{align*}
+A^S_{gg} & = 1+ a_s A^{S(1)}_{gg} + a_s^2 A^{S(2)}_{gg} \\
+A^S_{gq} & = 1+ a_s^2 A^{S(2)}_{gq} \\
+A^S_{gH} & = 1+ a_s A^{S(1)}_{gH} \\
+A^S_{qg} & = 1 \\
+A^{ns}_{qq} &= 1+ a_s^2 A^{ns(1)}_{qq}  \\
+A^S_{Hg} & = 1+ a_s A^{S(1)}_{Hg} + a_s^2 A^{S(2)}_{Hg} \\
+A^{ps}_{Hq} & = 1+ a_s^2 A^{ps(2)}_{Hq} \\
+A_{HH} & = 1+ a_s A^{(2)}_{HH}
+\end{align*}
+
+\section{Basis Rotation}
+
+After the matching we have to perform the basis rotation from the basis with $\Sigma_{(n_f)}^{(n_f+1)}$, $\Sigma_{\Delta(n_f)}^{(n_f+1)}$, $h_+^{(n_f+1)}$ to the basis with $\Sigma_{(n_f+1)}^{(n_f+1)}$, $\Sigma_{\Delta(n_f+1)}^{(n_f+1)}$, $T_i^{(n_f+1)}$ (all the considerations that we will do for this basis rotation apply identically to the components $V$, $V_\Delta$, $h_-$, $V_i$). Being all the PDFs in the $n_f+1$ flavor scheme, from now on we will drop the superscript $(n_f+1)$.
+
+\subsection{$\Sigma$}
+For the $\Sigma$ component the basis rotation is very simple, being
+\begin{equation*}
+\Sigma_{(n_f+1)}=\Sigma_{(n_f)}+h_+
+\end{equation*}
+\subsection{$\Sigma_{\Delta}$}
+This component requires a bit more work: starting from
+\begin{equation*}
+\begin{cases}
+\Sigma &= \Sigma_u+\Sigma_d \\
+\Sigma_{\Delta} &= \frac{n_d}{n_u}\Sigma_u-\Sigma_d
+\end{cases}
+\end{equation*}
+we obtain that
+\begin{equation*}
+\begin{cases}
+\Sigma_u &= \frac{n_u}{n_f}\Bigl(\Sigma+\Sigma_{\Delta} \Bigr)\\
+\Sigma_d &= \frac{n_d}{n_f}\Sigma-\frac{n_u}{n_f}\Sigma_{\Delta}
+\end{cases}
+\end{equation*}
+
+Therefore, we find
+\begin{align*}
+\Sigma_{\Delta(n_f+1)} &= \frac{n_d(n_f+1)}{n_u(n_f+1)}\Sigma_{u(n_f)} - \Sigma_{d(n_f)} + k(n_f) h_+
+\end{align*}
+where
+\begin{equation*}
+k(n_f) =
+\begin{cases}
+\frac{n_d(n_f+1)}{n_u(n_f+1)} \quad \text{if h=up-like}\\
+-1  \quad \text{if h=down-like}
+\end{cases}
+\end{equation*}
+
+In the end we find that
+\begin{align*}
+\Sigma_{\Delta(n_f+1)} &= \frac{1}{n_f}\Bigl(\frac{n_d(n_f+1)}{n_u(n_f+1)}n_u(n_f)-n_d(n_f)\Bigr)\Sigma_{(n_f)} + \frac{n_f+1}{n_u(n_f+1)}\frac{n_u(n_f)}{n_f}\Sigma_{\Delta(n_f)} + k(n_f) h_+
+\end{align*}
+
+\subsection{$T_i$}
+In the end, we have to find the rotation for the $T_i$ component: being
+\begin{align*}
+T_3^d &=d^+ - s^+ = \Sigma_{d(n_f)} - h_+\quad \text{for $n_f=3$}\\
+T_3^u &=u^+ - c^+ =\Sigma_{u(n_f)} - h_+\quad \text{for $n_f=4$}\\
+T_8^d &=d^+ + s^+ - 2b^+ =\Sigma_{d(n_f)} - 2h_+\quad \text{for $n_f=5$}\\
+T_8^u &=u^+ + c^+ - 2t^+ =\Sigma_{u(n_f)} - 2h_+\quad \text{for $n_f=6$}
+\end{align*}
+
+Using the expressions of $\Sigma_u$ and $\Sigma_d$ as a function of $\Sigma$ and $\Sigma_\Delta$, we can write that
+\begin{equation*}
+T_i = f_1(n_f) \Sigma_{(n_f)} + f_2(n_f)\Sigma_{\Delta(n_f)} + f_3(n_f)h_+
+\end{equation*}
+with
+\begin{align*}
+f_1(n_f) &=
+\begin{cases}
+&\frac{n_u(n_f)}{n_f} \quad \text{if $h$ is up-like ($n_f$=3,5)} \\
+&\frac{n_d(n_f)}{n_f} \quad \text{if $h$ is down-like ($n_f$=2,4)}
+\end{cases} \\
+f_2(n_f) &=
+\begin{cases}
+&\frac{n_u(n_f)}{n_f} \quad \text{if $h$ is up-like} \\
+&-\frac{n_u(n_f)}{n_f} \quad \text{if $h$ is down-like}
+\end{cases} \\
+f_3(n_f) &=
+\begin{cases}
+&-1\quad \text{if $h$ is $s$, $c$  ($n_f$=2,3)} \\
+&-2 \quad \text{if $h$ is $b$, $t$  ($n_f$=4,5)}
+\end{cases}
+\end{align*}
+
+\end{document}