Higher-Order Sudakov Resummation in Coupled Gauge Theories

TL;DR

This paper develops a method for higher-order resummation of Sudakov logarithms in theories with multiple coupled gauge interactions, providing precise calculations relevant for QCD, QED, and electroweak processes.

Contribution

It introduces strategies to solve coupled evolution equations and presents the complete three-loop Sudakov evolution factor for QCD and QED, highlighting the importance of exact calculations.

Findings

Complete three-loop (NNLL) QCD×QED Sudakov evolution factor obtained.

Approximate analytic expressions can cause several percent errors at low scales.

Numerical differences between methods can exceed perturbative precision, affecting high-precision analyses.

Abstract

We consider the higher-order resummation of Sudakov double logarithms in the presence of multiple coupled gauge interactions. The associated evolution equations depend on the coupled $\beta$ functions of two (or more) coupling constants $\alpha_a$ and $\alpha_b$, as well as anomalous dimensions that have joint perturbative series in $\alpha_a$ and $\alpha_b$. We discuss possible strategies for solving the system of evolution equations that arises. As an example, we obtain the complete three-loop (NNLL) QCD$\otimes$QED Sudakov evolution factor. Our results also readily apply to the joint higher-order resummation of electroweak and QCD Sudakov logarithms. As part of our analysis we also revisit the case of a single gauge interaction (pure QCD), and study the numerical differences and reliability of various methods for evaluating the Sudakov evolution factor at higher orders. We find that the approximations involved in deriving commonly used analytic expressions for the evolution kernel can induce noticeable numerical differences of several percent or more at low scales, exceeding the perturbative precision at N$^3$LL and in some cases even NNLL. Therefore, one should be cautious when using approximate analytic evolution kernels for high-precision analyses.

Equations

$$\alpha (\mu {)}^n{\ln }^m\frac{\mu }{{\mu }_0}\text{with}m\le 2n.$$

$$\displaystyle U(\mu_{0},\mu)=\exp\biggl{\{}\int_{\mu_{0}}^{\mu}\frac{d\mu^{\prime}}{\mu^{\prime}}\Gamma_{\mathrm{cusp}}[\alpha_{a}(\mu^{\prime}),\alpha_{b}(\mu^{\prime})]\ln\frac{Q}{\mu^{\prime}}+\gamma[\alpha_{a}(\mu^{\prime}),\alpha_{b}(\mu^{\prime})]\biggr{\}}\,.$$

$$U(\mu ,{\mu }_0)U({\mu }_0,\mu )\stackrel{!}{=}1.$$

$$\frac{\mathrm{d}\alpha (\mu )}{\mathrm{d}\ln \mu }\equiv \beta [\alpha (\mu )]$$

$$\displaystyle=-2\beta_{0}\,\frac{\alpha(\mu)^{2}}{4\pi}\,\biggl{[}1+\epsilon\,\frac{\alpha(\mu)}{4\pi}\,b_{1}+\epsilon^{2}\,\frac{\alpha(\mu)^{2}}{(4\pi)^{2}}\,b_{2}+\mathcal{O}(\epsilon^{3})\biggr{]}\,.$$

$$b_n=\frac{{\beta }_n}{{\beta }_0}.$$

$$\alpha (\mu )=\frac{\alpha ({\mu }_0)}{X},X\equiv X({\mu }_0,\mu )=1+\frac{\alpha ({\mu }_0)}{2\pi }{\beta }_0\ln \frac{\mu }{{\mu }_0},$$

$$\frac{\mathrm{d}\alpha (\mu )}{\mathrm{d}\ln \mu }$$

$$\displaystyle=-2\beta_{0}\frac{\alpha(\mu)^{2}}{4\pi}\biggl{[}1+\epsilon\,\frac{\alpha(\mu_{0})}{4\pi}\,\frac{b_{1}}{X(\mu_{0},\mu)}+\mathcal{O}(\epsilon^{2})\biggr{]}\,.$$

$$\frac{1}{\alpha (\mu )}$$

$$\Rightarrow \alpha (\mu )$$

$$\alpha (\mu )$$

$$\frac{\mathrm{d}\alpha (\mu )}{\mathrm{d}\ln \mu }$$

$$\displaystyle=-2\beta_{0}\,\frac{\alpha(\mu_{0})^{2}}{4\pi}\,\frac{1}{X^{2}}\biggl{\{}\biggl{[}1-2\epsilon\,\frac{\alpha(\mu_{0})}{4\pi}\,b_{1}\,\frac{\ln X}{X}\biggl{]}+\epsilon\,\frac{\alpha(\mu_{0})}{4\pi}\,\frac{b_{1}}{X}\biggr{\}}\,.$$

$$\frac{\mathrm{d}\alpha (\mu )}{\mathrm{d}\ln \mu }$$

$$\displaystyle=-2\beta_{0}\frac{\alpha(\mu)^{2}}{4\pi}\biggl{[}1+\epsilon\,\frac{\alpha(\mu_{0})}{4\pi}\,\frac{b_{1}}{X}+\epsilon^{2}\,\frac{\alpha(\mu_{0})^{2}}{(4\pi)^{2}}\frac{b_{2}-b_{1}^{2}\,\ln X}{X^{2}}+\mathcal{O}(\epsilon^{3})\biggr{]}\,.$$

$$\alpha (\mu )$$

$$\displaystyle\quad+\epsilon^{2}\,\frac{\alpha(\mu_{0})^{2}}{(4\pi)^{2}}\,\frac{1}{X^{2}}\Bigl{[}b_{2}(1-X)+b_{1}^{2}\,(\ln^{2}X-\ln X-1+X)\Bigr{]}+\mathcal{O}(\epsilon^{3})\biggr{\}}\,.$$

$$\frac{\mathrm{d}{\alpha }_a(\mu )}{\mathrm{d}\ln \mu }$$

$$\frac{\mathrm{d}{\alpha }_b(\mu )}{\mathrm{d}\ln \mu }$$

$$b_{nk}^a=\frac{{\beta }_{nk}^a}{{\beta }_{00}^a},b_{nk}^b=\frac{{\beta }_{nk}^b}{{\beta }_{00}^b},\text{and}{b}_0^a=\frac{{\beta }_{00}^a}{{\beta }_{00}^b},b_0^b=\frac{{\beta }_{00}^b}{{\beta }_{00}^a}.$$

$${\alpha }_a(\mu )$$

$${\alpha }_b(\mu )$$

$$\frac{\mathrm{d}{\alpha }_a(\mu )}{\mathrm{d}\ln \mu }$$

$$\displaystyle=-2\beta_{00}^{a}\frac{\alpha_{a}(\mu)^{2}}{4\pi}\biggl{[}1+\epsilon_{a}\,\frac{\alpha_{a}(\mu_{0})}{4\pi}\,\frac{b^{a}_{10}}{X_{a}(\mu_{0},\mu)}+\epsilon_{b}\,\frac{\alpha_{b}(\mu_{0})}{4\pi}\,\frac{b^{a}_{01}}{X_{b}(\mu_{0},\mu)}+\mathcal{O}(\epsilon_{a}^{2},\epsilon_{a}\epsilon_{b},\epsilon_{b}^{2})\biggr{]}\,.$$

$$\displaystyle\alpha_{a}(\mu)=\alpha_{a}(\mu_{0})\biggl{[}X_{a}+\frac{\alpha_{a}(\mu_{0})}{4\pi}\Bigl{(}\epsilon_{a}\,b_{10}^{a}\,\ln X_{a}+\epsilon_{b}\,b^{a}_{0}\,b_{01}^{a}\,\ln X_{b}\Bigr{)}\biggr{]}^{-1}\,.$$

$$\frac{{\alpha }_a({\mu }_0)}{{\alpha }_a(\mu )}$$

$$\displaystyle\quad+\epsilon_{b}\,\frac{\alpha_{a}(\mu_{0})}{4\pi}\,b^{a}_{0}\,b_{01}^{a}\ln X_{b}+\epsilon_{b}^{2}\,\frac{\alpha_{a}(\mu_{0})\alpha_{b}(\mu_{0})}{(4\pi)^{2}}b^{a}_{0}\biggl{(}b_{02}^{a}\,\frac{X_{b}-1}{X_{b}}+b_{01}^{a}b_{10}^{b}\,\frac{1-X_{b}+\ln X_{b}}{X_{b}}\biggr{)}$$

$$\displaystyle\quad+\epsilon_{a}\epsilon_{b}\,\frac{\alpha_{a}(\mu_{0})}{b_{0}^{a}\,\alpha_{a}(\mu_{0})-\alpha_{b}(\mu_{0})}\biggl{[}\frac{\alpha_{a}^{2}(\mu_{0})}{(4\pi)^{2}}(b^{a}_{0})^{2}b_{10}^{a}b_{01}^{a}\biggl{(}\frac{X_{b}}{X_{a}}\ln X_{b}-\frac{1-X_{b}}{1-X_{a}}\ln X_{a}\biggr{)}$$

$$\displaystyle\qquad-\frac{\alpha_{b}^{2}(\mu_{0})}{(4\pi)^{2}}\,b_{01}^{a}b_{01}^{b}\biggl{(}\frac{X_{a}}{X_{b}}\ln X_{a}-\frac{1-X_{a}}{1-X_{b}}\ln X_{b}\biggr{)}+\frac{\alpha_{a}(\mu_{0})\alpha_{b}(\mu_{0})}{(4\pi)^{2}}b^{a}_{0}\,b_{11}^{a}\ln\frac{X_{a}}{X_{b}}\biggr{]}\,.$$

