Higher-Order Sudakov Resummation in Coupled Gauge Theories

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.


Introduction
In perturbative quantum field theories it is well known that observables sensitive to physics at different scales µ µ 0 can have their perturbative expansions in the coupling constant α generically enhanced by Sudakov logarithms of the form α(µ) n ln m µ µ 0 with m ≤ 2n . (1.1) 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 α, 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 N 3 LL 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. α e ∼ α 2 s ). 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⊗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 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⊗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 O(α e α s ) mixed QCD⊗QED corrections and O(α 2 e ) 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⊗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 ⊗ 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 An evolution factor of this form necessarily appears in all formulations of higher-order Sudakov resummation. It implicitly depends on the β functions controlling the renormalization group evolution (RGE) of both gauge coupling constants α a (µ) and α b (µ), which in general are a coupled system of differential equations, as well as anomalous dimensions (Γ cusp , γ) whose perturbative expansions are themselves joint series in α 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⊗QED, i.e. G a ≡ SU (3) c and G b ≡ U (1) 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 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 µ-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 β functions for α a,b . We first review the standard approximate analytic solutions for α 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 N 3 LL. 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⊗QED coefficients. In appendix B, we discuss the extraction of the complete three-loop mixed QCD⊗QED β-function coefficients.

Iterative solutions to β-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 β-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.