$$\frac{\mathrm{d}F(\mu )}{\mathrm{d}\ln \mu }={\gamma }_F(\mu )\otimes F(\mu ),$$

$${\gamma }_F(\mu )=\frac{1}{ϵ}{\Gamma }_{\mathrm{cusp}}[\alpha (\mu )]\ln \frac{Q}{\mu }+\gamma [\alpha (\mu )],$$

$$\Gamma_{\mathrm{cusp}}(\alpha)=\sum_{n=0}^{\infty}\epsilon^{n+1}\,\Gamma_{n}\,\Bigl{(}\frac{\alpha}{4\pi}\Bigr{)}^{n+1}\,,\qquad\gamma(\alpha)=\sum_{n=0}^{\infty}\epsilon^{n+1}\,\gamma_{n}\,\Bigl{(}\frac{\alpha}{4\pi}\Bigr{)}^{n+1}\,.$$

$$F(\mu )=F({\mu }_0)U({\mu }_0,\mu ),$$

$$\displaystyle U(\mu_{0},\mu)=\exp\biggl{\{}\int_{\mu_{0}}^{\mu}\frac{\mathrm{d}\mu^{\prime}}{\mu^{\prime}}\gamma_{F}(\mu^{\prime})\biggr{\}}=\exp\biggl{\{}\int_{\mu_{0}}^{\mu}\frac{\mathrm{d}\mu^{\prime}}{\mu^{\prime}}\,\frac{1}{\epsilon}\,\Gamma_{\mathrm{cusp}}[\alpha(\mu^{\prime})]\ln\frac{Q}{\mu^{\prime}}+\gamma[\alpha(\mu^{\prime})]\biggr{\}}\,.$$

Full text

As part of our analysis we also revisit the case of a single gauge interaction (pure QCD), and study the numerical differences and reliability of various methods for evaluating the Sudakov evolution factor at higher orders. We find that the approximations involved in deriving commonly used analytic expressions for the evolution kernel can induce noticeable numerical differences of several percent or more at low scales, exceeding the perturbative precision at N3LL and in some cases even NNLL. Therefore, one should be cautious when using approximate analytic evolution kernels for high-precision analyses.

1 Introduction

In perturbative quantum field theories it is well known that observables sensitive to physics at different scales $\mu\ll\mu_{0}$ can have their perturbative expansions in the coupling constant $\alpha$ generically enhanced by Sudakov logarithms of the form

[TABLE]

For sufficiently separated scales, the logarithms grow large enough to dominate (and eventually deteriorate) the perturbative series. The reliability and precision of theoretical predictions can therefore be improved (or restored) by a reorganization of the perturbative series into a form that keeps the highest-power logarithms to all orders in $\alpha$, a procedure called resummation. Various different formalisms have been developed to achieve the resummation, with strongly interacting processes being the most widely studied due to their dominant contribution to many key collider processes. Indeed, multiple observables have been resummed to next-to-next-to-leading logarithmic (NNLL) and even N3LL accuracy within QCD, in an effort to match the ever increasing experimental precision.

Despite their smaller couplings in the Standard Model (SM), corrections from the emission of electroweak (EW) bosons can be comparable to those of QCD calculated at NNLO (cf. $\alpha_{e}\sim\alpha_{s}^{2}$). Furthermore, the exchange of massive virtual EW bosons in high-energy processes can generate EW Sudakov logarithms of the form of eq. (1.1), which can cause sizeable EW corrections. The resummation of EW Sudakov logarithms has been studied for many years, albeit typically at lower orders than in QCD, see e.g. refs. [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14]. As a result, achieving sufficiently precise predictions for many collider observables requires considering them in a joint QCD$\otimes$EW environment to fully capture all relevant effects. This is of course true when considering an extremely high energy future collider [15], where EW corrections can be $\mathcal{O}(1)$, but also for measurements at the LHC reaching percent-level precision, the prime example being the high-precision measurements of $W$ and $Z$ production [16, 17, 18, 19, 20, 21]. The recent literature reflects this: The impact of EW and mixed QCD$\otimes$EW corrections on the $W$-mass measurement have received much attention (see e.g. refs. [22, 23, 24, 25] and references therein). QED corrections to the evolution of parton distribution functions (PDFs) have also been obtained, see e.g. refs. [26, 27, 28, 29, 30, 31]. The full NNLO $\mathcal{O}(\alpha_{e}\alpha_{s})$ mixed QCD$\otimes$QED corrections and $\mathcal{O}(\alpha_{e}^{2})$ QED corrections for on-shell $Z$ production were calculated recently in ref. [32] and for Higgs production in bottom-quark annihilation in ref. [33]. The one-loop QED corrections to the Sudakov resummation were included in the high-precision analysis of thrust in $e^{+}e^{-}$ collisions [34]. The resummed $p_{T}$ spectrum of $Z$-boson production including QED corrections was obtained in ref. [35], capturing the pure QED logarithmic contributions at NLL and the mixed QCD$\otimes$QED contributions at LL. Closely related, the QED corrections to the two-loop anomalous dimensions of $p_{T}$-dependent distributions were obtained in ref. [36] and of the quark form factors in ref. [33].

In this paper, we step back and reevaluate some of the technical aspects of Sudakov resummation when the interactions of two gauge symmetries $G_{a}\otimes G_{b}$ are involved, staying agnostic as to the precise resummation formalism utilized. In particular, we analyze the integrand structure of the Sudakov evolution factor

[TABLE]

An evolution factor of this form necessarily appears in all formulations of higher-order Sudakov resummation. It implicitly depends on the $\beta$ functions controlling the renormalization group evolution (RGE) of both gauge coupling constants $\alpha_{a}(\mu)$ and $\alpha_{b}(\mu)$, which in general are a coupled system of differential equations, as well as anomalous dimensions ($\Gamma_{\mathrm{cusp}},\gamma$) whose perturbative expansions are themselves joint series in $\alpha_{a,b}$. We attempt to be as generic as possible in our discussion and therefore do not immediately specify $G_{a,b}$. To draw some phenomenological conclusions we will eventually consider the example of QCD$\otimes$QED, i.e. $G_{a}\equiv SU(3)_{c}$ and $G_{b}\equiv U(1)_{\rm em}$, for which we obtain the complete three-loop (NNLL) Sudakov evolution factor.

Interestingly, we also find that the approximations made in obtaining closed-form analytic expressions for eq. (1.2) that are commonly used in the literature can lead to nonnegligible numerical differences. When evolving to low scales, where the resummation becomes most important, the resulting effects can reach several percent or more. In other words, the use of different strategies for evaluating eq. (1.2) can be a source of nontrivial systematic differences between different resummation implementations that can potentially exceed the perturbative precision one is aiming for. As one manifestation, we find nontrivial violations of the consistency of the evolution, which we probe by testing the closure condition

[TABLE]

Since these issues already appear for a single gauge interaction we devote a substantial fraction of our paper to exploring them in this simpler case, using pure QCD as test case. We consider various combinations of numerical and approximate analytic treatments of the Sudakov factor and study their accuracy. Ultimately, we are forced to conclude that the commonly used approximate analytic expressions for eq. (1.2) are not sufficiently reliable when aiming for percent-level precision. We find that a seminumerical approach, where the $\mu$-integration in the exponent of eq. (1.2) is carried out numerically while using an approximate analytic solution for the running of the coupling, provides a good compromise, which also straightforwardly generalizes to multiple gauge interactions.

The structure of the paper is as follows. In section 2, we discuss the coupled $\beta$ functions for $\alpha_{a,b}$. We first review the standard approximate analytic solutions for $\alpha$ in the one-dimensional case, and then derive corresponding approximate analytic solutions for the coupled system up to three-loop (NNLO) running. In section 3, we move to analyzing the Sudakov evolution kernels for the one-dimensional case, discussing several methods for evaluating it and studying their numerical performance up to N3LL. In section 4 we then apply the lessons learned to the case of two coupled gauge interactions. We conclude in section 5. In appendix A, we provide the perturbative ingredients utilized in our numerical comparisons, including the complete set of three-loop QCD$\otimes$QED coefficients. In appendix B, we discuss the extraction of the complete three-loop mixed QCD$\otimes$QED $\beta$-function coefficients.

2 Iterative solutions to $\beta$-function RGEs

Before discussing the evaluation of the Sudakov evolution kernel, it will be important to analyze one of its key ingredients, namely the solution of the $\beta$-function RGE of the coupling constant. In subsection 2.1, we review the well-known case of a single gauge theory, with particular attention given to the numerical accuracy of different approximate analytic RGE solutions. In subsection 2.2 we then discuss the solution of the coupled RGE system for the coupling constants of two gauge theories.

2.1 Single gauge theory

We start from the well-known $\beta$-function RGE for the coupling constant $\alpha(\mu)$ of a generic gauge theory,

[TABLE]

Here, $\mu$ is the renormalization scale, and we introduced a formal expansion parameter $\epsilon\equiv 1$, which we use to keep track of the evolution order. The perturbative coefficients $b_{n}$ in the second line are defined by the ratio

[TABLE]

We take the viewpoint that eq. (2.1) defines the running order of the RGE. That is, keeping the terms up to $\mathcal{O}(\epsilon^{k})$ in eq. (2.1) defines the NkLO or ($k+1$)-loop running of $\alpha(\mu)$. At leading order, $\mathcal{O}(\epsilon^{0})$, eq. (2.1) has the well-known exact analytic solution,