Single gauge theory
We start from the well-known β-function RGE for the coupling constant α(µ) of a generic gauge theory, Here, µ is the renormalization scale, and we introduced a formal expansion parameter ≡ 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 We take the viewpoint that eq. (2.1) defines the running order of the RGE. That is, keeping the terms up to O( k ) in eq. (2.1) defines the N k LO or (k + 1)-loop running of α(µ). At leading order, O( 0 ), eq. (2.1) has the well-known exact analytic solution, where α(µ 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. 1 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 µ is not fixed but dynamical. This is precisely the case for the Sudakov evolution kernel, for which we will need to integrate α(µ) over µ.
In what follows, we review an iterative method to obtain an approximate analytic solution for eq. (2.1). At NLO, O( ), the β-function RGE reads In the second line we substituted the LO solution in eq. (2.3) for the µ dependence of α(µ) in the O( ) term. This induces an O( 2 ) error, since the difference between the LO and NLO µ dependence will itself be of O( ). Since the µ dependence in the term in square brackets is now explicit, it can be easily integrated, yielding the NLO solution We refer to this (and its higher-order analogues) as the "iterative" solution for α(µ). We can also expand the inverse in eq. (2.5) in to obtain 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 β-function RGE. For the iterative solution in eq. (2.5), it is trivial to verify that upon differentiation with respect to ln µ it reproduces the RGE as given in the second line of eq. (2.4). On the other hand, taking the ln µ derivative of the expanded NLO solution in eq. (2.6) yields Comparing this with eqs. (2.4) and (2.6), we see that the term in square brackets corresponds to expanding the overall α(µ) 2 in the β function to O( ). 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, O( 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 O( ) and O( 2 ) terms, keeping all terms up to O( 2 ) as well as the overall α(µ) 2 , The µ-dependence in the square brackets on the second line, encoded in X ≡ X(µ 0 , µ), is again simple enough to be integrated analytically. 2 This yields the NNLO iterative solution 2 This would not be the case had we used the unexpanded NLO solution from eq. (2.5) in the O( ) term. Expanding the inverse in eq. (2.9) in , we obtain the expanded NNLO solution, It is straightforward to extend the iterative solution to higher orders. To obtain the N k LO solution, one simply inserts the N k−1 LO solution for α(µ) into the N k LO β function and expands it to O( k ), while keeping the overall α(µ) 2 exact. The N 3 LO solution is given in appendix A.1.
To illustrate the numerical precision of the approximate analytic solutions, we take the QCD coupling constant α s (µ) as an example, using α 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 exact 3 numerical solution, which we refer to as the approximation error, is shown in figure 1 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 µ = 2 GeV, and at most 0.3% at µ = 1 GeV. For running above m Z , the approximation error is much smaller due to asymptotic freedom. Beyond the highest scale shown, µ = 10 4 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.

Two coupled gauge theories
We now consider the case of two gauge theories. Their β 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⊗QED running. For now, we keep the discussion general and consider the following set of coupled β-function RGEs for two gauge couplings α a (µ) and α b (µ), (2.11) We have again introduced formal expansion parameters a,b ≡ 1 to easily keep track of the evolution order. For future convenience, we also define the rescaled coefficients (2.12) Note that by definition β x n0 ≡ β n and b x n0 ≡ b n are the coefficients of the individual gauge theories in the absence of the second.
As before, the order in a,b to which eq. (2.11) is expanded is what defines the running order of the couplings. Generically, we consider a,b on equal footing and define the N n LO evolution by including in eq. (2.11) all terms of O( k a n−k b ) with 0 ≤ k ≤ 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 a and b that are included in eq. (2.11).
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 β-function RGE system by extending the iterative method in subsection 2.1. The key property of eq. (2.11) that allows us to do so is that, at leading order, O( 0 a , 0 b ), the RGE system decouples, yielding exact LO solutions: To obtain an approximate NLO solution, it then suffices to substitute the above LO solutions into the O( a ) and O( b ) terms of the NLO RGE system, which induces O( 2 a , a b , 2 b ) errors, (2.14) The terms in square brackets can now be explicitly integrated over ln µ to obtain the iterative NLO solution The solution for α b (µ) is given by replacing a ↔ b everywhere. Note that in eq. (2.15) the utility of the a,b counting parameters becomes evident. Naively, one might have expected that the O( a ) and O( b ) terms will be proportional to α a (µ 0 ) and α b (µ 0 ) respectively, which however is not the case, as both are proportional to α a (µ 0 ). Instead, the a,b actually keep track of the fact that the ln X a,b factors respectively resum a series of α n a,b ln n (µ/µ 0 ) terms. The iterative NNLO solution is obtained by substituting the NLO solutions into the β- (2.16) As before, the solution for α b (µ) is obtained by replacing a ↔ b. The terms in the first line correspond to the NNLO solution of α a in the absence of α b , while the remaining ones are the mixing contributions involving at least one power of b . The corresponding expanded solution is obtained by inverting eq. (2.16) and expanding it in i . Note that when doing so, it becomes essential to expand in terms of i and not α i (µ 0 ). The necessity of using i as expansion parameter is also evident in the mixed O( a 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⊗QED. The relevant coefficients for the coupled β-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 β e 11 . 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 α s (m Z ) = 0.118 and α e (m Z ) = 1/127 as boundary conditions, n f = 5 for the number of active quark flavors, and n = 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 2 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 α 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 β function. In figures 3 and 4 we show them for the case of QCD⊗QED as well as for a toy QCD⊗QED for which we set the boundary condition to α e (m Z ) = 1/10. In figure 3, as expected, the QED corrections ∼ O( 0 s , 1 e ) have almost negligible effect on α s , with only the higher order QCD corrections ∼ O( 1 s , 0 e ) 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⊗QED scenario shown in figure 4. 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 β-function coefficients β s 10 and β 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 β-function coefficients are of similar numerical size.

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.

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 where γ F (µ) is the anomalous dimension, and the function F (µ) is any one of the factorized ingredients. The Sudakov resummation is then performed by solving the resulting coupled system of RGE equations.
In general, γ F (µ) and F (µ) both depend on an additional external kinematic quantity (which appears as part of the argument of the Sudakov logarithms in the cross section). The ⊗ denotes the fact that γ 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 where Q denotes the above-mentioned kinematic quantity, Γ cusp is (proportional to) the cusp anomalous dimension, and γ is the noncusp anomalous dimension. They obey the perturbative expansions We have again introduced the formal expansion parameter ≡ 1, which we will use to define the resummation order. Since the µ dependence of γ F (µ) primarily enters via the coupling constant's µ dependence, the β RGE for α(µ) in eq. (2.1) is an integral part of the full RGE system to be solved. In particular, the 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 fundamentally defines the resummation order. That is, keeping terms up to O( k ) defines the Sudakov evolution at N k LL order. The explicit 1/ 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 β function in eq. (2.1) starts at O( 0 ), it enters at the same loop order as the cusp anomalous dimension. So as usual, at N k LL 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 where U (µ 0 , µ) is the Sudakov evolution factor given by It resums the Sudakov logarithms appearing in the perturbative series of F (µ)/F (µ 0 ). One can easily check that counting powers of in the anomalous dimension is equivalent to counting powers of logarithms in the Sudakov exponent, simply because the entire structure of U (µ 0 , µ) is fully encoded by the anomalous dimension. In the final resummed cross section, where F (µ) is combined with other ingredients, U (µ 0 , µ) eventually appears with both scales µ and µ 0 corresponding to two different kinematic quantities such that U (µ 0 , µ) 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 β-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 β-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 α(µ) 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 β function to turn the µ integration into an integration over α, which can be performed analytically after some expansion in . This is combined with the analytic iterative solution for α s (µ) 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 α(µ) in eq. (2.10) in the perturbative expansions of Γ(α) and γ(α) and fully expanding the integrand in . The integration can then be performed directly in terms of µ. 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 α(µ 0 ) in terms of α(µ R ) at a different reference scale µ R , assuming that µ 0 ∼ µ 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 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 µ 2 = µ 0 , which yields 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.

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 where we have written the generic evolution kernel U (µ 0 , µ) 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 , and evaluate the overall µ-integration numerically.
For the numerical method, we use the exact numerical result for the solution of the running coupling α(µ ), 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 .
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.

Unexpanded analytic method
We start from eq. (3.2) and split the logarithm at µ 0 , To integrate it over µ we use the standard method of exploiting the β function for α, to change the integration variable from µ to α. To replace the explicit logarithm of µ in eq. (3.9), we integrate eq. (3.10) once to obtain The evolution kernel in eq. (3.5) then takes the form where the individual functions are defined as 14) The integrals in eqs. (3.13), (3.14), and (3.15) yield an evolution kernel which manifestly depends on µ and µ 0 only via α(µ) and α(µ 0 ). To illustrate this, for K Γ at LL order, The integrations can be carried out easily to yield The approach followed in the unexpanded analytic method is to expand the denominators in keeping terms up to the O( ), which induces an approximation error of O( 2 ), but in turn allows one to easily perform the integration, The other two integrals, η Γ and K γ , are obtained in a completely analogous fashion. At NLL, for η Γ we find We stress, that in this method the anomalous dimensions in η Γ are kept to the same loop order as in K Γ , despite the fact that η Γ has an additional power of α(µ 0 ) compared to K Γ at each order in . From the above derivation it should be clear that the separation of the cusp term into K Γ and η Γ is arbitrary and merely a technical tool to perform the integration. In other words, the fact that we used µ 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 η Γ (as is sometimes done) amounts to multiplying the ln(Q/µ 0 ) in eq. (3.9) by , which introduces an artificial dependence on µ 0 into γ 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 γ can be obtained from η Γ by replacing Γ n → γ n and multiplying by an overall . Hence, it does not contribute at LL, and at NLL we have To fully specify the method, we also have to specify how to obtain the integration limits α(µ) and α(µ 0 ). To be consistent, we have to use a solution for the β RGE to the same order in to which we performed the integration. To render the method fully analytic, we use the unexpanded iterative solution to the corresponding order.

Expanded analytic method
If, instead of changing integration variables from µ to α via eq. (3.10), one inserts the expanded solution eq. (2.10) for α(µ) into the integrand of the evolution kernel eq. (3.5) and expands in , the resulting integral exhibits an explicit µ 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.2), (A.3), and (A.4)], substituting the expanded solution eq. (2.10) for α(µ) in terms of α(µ 0 ) and expanding everywhere in . We therefore refer to these kernels as "expanded analytic" as they involve a further expansion in compared to the unexpanded ones of the previous subsection.
The expanded analytic evolution kernels to NNLL, O( 2 ), are given by and where as before b i = β i /β 0 ,Γ i = Γ i /Γ 0 , and X ≡ X(µ 0 , µ) = 1 + α(µ 0 ) 2π β 0 ln µ µ 0 . For illustration, we kept explicit the additional terms compared to the unexpanded results that correspond to the expansion of r = α(µ)/α(µ 0 ) in terms of . The K γ kernel is again obtained from η Γ in eq. (3.23) by replacing γ n → Γ n and multiplying by an overall . The N 3 LL 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(µ 0 , µ). One can still choose how to obtain the value for α(µ 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 α(µ 0 ).

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 α. That is, starting from the expanded results of the previous subsection, all powers of α(µ 0 ) are reexpanded in terms of α(µ R ) evaluated at a reference scale µ R , which is typically chosen equal or proportional to the (hard) kinematic variable Q. 4 To illustrate this for the simplest case, we start from the LL expanded analytic result, which written out explicitly is given by One now reexpands α(µ 0 ) in terms of α(µ R ) using its fixed-order expansion In doing so one assumes that µ 0 ∼ µ R or, more precisely, one explicitly chooses that the logarithms of µ 0 /µ R are not resummed via the evolution of α but are instead treated at fixed order. Formally, this is implemented by multiplying any ln(µ 0 /µ R ) in the relation between α(µ 0 ) in terms of α(µ R ) by as in eq. (3.25) and expanding in . Substituting eq. (3.25) into eq. (3.24) and reexpanding to O( 0 ), we obtain the "reexpanded analytic" result at LL, At this order, the result only involves the resummed logarithms ln(µ 0 /µ). The results at higher orders also contain explicit fixed-order logarithms ln(µ 0 /µ R ) induced by the fixedorder expansion of α(µ 0 ) in terms of α(µ 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 Γ in eq. (3.26) is equivalent to the Lg (1) (αL) term in the notation there, and we also explicitly verified that we reproduce the NLL g (2) (αL) term. The N 3 LL expressions are taken from ref. [53] translated to our conventions.
Another important difference in this method is the treatment of the η Γ term, in the Sudakov exponent. Since µ R and Q are either identified or considered of similar size, the explicit ln(Q/µ 0 ) here is treated analogously to the ln(µ 0 /µ R ) in eq. (3.25) and multiplied by , such that the η Γ term is treated like the noncusp term and included at one order lower than K Γ . This is typically achieved by absorbing it into the noncusp term as For the specific choice µ 0 = µ 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 µ 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 (µ 0 , µ) are explicitly treated differently, and µ 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. [52,[54][55][56] and (presumably most) implementations following it, and also in the formalism of refs. [53,57,58].

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 qq vector current (corresponding to Drell-Yan production or e + e − → dijets) and the gg scalar current (corresponding to gg → Higgs production). The relevant hard evolution kernel is given by (3.29) where i = q, g denotes the quark and gluon cases, and we have separately defined the cusp U i Γ and noncusp U i γ evolution factors. All relevant anomalous dimension coefficients are given in appendix A.2.
In this case, the kinematic quantity Q appearing in eq. (3.29) is the invariant mass Q = q 2 of the momentum q flowing through the current, which we take to be Q = 100 GeV (i.e. typical of Drell-Yan or Higgs production). We consider the case of evolving from µ 0 ∼ Q to an arbitrary scale µ 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  Figure 5. Deviation from the closure condition U q (µ 0 , µ)U q (µ, µ 0 ) = 1 for the quark evolution kernels from NLL to N 3 LL for the unexpanded method (left) and expanded method (right). For the expanded method, the gray lines show the result of using the iterative instead of the expanded α s (µ 0 ).  Figure 6. Deviation from the closure condition U q (µ 0 , µ)U q (µ, µ 0 ) = 1 at NNLL (left) and N 3 LL (right) for the unexpanded method (solid) and expanded method (dashed). The central lines are for µ 0 = Q, while the bands show the variation when changing µ 0 to Q/2 and 2Q. The green dotted line shows the result of using the iterative α s (µ 0 ) in the expanded kernels (for µ 0 = Q only). 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 (µ 0 , µ) in this scenario has a direct correspondence in various resummation formalisms. In the formalisms of refs. [52,[54][55][56] and refs. [53,57,58] it corresponds to the Sudakov form factor or the Sudakov radiator with µ 0 ≡ Q res being the resummation scale. In the context of SCET, it constitutes the evolution factor for the qq and gg hard functions with µ 0 ≡ µ H being the renormalization scale of the hard function. In all cases, µ would then be associated with an appropriate low-energy quantity, e.g. µ ∼ b 0 /b T ∼ p T in the case of p T resummation.

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 µ to α, which means that the change of variables in the integration limits and the integrand are not in exact correspondence. 5 In figure 5, 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 5% when running to low scales. For the expanded kernels it even reaches ∼ 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 N 3 LL. For the expanded kernels, the non-closure reduces to below 1% at N 3 LL, but this is likely accidental, since the expanded kernels are very sensitive to numerical cancellations when evolving to scales µ 10 GeV. This is evident when comparing the effect of using the iterative instead of the expanded solution for α s (µ 0 ): Even a small change in α s (µ 0 ) causes large changes in the observed level of non-closure.
In figure 6, we compare the unexpanded and expanded methods to each other at NNLL (left panel) and N 3 LL (right panel). In addition we vary µ 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 η Γ , which contributes when µ 0 = Q, and the level at which the nonclosure may influence such uncertainty estimates. The (non-)closure of the unexpanded kernels tends to be less sensitive to the choice of µ 0 than the expanded ones.
Overall, we might say that the unexpanded kernels show a somewhat better closure behaviour. However, considering that at N 3 LL 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.

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 Γ , are shown in figure 7, and for the quark and gluon noncusp contributions, U q,g γ , in figures 8 and 9. In all cases, we show the results at NNLL (top rows) and N 3 LL (bottom rows), and at µ 0 = Q (left panels) and µ 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 N 3 LL, 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(µ) 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 N 3 LL 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 N 3 LL one would typically aim at a perturbative precision of order several to few percent.
For µ 0 = µ 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 µ 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 µ R = Q, which is a commonly used choice in resummation applications. The impact of treating the logarithms of µ 0 /µ R and Q/µ 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 N 3 LL, and not just for the cusp but even the noncusp contributions. For the cusp term, it exceeds 10% at the lowest scales. Note that roughly half of the observed difference is due to the reexpansion around α(µ R ) and half is due to the lower-order treatment of the η Γ term. This effect is not just due to the reduced amount of evolution from Q/2 to µ, 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 α(µ R ). Certainly, from the point of view of the evolution, the appropriate scale for α when starting the evolution from µ 0 is µ 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.

Structure of the two-dimensional evolution kernel
We consider the Sudakov resummation for the direct product of two gauge groups G a ⊗ G b . The extension to more groups is then straightforward. One key difference from the onedimensional case is that the β functions that govern the evolution of the couplings α a and α 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 α a and α 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 by 6 The bookkeeping parameters a,b ≡ 1 are the same as in the β functions in eq. (2.11), and ∼ a ∼ 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 , 7 The Sudakov evolution kernel is given by the two-dimensional analogue of eq. (3.5),

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 ∼ O(α a α b ) appear in conjunction with the coupled β functions.
The fully numerical method is of course still applicable, although it is even more computationally demanding now, since multiple coupled differential equations for α 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 β functions do not allow an analogous change of variables along the lines of of d ln µ → dα a /β a (α a , α b ), because of the dependence on the second coupling. Doing so would require one to express the µ dependence of α b in terms of α a , which in turn requires that one treats α 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⊗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 β RGE for α a and α b obtained from eq. (2.16), substituting it into the perturbative expansions of the anomalous dimensions, and then explicitly performing the integration in terms of ln µ. 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 ⊗ 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. 8 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: 1. We solve the coupled β functions in eq. (2.11) via the iterative method, up to the required order in , which yields the closed-form analytic expressions for α a (µ) and α b (µ) in eq. (2.16).
2. We then evaluate the evolution kernel U (µ, µ 0 ) by employing a numerical integration routine in eq. (4.5), using the analytic expressions for α a,b (µ) obtained in step 1 in the integrand. .

(4.6)
In what follows, we apply this method to the mixed gauge QCD⊗QED scenario. We consider the 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⊗QED Sudakov evolution.
In figure 10, 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 µ values).
In figure 11, we show the relative impact of the QED corrections by comparing the full QCD⊗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⊗QED model with α e (m Z ) = 0.1 in figure 12, 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 β 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  Figure 11. Impact of the QED corrections on the u-quark Sudakov evolution kernel at LL, NLL, and NNLL for the cusp (left panel) and noncusp (right panel) contributions. We show the relative difference of the full QCD⊗QED evolution to the pure QCD case at the corresponding order. receive any mixed contributions at three loops while its pure QCD and QED three-loop coefficients happen to be very small.

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⊗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 µ dependence of α(µ) 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 ≤ 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⊗QED Sudakov evolution factor.

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 β-function coefficients in the MS scheme up to 3 loops are [59,60] β s 00 ≡ β 0 = The four-loop coefficient for N c = 3 is given by [61,62]  The pure QED coefficients are given by the abelian QED limit of the QCD coefficients in eq. (A.9). The result for the mixed two-loop coefficient, γ q H(1,1) , is obtained by following the abelianization procedure in ref. [32], and agrees with that given in ref. [36].

B Extraction of three-loop mixed QCD⊗QED coefficients
In this section, we discuss our extraction of the three-loop coefficients of the mixed QCD⊗QED β functions from the results in ref. [37], which are required for the complete NNLO running.
Ref. [37] considers the coupled β-function RGE for a generic gauge group given as the product of n simple groups, G = G 1 ⊗ G 2 ⊗ · · · ⊗ G n . They explicitly calculate the case of three distinct simple groups G = G 1 ⊗ G 2 ⊗ 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⊗QED coefficients, we specify ourselves to the product group G = SU (3) c ⊗ U (1) EM ⊗ 1, where 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 β-function coefficients.
Since we only have two couplings, namely α i with i = s, e for QCD and QED respectively, the usual Casimir invariants for the SU (3) c gauge group are where F s and G s stand for the fundamental and adjoint representations, so C(F s ) ≡ C F and C(G s ) ≡ C A . For the case of U (1) EM , we take the abelian limit with f representing quarks and leptons. Note that sums over fermion species do appear in the β-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 → 2n f for the number of quarks and n → 2n for the number of leptons in their results. All of this considered, the β e 11 coefficient can be extracted from the generic results of ref. [37], finding 9 β e 11 = Observe that again only the sum over quarks contributes, since the leptonic sums include a C( s ) = T ( s ) = 0. Following the same procedure we also obtained the mixed QCD coefficients β s 11 and β s 02 as given in eq. (A.10), in agreement with the corresponding results given in ref. [76], where β s 11 was obtained from an explicit three-loop calculation.