[TABLE]

where $\alpha(\mu_{0})$ is a boundary condition for the coupling constant.

As is well known, eq. (2.1) does not admit an exact analytic solution at NLO and beyond.111More precisely, eq. (2.1) can still be integrated analytically at NLO and even NNLO. The resulting expressions, however, cannot be analytically solved for $\alpha(\mu)$ in terms of $\alpha(\mu_{0})$ anymore. While the exact solution can be easily obtained numerically using standard numerical differential-equation solvers, in practice it is often more convenient to have an approximate analytic solution that can be evaluated much faster than the numerical solution, which becomes important when the scale $\mu$ is not fixed but dynamical. This is precisely the case for the Sudakov evolution kernel, for which we will need to integrate $\alpha(\mu)$ over $\mu$.

In what follows, we review an iterative method to obtain an approximate analytic solution for eq. (2.1). At NLO, $\mathcal{O}(\epsilon)$, the $\beta$-function RGE reads

[TABLE]

In the second line we substituted the LO solution in eq. (2.3) for the $\mu$ dependence of $\alpha(\mu)$ in the $\mathcal{O}(\epsilon)$ term. This induces an $\mathcal{O}(\epsilon^{2})$ error, since the difference between the LO and NLO $\mu$ dependence will itself be of $\mathcal{O}(\epsilon)$. Since the $\mu$ dependence in the term in square brackets is now explicit, it can be easily integrated, yielding the NLO solution

[TABLE]

We refer to this (and its higher-order analogues) as the “iterative” solution for $\alpha(\mu)$. We can also expand the inverse in eq. (2.1) in $\epsilon$ to obtain

[TABLE]

We will refer to eq. (2.6) (and its higher-order analogues) as the “expanded” solution.

One might wonder whether the iterative or expanded solution provides a better approximation to the exact solution, and to this end it is instructive to see to what extent they satisfy the original $\beta$-function RGE. For the iterative solution in eq. (2.1), it is trivial to verify that upon differentiation with respect to $\ln\mu$ it reproduces the RGE as given in the second line of eq. (2.1). On the other hand, taking the $\ln\mu$ derivative of the expanded NLO solution in eq. (2.6) yields

[TABLE]

Comparing this with eqs. (2.1) and (2.6), we see that the term in square brackets corresponds to expanding the overall $\alpha(\mu)^{2}$ in the $\beta$ function to $\mathcal{O}(\epsilon)$. This clearly amounts to a further approximation, so we can expect the expanded solution in general to provide a worse approximation, which is indeed what we will find numerically below.

The iterative solution at NNLO, $\mathcal{O}(\epsilon^{2})$, follows analogously. Starting from the exact NNLO RGE, we substitute the (expanded) NLO and LO solutions, eqs. (2.6) and (2.3), in the $\mathcal{O}(\epsilon)$ and $\mathcal{O}(\epsilon^{2})$ terms, keeping all terms up to $\mathcal{O}(\epsilon^{2})$ as well as the overall $\alpha(\mu)^{2}$,

[TABLE]

The $\mu$-dependence in the square brackets on the second line, encoded in $X\equiv X(\mu_{0},\mu)$, is again simple enough to be integrated analytically.222This would not be the case had we used the unexpanded NLO solution from eq. (2.1) in the $\mathcal{O}(\epsilon)$ term. This yields the NNLO iterative solution

[TABLE]

Expanding the inverse in eq. (2.9) in $\epsilon$, we obtain the expanded NNLO solution,

[TABLE]

It is straightforward to extend the iterative solution to higher orders. To obtain the NkLO solution, one simply inserts the Nk-1LO solution for $\alpha(\mu)$ into the NkLO $\beta$ function and expands it to $\mathcal{O}(\epsilon^{k})$, while keeping the overall $\alpha(\mu)^{2}$ exact. The N3LO solution is given in appendix A.1.

To illustrate the numerical precision of the approximate analytic solutions, we take the QCD coupling constant $\alpha_{s}(\mu)$ as an example, using $\alpha_{s}(m_{Z})=0.118$ as our boundary condition. As we are primarily interested in the numerical precision of the solution, we always use $n_{f}=5$ massless flavors and do not consider any flavor thresholds. The difference of the iterative and expanded solutions to the exact333We always perform the numerical solution with sufficiently high numerical precision such that the numerical error is completely negligible for our purposes. numerical solution, which we refer to as the approximation error, is shown in figure -899 for different running orders.

The approximation error decreases as the order increases, as expected. The approximation error is largest for running from $m_{Z}$ down to lower scales, since here the running increases the coupling. The iterative solution still provides an excellent approximation, with the error at NNLO and beyond reaching at most 0.1% when running down to $\mu=2\,\mathrm{GeV}$, and at most 0.3% at $\mu=1\,\mathrm{GeV}$. For running above $m_{Z}$, the approximation error is much smaller due to asymptotic freedom. Beyond the highest scale shown, $\mu=10^{4}\,\mathrm{GeV}$, the error stops growing and at some point starts decreasing again. We also observe that the approximation error for the expanded solution (gray lines) is always 2-3 times larger than for the iterative solution.

2.2 Two coupled gauge theories

We now consider the case of two gauge theories. Their $\beta$ functions become coupled as soon as there are matter fields that are charged under both gauge interactions, since loops of matter particles can exchange the gauge bosons from both theories. The example we will eventually consider is the mixed QCD$\otimes$QED running. For now, we keep the discussion general and consider the following set of coupled $\beta$-function RGEs for two gauge couplings $\alpha_{a}(\mu)$ and $\alpha_{b}(\mu)$,

[TABLE]

We have again introduced formal expansion parameters $\epsilon_{a,b}\equiv 1$ to easily keep track of the evolution order. For future convenience, we also define the rescaled coefficients

[TABLE]

Note that by definition $\beta^{x}_{n0}\equiv\beta_{n}$ and $b^{x}_{n0}\equiv b_{n}$ are the coefficients of the individual gauge theories in the absence of the second.

As before, the order in $\epsilon_{a,b}$ to which eq. (2.2) is expanded is what defines the running order of the couplings. Generically, we consider $\epsilon_{a,b}$ on equal footing and define the NnLO evolution by including in eq. (2.2) all terms of $\mathcal{O}(\epsilon_{a}^{k}\epsilon_{b}^{n-k})$ with $0\leq k\leq n$. Note that this also makes it easy to have well-defined mixed orders, e.g., when there is a hierarchy between the two couplings, as is the case for QCD and QED. For this purpose, one can simply specify the explicit combinations of powers of $\epsilon_{a}$ and $\epsilon_{b}$ that are included in eq. (2.2).

As in the single gauge scenario, an exact solution to the coupled system of differential equations can be obtained straightforwardly at any given order by solving it numerically. Our goal is to derive an approximate analytic solution for the coupled $\beta$-function RGE system by extending the iterative method in subsection 2.1. The key property of eq. (2.2) that allows us to do so is that, at leading order, $\mathcal{O}(\epsilon_{a}^{0},\epsilon_{b}^{0})$, the RGE system decouples, yielding exact LO solutions:

[TABLE]

To obtain an approximate NLO solution, it then suffices to substitute the above LO solutions into the $\mathcal{O}(\epsilon_{a})$ and $\mathcal{O}(\epsilon_{b})$ terms of the NLO RGE system, which induces $\mathcal{O}(\epsilon_{a}^{2},\epsilon_{a}\epsilon_{b},\epsilon_{b}^{2})$ errors,

[TABLE]

The terms in square brackets can now be explicitly integrated over $\ln\mu$ to obtain the iterative NLO solution

[TABLE]

The solution for $\alpha_{b}(\mu)$ is given by replacing $a\leftrightarrow b$ everywhere.

Note that in eq. (2.15) the utility of the $\epsilon_{a,b}$ counting parameters becomes evident. Naively, one might have expected that the $\mathcal{O}(\epsilon_{a})$ and $\mathcal{O}(\epsilon_{b})$ terms will be proportional to $\alpha_{a}(\mu_{0})$ and $\alpha_{b}(\mu_{0})$ respectively, which however is not the case, as both are proportional to $\alpha_{a}(\mu_{0})$. Instead, the $\epsilon_{a,b}$ actually keep track of the fact that the $\ln X_{a,b}$ factors respectively resum a series of $\alpha_{a,b}^{n}\ln^{n}(\mu/\mu_{0})$ terms.

The iterative NNLO solution is obtained by substituting the NLO solutions into the $\beta$-function RGE system and expanding it to $\mathcal{O}(\epsilon_{a}^{2},\epsilon_{a}\epsilon_{b},\epsilon_{b}^{2})$, while keeping the overall $\alpha_{a,b}(\mu)^{2}$ exact. We find

[TABLE]

As before, the solution for $\alpha_{b}(\mu)$ is obtained by replacing $a\leftrightarrow b$. The terms in the first line correspond to the NNLO solution of $\alpha_{a}$ in the absence of $\alpha_{b}$, while the remaining ones are the mixing contributions involving at least one power of $\epsilon_{b}$. The corresponding expanded solution is obtained by inverting eq. (2.2) and expanding it in $\epsilon_{i}$. Note that when doing so, it becomes essential to expand in terms of $\epsilon_{i}$ and not $\alpha_{i}(\mu_{0})$. The necessity of using $\epsilon_{i}$ as expansion parameter is also evident in the mixed $\mathcal{O}(\epsilon_{a}\epsilon_{b})$ term in square brackets, which involves a nontrivial rational function of both couplings. The above iterative NNLO solution will be a key ingredient in the seminumerical evaluation of the Sudakov evolution factor in section 4.

To illustrate the approximation error for the two analytic solutions, we consider the case of QCD$\otimes$QED. The relevant coefficients for the coupled $\beta$-function RGE to NNLO are given in appendices A.2 and A.3. Curiously, we were not able to find explicit expressions in the literature for the mixed three-loop QED coefficient $\beta_{11}^{e}$. We therefore performed an explicit extraction of all mixed three-loop coefficients from the general results for a generic product group given in ref. [37], as discussed in more detail in appendix B.

For the numerical results, we use $\alpha_{s}(m_{Z})=0.118$ and $\alpha_{e}(m_{Z})=1/127$ as boundary conditions, $n_{f}=5$ for the number of active quark flavors, and $n_{\ell}=3$ for the number of active charged leptons. As before we do not consider any flavor thresholds.

The approximation error of the iterative (expanded) solution relative to the exact numerical solution is shown in figure -898 by the colored (gray) lines at NLO (dashed) and NNLO (solid). For the strong coupling constant (left panel), the approximation error for the iterative solution does not exceed 1%, and it is again 2-3 times larger for the expanded solution. For the QED coupling constant (right panel), the approximation error is much smaller owing to the fact that $\alpha_{e}$ is much smaller, and therefore both the iterative and expanded solutions yield equally good approximations.

It is interesting to see the individual effects of various terms in the coupled $\beta$ function. In figures -897 and -896 we show them for the case of QCD$\otimes$QED as well as for a toy QCD$\otimes$QED*′* for which we set the boundary condition to $\alpha_{e}^{\prime}(m_{Z})=1/10$. In figure -897, as expected, the QED corrections $\sim\mathcal{O}(\epsilon^{0}_{s},\epsilon^{1}_{e})$ have almost negligible effect on $\alpha_{s}$, with only the higher order QCD corrections $\sim\mathcal{O}(\epsilon_{s}^{1},\epsilon_{e}^{0})$ being relevant. This is also the case for the QED coupling constant, whose evolution is almost entirely dictated by the QCD corrections. But this is not the case for the toy QCD$\otimes$QED*′* scenario shown in figure -896. The fact that the QCD coupling constant is still almost unaffected by the QED*′* corrections is somewhat accidental and due to the fact that the $\beta$-function coefficients $\beta^{s}_{10}$ and $\beta^{s}_{01}$ differ numerically by an order of magnitude. On the other hand, the QED*′* coupling constant now has comparable higher-order QCD and QED*′* corrections since here both $\beta$-function coefficients are of similar numerical size.

3 Sudakov evolution kernels with a single gauge interaction

In this section, we examine different strategies for evaluating the Sudakov evolution kernel for the case of a single gauge interaction, with particular emphasis on their numerical accuracy and reliability. This will then serve as a guide when considering the two-dimensional case in section 4.

In subsection 3.1, we start with a general discussion and give an overview of the different methods we will study for evaluating the Sudakov evolution factor, which are then elaborated on in subsection 3.2 to subsection 3.5. We then provide a detailed numerical comparison in subsection 3.6.

3.1 General overview

One way to systematically resum Sudakov logarithms is based on factorizing the perturbative series. The relevant factorization can be derived diagrammatically or using effective field theories (EFTs). All ingredients of the factorized cross section obey (renormalization group) evolution equations of the form

[TABLE]

where $\gamma_{F}(\mu)$ is the anomalous dimension, and the function $F(\mu)$ is any one of the factorized ingredients. The Sudakov resummation is then performed by solving the resulting coupled system of RGE equations.

In general, $\gamma_{F}(\mu)$ and $F(\mu)$ both depend on an additional external kinematic quantity (which appears as part of the argument of the Sudakov logarithms in the cross section). The $\otimes$ denotes the fact that $\gamma_{F}$ and $F$ are not necessarily multiplied but could be convolved in said kinematic variable. It is worth nothing that the convolution structure can play a significant role in the solution of eq. (3.1) (for a recent detailed discussion see e.g. ref. [38]). However, it does not play any role for the purpose of our discussion, so we consider only the simplest multiplicative case. In case of a convolution, one can always transform to a suitable conjugate space (e.g. Fourier or Laplace space), where the convolution turns into a simple product. The Sudakov evolution factor in that conjugate space then has the same general form we discuss here and all our conclusions apply equally.

The all-order expansion of the anomalous dimension is given by

[TABLE]

where $Q$ denotes the above-mentioned kinematic quantity, $\Gamma_{\mathrm{cusp}}$ is (proportional to) the cusp anomalous dimension, and $\gamma$ is the noncusp anomalous dimension. They obey the perturbative expansions

[TABLE]

We have again introduced the formal expansion parameter $\epsilon\equiv 1$, which we will use to define the resummation order. Since the $\mu$ dependence of $\gamma_{F}(\mu)$ primarily enters via the coupling constant’s $\mu$ dependence, the $\beta$ RGE for $\alpha(\mu)$ in eq. (2.1) is an integral part of the full RGE system to be solved. In particular, the $\epsilon$ parameter in eqs. (3.2) and (3.3) is the same that appears in eq. (2.1).

As for the case of the coupling constant before, the truncation of eq. (3.2) together with eq. (2.1) to a certain order in $\epsilon$ fundamentally defines the resummation order. That is, keeping terms up to $\mathcal{O}(\epsilon^{k})$ defines the Sudakov evolution at NkLL order. The explicit $1/\epsilon$ factor for the cusp term accounts for the fact that it comes with an additional explicit logarithm relative to the noncusp term. As a result, the noncusp term always enters at one lower order in perturbation theory than the cusp term. And since the $\beta$ function in eq. (2.1) starts at $\mathcal{O}(\epsilon^{0})$, it enters at the same loop order as the cusp anomalous dimension. So as usual, at NkLL order, we require the $k+1$-loop cusp and beta function coefficients and the $k$-loop noncusp coefficients.

Solving the RGE in eq. (3.1), one finds

[TABLE]

where $U(\mu_{0},\mu)$ is the Sudakov evolution factor given by

[TABLE]

It resums the Sudakov logarithms appearing in the perturbative series of $F(\mu)/F(\mu_{0})$. One can easily check that counting powers of $\epsilon$ in the anomalous dimension is equivalent to counting powers of logarithms in the Sudakov exponent, simply because the entire structure of $U(\mu_{0},\mu)$ is fully encoded by the anomalous dimension. In the final resummed cross section, where $F(\mu)$ is combined with other ingredients, $U(\mu_{0},\mu)$ eventually appears with both scales $\mu$ and $\mu_{0}$ corresponding to two different kinematic quantities such that $U(\mu_{0},\mu)$ resums (part of) the Sudakov logarithms of the ratio of these quantities.

We stress that while we obtained eq. (3.5) starting from the RGE in eq. (3.1), a Sudakov evolution factor of the same structure (necessarily) appears in all the various approaches for performing Sudakov resummation that exist in the literature. This includes EFT-based and non-EFT-based approaches, and both analytic as well as numerical Monte-Carlo techniques such as parton showers.

On the other hand, different implementations tend to follow different strategies for evaluating the integral in the Sudakov exponent. In the following sections we investigate several methods for doing so, paying close attention to where additional assumptions and/or approximations are made. As we will see, additional approximations that may appear mathematically justified can still conspire to yield results for eq. (3.5) that exhibit nontrivial numerical differences. We investigate the following methods in the sections that follow:

• •

Numerical: In this method, both the $\beta$-function RGE eq. (2.1) and the evolution kernel eq. (3.5) are evaluated fully numerically (with sufficiently high numerical precision that numerical integration errors are negligible). This provides the exact solution of the complete Sudakov RGE system at a given resummation order as defined above, and we will use it as the benchmark to compare the other methods against. As this method can be computationally expensive, it is often not very suitable for practical purposes.

• •

Seminumerical: One option to speed up the fully numerically method is to employ an approximate analytic solution to the $\beta$-function RGE, but to perform the kernel integration numerically. In other words, we numerically integrate eq. (3.5) but with the iterative analytic solution for $\alpha(\mu)$ used in the perturbative expansion of the anomalous dimensions. This is described along with the fully numerical method in subsection 3.2.

• •

Unexpanded analytic: Approximate but fully analytic, closed-form expressions for eq. (3.5) can be obtained. One way to achieve this is to exploit the $\beta$ function to turn the $\mu$ integration into an integration over $\alpha$, which can be performed analytically after some expansion in $\epsilon$. This is combined with the analytic iterative solution for $\alpha_{s}(\mu)$ as input. We give details of the derivation and the resulting forms explicitly in subsection 3.3.

• •

Expanded analytic: Another way to obtain an approximate closed-form result for eq. (3.5) is to insert the analytic solution for $\alpha(\mu)$ in eq. (2.1) in the perturbative expansions of $\Gamma(\alpha)$ and $\gamma(\alpha)$ and fully expanding the integrand in $\epsilon$. The integration can then be performed directly in terms of $\mu$. This method is described in more detail in subsection 3.4.

• •

Reexpanded analytic: Finally, to connect to some of the available literature, we elaborate another analytic approach where, upon achieving the expanded analytic evolution kernels, one further expands $\alpha(\mu_{0})$ in terms of $\alpha(\mu_{R})$ at a different reference scale $\mu_{R}$, assuming that $\mu_{0}\sim\mu_{R}$, i.e., that there is no hierarchy between them. This method is discussed in more detail in subsection 3.5.

Given these five methods of integration, we will probe their reliability in two different ways in subsection 3.6: closure tests and approximation errors. The latter simply tests the absolute difference between any one method and the numerically exact method described above. The closure test provides a test of the mathematical consistency of the evolution factor. That is, it should satisfy the renormalization group property

[TABLE]

which is obvious from its definition in eq. (3.5). A simple way to test that eq. (3.6) is satisfied is to consider the special case of $\mu_{2}=\mu_{0}$, which yields

[TABLE]

and simply expresses the fact that the RG evolution should close on itself. Since the intermediate scale here is completely arbitrary, this property should be satisfied exactly at any given order. However, deviations from unity can arise due to simplifying assumptions or approximations made in evaluating the integral in the exponent.

3.2 Numerical and seminumerical methods

The most accurate method of integration is to perform the integration fully numerically. The error introduced by the numerical integration routine can be made arbitrarily small at the expense of computing time. We always use a sufficiently high integration precision that the numerical integration error is completely negligible.

Our strategy in the numerical and seminumerical approaches is summarized schematically as

[TABLE]

where we have written the generic evolution kernel $U(\mu_{0},\mu)$ explicitly in terms of the perturbative series of the anomalous dimensions from eq. (3.3). In both the numerical and seminumerical methods we truncate the series in parentheses at the desired order in $\epsilon$, and evaluate the overall $\mu$-integration numerically.

For the numerical method, we use the exact numerical result for the solution of the running coupling $\alpha(\mu^{\prime})$, as indicated in red in eq. (3.8). For the seminumerical method we instead insert the iterative solution for the running of the coupling, eq. (2.9), as indicated in blue in eq. (3.8). In both cases, the running of the coupling is truncated at the appropriate order in $\epsilon$.

The seminumerical method turns out to provide a very good approximation of the integral, since as we have seen in section 2, the iterative solution for the running coupling provides a very accurate solution. At the same time, it is much faster than the fully numerical method because it avoids calling the computationally costly numerical solution for the running coupling in each integrand call.

3.3 Unexpanded analytic method

We start from eq. (3.2) and split the logarithm at $\mu_{0}$,

[TABLE]

To integrate it over $\mu$ we use the standard method of exploiting the $\beta$ function for $\alpha$,

[TABLE]

to change the integration variable from $\mu$ to $\alpha$. To replace the explicit logarithm of $\mu$ in eq. (3.9), we integrate eq. (3.10) once to obtain

[TABLE]

The evolution kernel in eq. (3.5) then takes the form

[TABLE]

where the individual functions are defined as

[TABLE]

The integrals in eqs. (3.13), (3.14), and (3.15) yield an evolution kernel which manifestly depends on $\mu$ and $\mu_{0}$ only via $\alpha(\mu)$ and $\alpha(\mu_{0})$. To illustrate this, for $K_{\Gamma}$ at LL order, $\sim\mathcal{O}(\epsilon^{0})$, we have

[TABLE]

The integrations can be carried out easily to yield

[TABLE]

At NLL, $\mathcal{O}(\epsilon)$, we have

[TABLE]

The approach followed in the unexpanded analytic method is to expand the denominators in $\epsilon$ keeping terms up to the $\mathcal{O}(\epsilon)$, which induces an approximation error of $\mathcal{O}(\epsilon^{2})$, but in turn allows one to easily perform the integration,

[TABLE]

The other two integrals, $\eta_{\Gamma}$ and $K_{\gamma}$, are obtained in a completely analogous fashion. At NLL, for $\eta_{\Gamma}$ we find

[TABLE]

We stress, that in this method the anomalous dimensions in $\eta_{\Gamma}$ are kept to the same loop order as in $K_{\Gamma}$, despite the fact that $\eta_{\Gamma}$ has an additional power of $\alpha(\mu_{0})$ compared to $K_{\Gamma}$ at each order in $\epsilon$. From the above derivation it should be clear that the separation of the cusp term into $K_{\Gamma}$ and $\eta_{\Gamma}$ is arbitrary and merely a technical tool to perform the integration. In other words, the fact that we used $\mu_{0}$ on the right-hand side of eq. (3.9) was merely for convenience and we could have used any other arbitrary scale. Dropping the highest term in $\eta_{\Gamma}$ (as is sometimes done) amounts to multiplying the $\ln(Q/\mu_{0})$ in eq. (3.9) by $\epsilon$, which introduces an artificial dependence on $\mu_{0}$ into $\gamma_{F}$. In particular, doing so would lead to an explicit violation of the RGE consistency in eq. (3.6).

From the expressions in eqs. (3.14) and (3.15) it is clear that $K_{\gamma}$ can be obtained from $\eta_{\Gamma}$ by replacing $\Gamma_{n}\to\gamma_{n}$ and multiplying by an overall $\epsilon$. Hence, it does not contribute at LL, and at NLL we have

[TABLE]

The calculation at higher orders proceeds in exactly the same fashion. The results up to N3LL are given in appendix A.1 in eqs. (A.1), (A.1), and (A.4).

To fully specify the method, we also have to specify how to obtain the integration limits $\alpha(\mu)$ and $\alpha(\mu_{0})$. To be consistent, we have to use a solution for the $\beta$ RGE to the same order in $\epsilon$ to which we performed the integration. To render the method fully analytic, we use the unexpanded iterative solution to the corresponding order.

The Sudakov evolution factor expressed in terms of $K_{\Gamma}$, $\eta_{\Gamma}$, $K_{\gamma}$ obtained as described above and combined with the iterative solution for $\alpha$ is precisely the form that is commonly used to perform Sudakov resummation in QCD in much of the SCET literature; a few examples from a variety of applications include e.g. refs. [39, 40, 34, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51].

3.4 Expanded analytic method

If, instead of changing integration variables from $\mu$ to $\alpha$ via eq. (3.10), one inserts the expanded solution eq. (2.1) for $\alpha(\mu)$ into the integrand of the evolution kernel eq. (3.5) and expands in $\epsilon$, the resulting integral exhibits an explicit $\mu$ dependence which can be integrated analytically.

Equivalently, the same result is obtained by starting from the unexpanded analytic integrals of the previous subsection [eqs. (A.1), (A.1), and (A.4)], substituting the expanded solution eq. (2.1) for $\alpha(\mu)$ in terms of $\alpha(\mu_{0})$ and expanding everywhere in $\epsilon$. We therefore refer to these kernels as “expanded analytic” as they involve a further expansion in $\epsilon$ compared to the unexpanded ones of the previous subsection.

The expanded analytic evolution kernels to NNLL, $\mathcal{O}(\epsilon^{2})$, are given by

[TABLE]

and

[TABLE]

where as before $b_{i}=\beta_{i}/\beta_{0}$, $\hat{\Gamma}_{i}=\Gamma_{i}/\Gamma_{0}$, and $X\equiv X(\mu_{0},\mu)=1+\frac{\alpha(\mu_{0})}{2\pi}\beta_{0}\ln\frac{\mu}{\mu_{0}}$. For illustration, we kept explicit the additional terms compared to the unexpanded results that correspond to the expansion of $r=\alpha(\mu)/\alpha(\mu_{0})$ in terms of $\epsilon$. The $K_{\gamma}$ kernel is again obtained from $\eta_{\Gamma}$ in eq. (3.4) by replacing $\gamma_{n}\rightarrow\Gamma_{n}$ and multiplying by an overall $\epsilon$. The N3LL results are obtained analogously. Since they are rather lengthy and not very illuminating we do not give them explicitly here.

As one can see, in this method the dependence on the two scales appears explicitly in the argument of the logarithm in $X(\mu_{0},\mu)$. One can still choose how to obtain the value for $\alpha(\mu_{0})$, for which we use by default the expanded solution as it is closer in spirit to this method. In our numerical results we will however also show the effect of using the iterative solution for $\alpha(\mu_{0})$.

3.5 Reexpanded analytic method

This method is related to the expanded analytic method of the previous subsection but uses a different treatment for the remaining dependence on $\alpha$. That is, starting from the expanded results of the previous subsection, all powers of $\alpha(\mu_{0})$ are reexpanded in terms of $\alpha(\mu_{R})$ evaluated at a reference scale $\mu_{R}$, which is typically chosen equal or proportional to the (hard) kinematic variable $Q$.444In our notation, taking Drell-Yan as an example, the kinematic variable $Q$ would be equivalent to the $Z$-boson or dilepton invariant mass. It should not be confused with what is sometimes called the resummation scale and also denoted as $Q\equiv Q_{\rm res}$, and which is the same as our $\mu_{0}$, i.e. $\mu_{0}\equiv Q_{\rm res}$.

To illustrate this for the simplest case, we start from the LL expanded analytic result, which written out explicitly is given by

[TABLE]

One now reexpands $\alpha(\mu_{0})$ in terms of $\alpha(\mu_{R})$ using its fixed-order expansion

[TABLE]

In doing so one assumes that $\mu_{0}\sim\mu_{R}$ or, more precisely, one explicitly chooses that the logarithms of $\mu_{0}/\mu_{R}$ are not resummed via the evolution of $\alpha$ but are instead treated at fixed order. Formally, this is implemented by multiplying any $\ln(\mu_{0}/\mu_{R})$ in the relation between $\alpha(\mu_{0})$ in terms of $\alpha(\mu_{R})$ by $\epsilon$ as in eq. (3.25) and expanding in $\epsilon$. Substituting eq. (3.25) into eq. (3.24) and reexpanding to $\mathcal{O}(\epsilon^{0})$, we obtain the “reexpanded analytic” result at LL,

[TABLE]

At this order, the result only involves the resummed logarithms $\ln(\mu_{0}/\mu)$. The results at higher orders also contain explicit fixed-order logarithms $\ln(\mu_{0}/\mu_{R})$ induced by the fixed-order expansion of $\alpha(\mu_{0})$ in terms of $\alpha(\mu_{R})$. Since the higher-order results are very lengthy we do not give them here. Their calculation is discussed for example in appendix C of ref. [52]. The LL result for $K_{\Gamma}$ in eq. (3.26) is equivalent to the $Lg^{(1)}(\alpha L)$ term in the notation there, and we also explicitly verified that we reproduce the NLL $g^{(2)}(\alpha L)$ term. The N3LL expressions are taken from ref. [53] translated to our conventions.

Another important difference in this method is the treatment of the $\eta_{\Gamma}$ term,

[TABLE]

in the Sudakov exponent. Since $\mu_{R}$ and $Q$ are either identified or considered of similar size, the explicit $\ln(Q/\mu_{0})$ here is treated analogously to the $\ln(\mu_{0}/\mu_{R})$ in eq. (3.25) and multiplied by $\epsilon$, such that the $\eta_{\Gamma}$ term is treated like the noncusp term and included at one order lower than $K_{\Gamma}$. This is typically achieved by absorbing it into the noncusp term as

[TABLE]

For the specific choice $\mu_{0}=\mu_{R}=Q$, this method reduces to the expanded analytic method, since all logarithms that are treated differently vanish exactly. However, as soon as these scales are chosen not to coincide the different treatment of the $\mu_{0}$ dependence can have a sizeable numerical effect on the evolution kernel, as we will see below. It is also important to note that since the first and second arguments of $U(\mu_{0},\mu)$ are explicitly treated differently, and $\mu_{0}$ is assumed to be of order $Q$, the group property in eq. (3.6) is lost, and similarly the closure condition in eq. (3.7) becomes meaningless. In other words, with the above modifications the evolution factor explicitly targets a specific situation and cannot (and is not meant to) be used to evolve between two arbitrary scales.

This reexpanded analytic expression for the Sudakov evolution factor is also commonly used, in particular in the formalism of refs. [54, 52, 55, 56] and (presumably most) implementations following it, and also in the formalism of refs. [57, 58, 53].

3.6 Numerical analysis of the evolution kernel

Having elaborated various methods for evaluating the Sudakov evolution factor, we now study their numerical behaviour. We consider the relative deviation from the exact closure condition in eq. (3.7) as well as the approximation error as given by the relative difference of each method to the exact numerical method.

In general, the Sudakov evolution factor depends on the process and observable under consideration, so we have to specify a concrete example. Here, we consider the hard evolution for two colored partons in QCD, namely the $q\bar{q}$ vector current (corresponding to Drell-Yan production or $e^{+}e^{-}\to$ dijets) and the $gg$ scalar current (corresponding to $gg\to$ Higgs production). The relevant hard evolution kernel is given by

[TABLE]

where $i=q,g$ denotes the quark and gluon cases, and we have separately defined the cusp $U^{i}_{\Gamma}$ and noncusp $U^{i}_{\gamma}$ evolution factors. All relevant anomalous dimension coefficients are given in appendix A.2.

In this case, the kinematic quantity $Q$ appearing in eq. (3.6) is the invariant mass $Q=\sqrt{q^{2}}$ of the momentum $q$ flowing through the current, which we take to be $Q=100\,\mathrm{GeV}$ (i.e. typical of Drell-Yan or Higgs production). We consider the case of evolving from $\mu_{0}\sim Q$ to an arbitrary scale $\mu$ over two decades above and below. As before, we use $n_{f}=5$ and ignore any flavor thresholds.

Of course, the takeaway message of our analysis would be the same if we were to use anomalous dimension coefficients associated with soft or collinear quantities. The advantage of considering the hard evolution is that it is multiplicative and independent of the low-energy observable, so it provides a simple and generic use case, while the only process dependence of the evolution is via the color channel. Furthermore, the hard evolution factor $U^{i}(\mu_{0},\mu)$ in this scenario has a direct correspondence in various resummation formalisms. In the formalisms of refs. [54, 52, 55, 56] and refs. [57, 58, 53] it corresponds to the Sudakov form factor or the Sudakov radiator with $\mu_{0}\equiv Q_{\rm res}$ being the resummation scale. In the context of SCET, it constitutes the evolution factor for the $q\bar{q}$ and $gg$ hard functions with $\mu_{0}\equiv\mu_{H}$ being the renormalization scale of the hard function. In all cases, $\mu$ would then be associated with an appropriate low-energy quantity, e.g. $\mu\sim b_{0}/b_{T}\sim p_{T}$ in the case of $p_{T}$ resummation.

3.6.1 Closure tests

By construction, the seminumerical and numerical methods satisfy the closure condition exactly, since they treat the integral and integration limits exactly. As already mentioned in subsection 3.5, the closure condition cannot be applied to the reexpanded method. We therefore only consider the unexpanded and expanded methods here. The reason these methods do not satisfy the closure condition exactly is due to the expansions in the integrand after the variable transformation from $\mu$ to $\alpha$, which means that the change of variables in the integration limits and the integrand are not in exact correspondence.555Note that in ref. [49] modified unexpanded kernels were constructed that restore exact closure.

In figure -895, we compare the non-closure for the full quark evolution factor at different orders for the unexpanded (left panel) and expanded (right panel) methods. At LL, the integrals are trivially exact and satisfy exact closure, so we do not show them. At NLL, the non-closure can be quite sizeable, exceeding $\gtrsim 5\%$ when running to low scales. For the expanded kernels it even reaches $\sim 1\%$ when running to high scales. While the non-closure effect is reduced at higher orders, it can still reach $1-2\%$ at the lowest scales, and for the unexpanded kernels it does not reduce from NNLL to N3LL. For the expanded kernels, the non-closure reduces to below $1\%$ at N3LL, but this is likely accidental, since the expanded kernels are very sensitive to numerical cancellations when evolving to scales $\mu\lesssim 10\,\mathrm{GeV}$. This is evident when comparing the effect of using the iterative instead of the expanded solution for $\alpha_{s}(\mu_{0})$: Even a small change in $\alpha_{s}(\mu_{0})$ causes large changes in the observed level of non-closure.

In figure -894, we compare the unexpanded and expanded methods to each other at NNLL (left panel) and N3LL (right panel). In addition we vary $\mu_{0}$ away from $Q$ to $Q/2$ and $2Q$, as one would do in practical applications to estimate a resummation uncertainty. Here, this should however not be considered as an uncertainty estimate on the non-closure. Rather, it illustrates the effect of $\eta_{\Gamma}$, which contributes when $\mu_{0}\neq Q$, and the level at which the non-closure may influence such uncertainty estimates. The (non-)closure of the unexpanded kernels tends to be less sensitive to the choice of $\mu_{0}$ than the expanded ones.

Overall, we might say that the unexpanded kernels show a somewhat better closure behaviour. However, considering that at N3LL one is aiming for perturbative precision in the several percent range, their non-closure at this order is uncomfortably large.

For brevity we have only shown results for the quark evolution kernels here. The gluon evolution kernels have the same qualitative behaviour. The only difference is that the overall non-closure effect is about a factor of two larger for gluons than for quarks, corresponding to their larger color factor.

3.6.2 Approximation errors

We now study the approximation errors of all four approximate methods relative to the exact numerical method. The results for the quark cusp contribution, $U^{q}_{\Gamma}$, are shown in figure -893, and for the quark and gluon noncusp contributions, $U^{q,g}_{\gamma}$, in figures -892 and -891. In all cases, we show the results at NNLL (top rows) and N3LL (bottom rows), and at $\mu_{0}=Q$ (left panels) and $\mu_{0}=Q/2$ (right panels).

At NNLL, we see a clear hierarchy, with the seminumerical method having the smallest approximation errors, followed by the unexpanded, and then the expanded methods, which is as expected given the increasing level of approximation involved in each method. At N3LL, the seminumerical method again performs best. The unexpanded and expanded kernels have similar errors for the cusp term, while the unexpanded ones fare better for the noncusp terms. The approximation error for the cusp term is always much larger than for the noncusp term, which is not surprising due to the additional $\ln(\mu)$ in its integrand. The errors for the gluon cusp contribution are about a factor of two larger than for quarks.

We find, somewhat surprisingly, that at NNLL the approximation error for the expanded kernels can exceed several percent when evolving to low scales, while at N3LL it still exceeds the percent level for both the unexpanded and expanded kernels alike. Overall, the picture is quantitatively quite similar to what we observed with the closure test.

We stress that the systematic differences we observe are solely due to the method of integrating the RGE for identical perturbative inputs. Typically, we would want such systematic effects to be much smaller than the perturbative precision we are aiming at, as was the case for the running of the coupling. In other words, we clearly want to avoid the method of integration to bias the result in any way. This is clearly not the case here, since at NNLL and N3LL one would typically aim at a perturbative precision of order several to few percent.

For $\mu_{0}=\mu_{R}=Q$, the reexpanded and expanded kernels are still equivalent. In the right panels, we therefore also show the results when reducing the hard scale of the evolution to $\mu_{0}=Q/2$. This has essentially no effect on the approximation error of the seminumerical, unexpanded, and expanded kernels. For the reexpanded kernels, we keep $\mu_{R}=Q$, which is a commonly used choice in resummation applications. The impact of treating the logarithms of $\mu_{0}/\mu_{R}$ and $Q/\mu_{0}$ at fixed-order in the reexpanded kernels compared to the expanded ones now becomes visible and turns out to be very large, which is quite unexpected. It easily exceeds the percent level even at N3LL, and not just for the cusp but even the noncusp contributions. For the cusp term, it exceeds $\gtrsim 10\%$ at the lowest scales. Note that roughly half of the observed difference is due to the reexpansion around $\alpha(\mu_{R})$ and half is due to the lower-order treatment of the $\eta_{\Gamma}$ term. This effect is not just due to the reduced amount of evolution from $Q/2$ to $\mu$, as that is present for all methods including the exact numerical result. Given the large numerical impact it has, it might be worthwhile to reconsider the reasons for performing this additional reexpansion around $\alpha(\mu_{R})$. Certainly, from the point of view of the evolution, the appropriate scale for $\alpha$ when starting the evolution from $\mu_{0}$ is $\mu_{0}$.

Of course, the differences due to the various approximations involved in all methods, including the reexpanded method, might be considered as higher-order effects. Nevertheless, given that they are not negligible or even exceed the perturbative precision they have to be addressed and accounted for. One option would be to include them as an additional systematic uncertainty in the theoretical uncertainty estimate. On the other hand, there is no fundamental theoretical reason for using any specific approximate solution. Hence, the best option would be to avoid incurring this additional uncertainty by using the unique exact solution of the defining RGE system resulting from the truncation of the perturbative series of the anomalous dimensions. If that is computationally prohibitive (or technically inconvenient), the seminumerical method offers a good compromise, since its approximation error is always well below the percent level, and so it is sufficiently accurate even for high-precision predictions.

4 Sudakov evolution kernels with two gauge interactions

Having investigated the integration of the one-dimensional kernels, we now consider the extension to two gauge interactions, in which case also mixed effects involving both gauge couplings $\alpha_{a,b}$ require resummation. In section 4.1, we review the general analytic structure for this case and the methods for evaluating them, based on what we learned from our exhaustive analysis of the one-dimensional case. We then present numerical results for the example of QCD$\otimes$QED in section 4.2.

4.1 Structure of the two-dimensional evolution kernel

We consider the Sudakov resummation for the direct product of two gauge groups $G_{a}\otimes G_{b}$. The extension to more groups is then straightforward. One key difference from the one-dimensional case is that the $\beta$ functions that govern the evolution of the couplings $\alpha_{a}$ and $\alpha_{b}$ now become a set of nonlinear coupled differential equations, as discussed in subsection 2.2.

The generic Sudakov RGE structure remains the same as in eq. (3.1), except that the perturbative expansions for all quantities now involve a double series in $\alpha_{a}$ and $\alpha_{b}$, including mixed terms corresponding to the emissions of two distinct gauge bosons. Hence, the all-order structure of the anomalous dimension is now given by666For the general EW case, the complete factorized structure of the cross section can become more involved, because the masses of the EW gauge bosons introduce an additional scale. The generic Sudakov RGE however still has the form of eq. (4.1) with at most a single logarithm $\ln(\mu)$ [6].

[TABLE]

with

[TABLE]

The bookkeeping parameters $\epsilon_{a,b}\equiv 1$ are the same as in the $\beta$ functions in eq. (2.2), and $\epsilon\sim\epsilon_{a}\sim\epsilon_{b}$, i.e. we make no assumption about the relative hierarchy of the two coupling constants. It is also important to note the correspondence to the typical notation for a single gauge theory $G_{a}$,777In contrast to the $\beta$ functions, $\alpha_{a}$ and $\alpha_{b}$ appear on equal footing in eqs. (4.2) and (4.3), so we have no choice but to increment the meaning of $n$ in the subscript compared to the one-dimensional case.

[TABLE]

The Sudakov evolution kernel is given by the two-dimensional analogue of eq. (3.5),

[TABLE]

It is clear from eqs. (4.2) and (4.3) that the sums implicit in eq. (4.5) include pure $G_{a}$ terms $\mathcal{O}(\epsilon_{a}^{n})$, pure $G_{b}$ terms $\mathcal{O}(\epsilon^{m}_{b})$ terms, and mixed terms $\mathcal{O}(\epsilon^{n}_{a}\epsilon^{m}_{b})$.

4.2 Evaluation of the two-dimensional evolution kernel

Evaluating eq. (4.5) does not correspond to a simple extension of the single gauge interaction scenario, since mixed terms $\sim\mathcal{O}(\alpha_{a}\alpha_{b})$ appear in conjunction with the coupled $\beta$ functions.

The fully numerical method is of course still applicable, although it is even more computationally demanding now, since multiple coupled differential equations for $\alpha_{a,b}$ must be solved. As before we use the fully numerical method to provide the exact reference result to which other methods are compared.

The unexpanded analytic method is not easily extendable, since the coupled $\beta$ functions do not allow an analogous change of variables along the lines of of $\mathrm{d}\ln\mu\rightarrow\mathrm{d}\alpha_{a}/\beta^{a}(\alpha_{a},\alpha_{b})$, because of the dependence on the second coupling. Doing so would require one to express the $\mu$ dependence of $\alpha_{b}$ in terms of $\alpha_{a}$, which in turn requires that one treats $\alpha_{b}$ as in the expanded analytic method, at which point the advantage of the unexpanded method is lost. In principle, this could still be an option for cases where there is a clear hierarchy between two couplings to justify treating them on unequal footing, as would be the case for QCD$\otimes$QED. However, we do not pursue this option further here for the reasons given below.

The expanded analytic method can still be applied to evaluate eq. (4.5). This is achieved by using the expanded solution of the coupled $\beta$ RGE for $\alpha_{a}$ and $\alpha_{b}$ obtained from eq. (2.2), substituting it into the perturbative expansions of the anomalous dimensions, and then explicitly performing the integration in terms of $\ln\mu$. This was done in ref. [35].

In either case, the obtained analytic kernels will inevitably suffer from the same-sized approximation errors already seen in the pure QCD case in subsection 3.6. The generalization to the product group $G_{a}\otimes G_{b}$ does not alter the ultimate source, which are the additional approximation(s) made in the integrand. By including EW or QED corrections, one is looking for percent-level effects, and based on our findings and discussion in the previous section, the analytic methods do not appear to provide sufficient numerical accuracy.888One might consider using a more accurate (semi)numerical method for the dominant QCD contributions, while including the smaller mixed and pure EW corrections via an analytic approximation. However, we do not see any gain in doing so compared to e.g. using the seminumerical method everywhere.

We therefore take the seminumerical method as our method of choice for evaluating the two-dimensional evolution kernels. As already seen in subsection 3.6, it features exact closure and very small approximation errors (well below the percent-level at NNLL), at a reasonable computational cost. The extension to the two-dimensional case eq. (4.5) only requires two steps:

We solve the coupled $\beta$ functions in eq. (2.2) via the iterative method, up to the required order in $\epsilon$, which yields the closed-form analytic expressions for $\alpha_{a}(\mu)$ and $\alpha_{b}(\mu)$ in eq. (2.2). 2. 2.

We then evaluate the evolution kernel $U(\mu,\mu_{0})$ by employing a numerical integration routine in eq. (4.5), using the analytic expressions for $\alpha_{a,b}(\mu)$ obtained in step 1 in the integrand.

Schematically, this procedure is illustrated by

[TABLE]

In what follows, we apply this method to the mixed gauge QCD$\otimes$QED scenario. We consider the $u\bar{u}$ hard function as a concrete example. The relevant coefficients up to three loops can be found in appendices A.2 and A.3, allowing us to obtain the complete NNLL joint QCD$\otimes$QED Sudakov evolution.

In figure -890, we show the approximation errors for the seminumerical method at NLL and NNLL for the cusp (left panel) and noncusp (right panel) contributions. As for the pure QCD case, the former has a larger approximation error, while overall it remains well below 1% everywhere (except at NLL for the cusp term at the very lowest $\mu$ values).

In figure -889, we show the relative impact of the QED corrections by comparing the full QCD$\otimes$QED to the pure QCD resummation kernels at each order. For the cusp piece (left panel) the impact reaches $4-5\%$ over two decades of evolution, while for the noncusp piece (right panel) it reaches close to $1\%$. The overall effect for QED is expectedly small due to the smallness of the electromagnetic coupling. Considering a toy QCD$\otimes$QED*′* model with $\alpha_{e}^{\prime}(m_{Z})=0.1$ in figure -888, the corrections become much larger, $30-50\%$ for the cusp piece and $5-15\%$ for the noncusp piece. Note that the impact is driven not only by pure QED and mixed corrections in the expansion of the anomalous dimensions, but also by the effects induced by the mixed $\beta$ functions in the coupling evolution. Note also that the impact for the cusp term is practically identical at NLL and NNLL even for QED*′*. This is somewhat accidental and due to the fact that the cusp anomalous dimension does not yet receive any mixed contributions at three loops while its pure QCD and QED three-loop coefficients happen to be very small.

5 Conclusions

High-precision experimental measurements of (e.g.) $W$ and $Z$ production at the LHC at the (sub-)percent level, as well as future measurements at the high-luminosity LHC or future high-energy colliders, demand equally precise resummed predictions in QCD$\otimes$EW. In this paper we studied the technical aspects of achieving this joint resummation at higher orders in generic coupled gauge environments.

In particular, at the level of the Sudakov evolution kernel appearing in any resummation formalism, we showed that the commonly used methods of evaluating the associated integrals in the Sudakov exponent via analytic approximations can cause numerical differences of the same size as the contributions coming from moving to a higher logarithmic accuracy or including EW radiation. In other words, systematic integration errors that are typically assumed to be subleading cannot be overlooked when attempting percent-level precision.

To show this we first studied five methods for integrating the evolution kernels in the case of a single gauge interaction: numerical, seminumerical, unexpanded analytic, expanded analytic, and reexpanded analytic. Their main difference is in the treatment of the $\mu$ dependence of $\alpha(\mu)$ that one has to integrate over. Although all methods employ a priori justifiable assumptions, we showed that the latter three analytic methods can introduce errors at and above the percent-level compared to the exact result (obtained via a fully numerical treatment at a sufficiently high numerical precision). One should therefore be cautious in using them, since their approximation errors can be nonnegligible compared to the perturbative precision one is aiming for. The reexpanded kernels, which are often used in the literature, are particularly notable; differences greater than $10\%$ are possible over only two decades of evolution.

On the other hand, we found that a seminumerical method, which combines an accurate analytic approximation for the running of the couplings with a numerical integration, yielded minimal errors (typically $\leq 0.1\%$), while still maintaining reasonable computing times. We therefore advocate its use for phenomenological studies at high precision and in joint resummation environments. As an example, we have used it to obtain the complete NNLL QCD$\otimes$QED Sudakov evolution factor.

Acknowledgments

We would like to thank Iain Stewart, Peter Marquard, and Johannes Michel for discussions.

Appendix A Perturbative results

In this appendix we summarize the results for the higher-order analytic RGE solutions, and the coefficients of the QCD and QED $\beta$ functions and anomalous dimensions.

A.1 RGE solutions

The iterative solution of the $\beta$-function RGE up to N3LO, $\mathcal{O}(\epsilon^{3})$, reads

[TABLE]

where $b_{n}=\beta_{n}/\beta_{0}$ and $X=1+\frac{\alpha(\mu_{0})}{2\pi}\,\beta_{0}\,\ln(\mu/\mu_{0})$ as in eq. (2.3).

The functions $K_{\Gamma}$, $\eta_{\Gamma}$, and $K_{\gamma}$ appearing in the unexpanded analytic Sudakov exponents in eq. (3.12) are given up to N3LL, $\mathcal{O}(\epsilon^{3})$, by

[TABLE]

where $b_{n}=\beta_{n}/\beta_{0}$, $\hat{\Gamma}_{n}=\Gamma_{n}/\Gamma_{0}$, $\hat{\gamma}_{n}=\gamma_{n}/\gamma_{0}$, and $r=\alpha_{s}(\mu)/\alpha_{s}(\mu_{0})$.

The corresponding kernels for the expanded analytic method are obtained by inserting eq. (A.1) into the above and expanding the results to $\mathcal{O}(\epsilon^{3})$.

A.2 QCD anomalous dimensions

Here we collect the pure QCD anomalous dimensions. For clarity and to avoid any confusion we use the two-dimensional notation everywhere. The QCD $\beta$-function coefficients in the $\overline{\text{MS}}$ scheme up to 3 loops are [59, 60]

[TABLE]

The four-loop coefficient for $N_{c}=3$ is given by [61, 62]

[TABLE]

The $\overline{\text{MS}}$ quark cusp anomalous dimension coefficients up to three loops are [63, 64]

[TABLE]

The four-loop coefficient entering at N3LL has become available during recent years, see e.g. refs. [65, 66, 67, 68, 69]. We use the numerical result for $N_{c}=3$ and $n_{f}=5$ obtained from ref. [69],

[TABLE]

Finally, the $\overline{\text{MS}}$ noncusp quark anomalous dimension coefficients for the hard function up to three loops are [70, 71, 72]

[TABLE]

The corresponding gluon noncusp anomalous dimensions can be found e.g. in appendix A.2 of ref. [45].

A.3 QED and mixed QCD$\otimes$QED anomalous dimensions

Here we collect all QED and QCD$\otimes$QED coefficients that are required for the full NNLL hard evolution, and which enter in the numerical results in subsection 2.2 and section 4. The mixed QCD $\beta$-function coefficients to three loops are

[TABLE]

while the pure and mixed QED $\beta$-function coefficients up to three loops are

[TABLE]

Here we defined

[TABLE]

where the sum runs over the active quark flavors with $Q_{u,c}=2/3$ and $Q_{d,s,b}=-1/3$ the quark charges, $N_{c}=3$ the number of colors, $Q_{\ell}=-1$ the lepton charge, and $n_{\ell}=3$ the number of charged leptons. The extraction of the three-loop mixed coefficients is discussed in appendix B.

The QED coefficients for the cusp anomalous dimension are obtained straightforwardly by taking the abelian QED limit of eq. (A.2),

[TABLE]

Up to three loops, the mixed coefficients vanish, as was noted in ref. [7],

[TABLE]

A nonzero mixed contribution is expected to first appear at four loops.

The noncusp anomalous dimensions are

[TABLE]

The pure QED coefficients are given by the abelian QED limit of the QCD coefficients in eq. (A.2). The result for the mixed two-loop coefficient, $\gamma^{q}_{H(1,1)}$, is obtained by following the abelianization procedure in ref. [32], and agrees with that given in ref. [36].

Appendix B Extraction of three-loop mixed QCD$\otimes$QED coefficients

In this section, we discuss our extraction of the three-loop coefficients of the mixed QCD$\otimes$QED $\beta$ functions from the results in ref. [37], which are required for the complete NNLO running.

Ref. [37] considers the coupled $\beta$-function RGE for a generic gauge group given as the product of $n$ simple groups, $G=G_{1}\otimes G_{2}\otimes\dots\otimes G_{n}$. They explicitly calculate the case of three distinct simple groups $G=G_{1}\otimes G_{2}\otimes G_{3}$, since at most three gauge bosons can propagate simultaneously at three loops. They consider the case where each simple group $G_{i}$ is nonabelian, but it is straightforward to apply their results to the abelian case by a proper modification of the Casimir invariants.

To obtain the QCD$\otimes$QED coefficients, we specify ourselves to the product group $G=SU(3)_{c}\otimes U(1)_{EM}\otimes\mathbf{1}$, where $\mathbf{1}$ is the trivial identity group. In addition, we only keep fermionic matter couplings, setting the Yukawa and scalar (quartic) couplings and the Casimir invariants of the scalar representation $S$ to zero, as they appear in eq. (3.1) of ref. [37]. We have explicitly checked that with this procedure we reproduce the pure QCD $\beta$-function coefficients.

Since we only have two couplings, namely $\alpha_{i}$ with $i=s,e$ for QCD and QED respectively, the usual Casimir invariants for the $SU(3)_{c}$ gauge group are

[TABLE]

where $F_{s}$ and $G_{s}$ stand for the fundamental and adjoint representations, so $C(F_{s})\equiv C_{F}$ and $C(G_{s})\equiv C_{A}$. For the case of $U(1)_{EM}$, we take the abelian limit

[TABLE]

with $f$ representing quarks and leptons. Note that sums over fermion species do appear in the $\beta$-function coefficients that each gauge group involves, i.e. quarks for $SU(3)_{c}$ and both quarks and leptons for $U(1)_{EM}$. Also, since ref. [37] decomposes Dirac fermions into chiral fermions, we have to substitute $n_{f}\rightarrow 2n_{f}$ for the number of quarks and $n_{\ell}\rightarrow 2n_{\ell}$ for the number of leptons in their results.

All of this considered, the $\beta^{e}_{11}$ coefficient can be extracted from the generic results of ref. [37], finding999We found a typo in this term in eq. (3.1) of ref. [37], which is missing the last $D(F_{i})$ factor. This is confirmed by comparing this term for the full SM gauge group with the results of refs. [73, 74].

[TABLE]

where $i=e$, since we are considering the $U(1)_{EM}$ $\beta$ function, and $j=s$. These specifications then give

[TABLE]

where the fermion sum over the leptons does not contribute due to the presence of $C(\ell_{s})=0$ and the multiplicity factors

[TABLE]

in our case correspond to the dimension of the trivial group factor and of the fundamental $SU(3)_{c}$ representation.

We can also read off the $\beta^{e}_{02}$ coefficient101010We also found a typo in these terms in eq. (3.1) of ref. [37], which are missing the last $D(F_{i})$ and $D(F_{m,i})$ factors. This is confirmed by comparing the result for $\beta^{e}_{02}$ with the results in ref. [75].

[TABLE]

where again $i=e$ and $j=s$. We therefore find

[TABLE]

Observe that again only the sum over quarks contributes, since the leptonic sums include a $C(\ell_{s})=T(\ell_{s})=0$.

Following the same procedure we also obtained the mixed QCD coefficients $\beta^{s}_{11}$ and $\beta^{s}_{02}$ as given in eq. (A.3), in agreement with the corresponding results given in ref. [76], where $\beta_{11}^{s}$ was obtained from an explicit three-loop calculation.