Light-cone distribution amplitudes of light mesons with QED effects

We discuss the generalization of the leading-twist light-cone distribution amplitude for light mesons including QED effects. This generalization was introduced to describe virtual collinear photon exchanges above the strong-interaction scale $\Lambda_{\rm QCD}$ in the factorization of QED effects in non-leptonic $B$-meson decays. In this paper we study the renormalization group evolution of this non-perturbative function. For charged mesons, in particular, this exhibits qualitative differences with respect to the well-known scale evolution in QCD only, especially regarding the endpoint-behaviour. We analytically solve the evolution equation to first order in the electromagnetic coupling $\alpha_{\rm em}$, which resums large logarithms in QCD on top of a fixed-order expansion in $\alpha_{\rm em}$. We further provide numerical estimates for QED corrections to Gegenbauer coefficients as well as inverse moments relevant to (QED-generalized) factorization theorems for hard exclusive processes.


Introduction
The study of QED effects in B-meson decays has become an active field of research in recent years [1][2][3][4][5][6][7], as such effects will become increasingly more important when higher experimental precision is reached. The standard treatment of QED effects assumes the mesons to be point-like to very small distances of order 1/m B . However, this assumption neglects structure-dependent contributions; photons with wavelength of the order of, or smaller than the typical size 1/Λ QCD of a meson can resolve its partonic substructure. A way to systematically incorporate these effects is to disentangle the relevant energy scales below m B using effective field theories. This allows for a clear separation of very lowenergetic photons with a point-like coupling to mesons from photons with energy of order Λ QCD up to m B , at which scale one can use the 1/m B expansion. First studies following this path were made in [1][2][3][4]. In particular, in [3] we have shown that the QCD factorization formula [8,9] for non-leptonic charmless B decays into two light mesons,B → M 1 M 2 , can be extended to include QED corrections to all orders in α em . Similar to the well-known QCD case, the relevant hadronic matrix elements can then be expressed as [3] 1 This generalized factorization theorem expresses the matrix elements of weak effective operators Q i in the heavy-quark limit as convolutions of hard-scattering kernels T i with lightcone distribution amplitudes (LCDAs) Φ of heavy and light mesons. Besides the perturbatively calculable hard-scattering kernels discussed in detail in [3], a better understanding of adequately generalized non-perturbative objects is required. Therefore, the purpose of the present paper is to study the QED-generalized LCDA Φ M (u; µ) for light mesons M = π, K, . . .. Compared to their definition in QCD-only, LCDAs for electrically charged mesons M ± in QED exhibit some qualitatively new features which can be partly attributed to the non-decoupling of soft photons from the net charge of the meson, but also to the different electric charges of its constituents. First, this concerns the definition of a renormalizable function Φ M (u; µ), as its naive ultraviolet (UV) scale evolution is plagued by infrared (IR) sensitivity that needs to be removed by a "soft rearrangement". The thus properly defined LCDA can be used to derive its anomalous dimension and study the properties under scale variation using the renormalization group (RG). We find that the evolution kernel including QED effects contains, in addition to the QCD ERBL kernel [10][11][12], new local logarithmic terms which have important consequences. For example, the evolution of Gegenbauer coefficients is no longer diagonal, and even the norm is no longer conserved, which spoils the interpretation of Φ M (u; µ) as a probability distribution. Although we do not aim at an analytic solution of the full QCD×QED renormalization group equation (RGE), we provide numerical solutions and an analytic expression at O(α em ) that resums large logarithms in QCD on top of a fixed-order expansion in the electromagnetic coupling.
QED corrections induce isospin-symmetry violation due to the different electric charges of up-and down-type quarks. For non-leptonic B decays, this is of particular interest as it mimics short-distance electroweak penguin contributions which serve as a probe for new physics [3]. In the present work, we find that the distribution Φ M (u; µ) favours larger momenta of the u quark, due to its larger electromagnetic coupling, resulting in a slightly asymmetric function even in the limit of massless quarks. This also modifies the endpoint behaviour of the LCDA, which eventually leads to an ill-defined evolution in the formal limit µ → ∞. Despite these new features, and contrary to the QED generalization of the B-meson LCDA [13], the light-meson LCDAs (almost) retain their universality and are thus relevant to a variety of different hard exclusive processes. In this paper, we study these new properties qualitatively, and estimate them quantitatively for realistic applications.
The outline of the paper is as follows. In Section 2 we give the relevant basic definitions, including a brief review of the soft rearrangement required to define a renormalizable LCDA. Section 3 then states the evolution kernel, for which we study the endpoint behaviour of the LCDA in Section 4. The analytic first-order O(α em ) solution in Gegenbauer moment space is presented in Section 5. We provide some numerical estimates of QED corrections to the Gegenbauer coefficients as well as inverse LCDA moments relevant in factorization theorems for hard exclusive processes in Section 6. We conclude in Section 7. We give technical details on the soft rearrangement and the endpoint behavior in two appendices.

Basic definitions
In QCD, light-cone distribution amplitudes for light mesons M are well-established nonperturbative but universal objects, which appear in the theoretical description of hard exclusive processes at large energies. They are defined as hadronic matrix elements of non-local operators composed of two light-like separated quark fields. The twist-2 LCDA φ M (u; µ), M (p)|q 1 (tn + )[tn + , 0] / n + 2 (1 − γ 5 )q 2 (0)|0 = in + p 2 1 0 du e iu(n + p)t f M φ M (u; µ) , is usually the leading contribution in the twist expansion. In the above definition, the lightlike reference vectors n µ ± , obeying n 2 ± = 0, n + · n − = 2, are conveniently defined through the meson's momentum p µ = En µ − + m 2 M /(4E) n µ + in the frame where E m M is the large energy of order the hard scale of the process. The function φ M (u; µ) is, however, boost-invariant. The displaced fields are connected by a straight light-like Wilson line [tn + , 0] of finite length to ensure gauge invariance of the definition in (2). It is sometimes convenient to express the finite-distance Wilson line in terms of infinite Wilson lines as [tn Finally, f M is the scale-independent meson decay constant in QCD, defined through the local limit t = 0 of (2). This implies the normalization condition 1 0 du φ M (u; µ) = 1.
In the present work, we study the QED generalization Φ M (u; µ) of the leading-twist LCDA φ M (u; µ). By this we mean that the matrix element is computed with L QCD+QED , which accounts for an arbitrary number of virtual collinear photon exchanges between the (electrically charged) constituents of M on top of the strong interaction. For neutral mesons M 0 , one has Q q 1 = Q q 2 = Q q , and the definition (2) remains valid after modifying the Wilson line to include the photon field A µ : where Q q denotes the electric charge of the quark field q in units of e = √ 4πα em . The situation is, however, different for electrically charged mesons M ± . For definiteness, we consider q 1 = D = d, s and q 2 = u, such that M has total charge Q M = Q q 1 − Q q 2 = −1. The corresponding case for M + is related by CP invariance of QCD and QED. The gaugeinvariant bilinear non-local operator now takes the form We first notice that gauge invariance dictates the operator to extend on the infinite lightray, instead of being localized on a finite interval [tn + , 0]. This can be seen by combining the Wilson lines associated with the quark fields of different electric charge to . In addition to the gauge-link [tn + , 0] (d) associated with the charge Q d , the operator contains a Wilson line extending from −∞ to 0 with the total electric charge of the meson Q M = Q d − Q u , defined as where the QCD part of the Wilson line cancels in W (Q M ) (x) due to the meson M − being a colour-singlet. As discussed in detail in [2,3], due to the non-decoupling of soft photons from electrically charged mesons, the operator (5) itself is no longer renormalizable in the sense that its anomalous dimension is IR divergent. This is due to a non-trivial overlap between the soft and collinear sector in QED-generalized collinear factorization theorems. This overlap renders the UV divergences of the collinear part of the operator dependent on the IR regulator. To remove this overlap and make the operator renormalizable, it is sufficient to multiply the collinear operator with a remnant of the soft function of the process which must be present if the entire process is to be IR safe. One particular choice-inspired by the soft function for a decay of a neutral particle into two back-to-back charged particles-is to define subtraction factors R c and Rc through the following vacuum matrix element of soft Wilson lines: The soft Wilson lines originate from the coupling of soft photons to the electrically charged constituents of the particles in the process. For outgoing antiquarks with electric charge Q q , the soft Wilson line reads n ± must be used for outgoing quarks with electric charge Q q . As for the collinear case, the soft Wilson line depending on the total electric charge Q M in (7) is defined by The factors R are defined such that their UV divergences only depend on the IR-regulator associated with the corresponding collinear direction. This cancels the regulator-dependent terms in the anomalous dimension of the collinear operator, leaving a renormalizable operator. We have taken the absolute value of the matrix element (7) to avoid spurious imaginary terms due to soft rescattering phases in the collinear sector. More details and explicit expressions for the subtraction factors in an off-shell regularization scheme can be found in Appendix A as well as the next subsection. Hence, we define the QED-generalized LCDA for an electrically charged light meson as (10) We emphasize that we choose to normalize Φ M (u; µ) with respect to the renormalizationscale independent QCD decay constant f M in the absence of QED. We do not pull out the local limit of the operator in QED, since it would mix into higher logarithmic moments under renormalization group evolution. Hence, also the normalization condition is no longer fulfilled at any scale µ, 1 0 du Φ(u; µ) = 1. The QED generalized LCDA defined in (10) is renormalizable and has a well-defined UV scale evolution. It appears as part of the factorization formula (1) for the IR divergent non-radiative amplitude that describes virtual photon (and gluon) interactions above the strong interaction scale Λ QCD . The complete process is IR finite once real photon radiation of undetected photons with energies below a sufficiently small resolution ∆E is accounted for. The above LCDA is relevant when the energy of the radiated photons is much smaller than Λ QCD ("ultrasoft"), or of the order of Λ QCD ("soft"), in a frame where the electrically charged particles are ultrarelativistic. The LCDAs are defined through exclusive matrix elements and are hence themselves IR divergent. Their proper non-perturbative definition contains a prescription for subtracting these divergences. In the complete description of the process they then appear as IR-finite matching coefficients for the effective theory of ultrasoft radiation. In the application to non-leptonic decays [3], the matching was performed in dimensional regularization, and hence the IR divergences in the LCDA were minimal subtracted, resulting in an IR-subtraction scale dependence, which is not the subject of the present paper (see also discussion at the end of Sec. 5.2. of [3]). An important remark is to be made about the universality of Φ M (u; µ). Although we defined the subtraction factors R c and Rc in (7) by implicitly assuming one additional anticollinear direction n µ + , Φ M (u; µ) remains a universal object relevant to various two-and multi-body exclusive processes. This can be seen by performing the soft rearrangement for a generic n-jet SCET operator: For every light-like direction n i− one defines a corresponding back-to-back vector n i+ with n i+ n i− = 2. The vector n i+ does not have to coincide with the direction of flight of any particle involved in the process. Then the subtraction factors R c i and Rc i are defined through vacuum matrix elements of two soft Wilson lines along n i− and n i+ as in (7). The sodefined collinear sectors are individually renormalizable and so is the left-hand side of (11). Thus, by consistency of SCET as an effective theory, also the rearranged soft operator is renormalizable. The definition implies the choice of a "soft reference frame", which for B meson decays is naturally the B meson rest frame. An unavoidable consequence of the QED effects is the breaking of boost invariance in the LCDA definition, which results in a dependence of Φ M (u; µ) on the large energy E of the meson measured in the soft reference frame. Finally, we remark that once isospin breaking effects are considered, neutral π mesons are described by two distinct quark LCDAs. These can either be defined by Φ M 0 , or alternatively as the SU (2) singlet and triplet linear combinations. In general, the different LCDAs for neutral mesons will mix under renormalization, including a two-gluon LCDA. We do not compute this mixing in this paper, because we are mainly interested in electrically charged mesons, where the non-decoupling of soft photons is important.

Renormalization
The renormalized non-local operator O ren (u; µ) is related to the bare operator through To derive the Z-factor of the operator in (5), we calculate the diagrams in Fig. 1 with slightly off-shell external quark states, and add the MS renormalization factors for the quark fields. The resulting UV poles in dimensional regularization with space-time dimensions d = 4 − 2 have already been given in (21) and (22) in [3]. For completeness, we repeat the O(α em , α 0 s ) result (with slightly adapted notation) wherex ≡ 1 − x and for simplicity we dropped the µ in the argument of Z. Here and below, α em (α s ) denotes the electromagnetic (strong) coupling in the MS scheme at the scale µ.
The plus-distribution is defined in the variable u: As discussed earlier, the UV divergences depend on the off-shellness of the external quark states, k 2 q 1 and k 2 q 2 , inconsistent with renormalization. This dependence is removed by multiplying with R c , which at O(α em ) reads: Here δ c is the remnant of the off-shell regulator in the soft matrix element (see Appendix A for more details). Since soft Wilson lines obey the multiplication law S †(q 1 ) , the constraint δ c = k 2 q 1 /(n + k q 1 ) = k 2 q 2 /(n + k q 2 ) must be imposed on a consistent off-shell regularization prescription.
After this rearrangement the operator is renormalizable and we can compute its anomalous dimension via Then Φ(u; µ) obeys the renormalization group equation (RGE) From the one-loop Z-factor in (13), including now the standard one-loop QCD renormalization, after multiplication with (15), we obtain from (16) the one-loop anomalous dimension Note that this result holds for both charged (q 1 = d, q 2 = u, Q M = −1) and neutral (Q q 1 = Q q 2 = Q q , Q M = 0) mesons. In the latter case the local contributions in the second line of (18) vanish and we recover the standard ERBL evolution kernel [10][11][12] with a modified one-loop coefficient For charged mesons, however, the second line in (18) implies additional logarithmically enhanced local terms. In particular, the anomalous dimension now depends on the large energy E of the meson, which is the hard scale of the process. It is important to emphasize that this energy is measured in a frame that defines the soft modes of the process, i.e. a frame in which the meson is ultrarelativistic. This can be viewed as a manifestation of the factorization anomaly [14] (or collinear anomaly [15]) in SCET, which requires the explicit breaking of boost invariance, in this case through the soft rearrangement. The energy dependence is independent of the meson constituents and only related to its overall charge, which can be seen by comparing to the corresponding collinear matrix element Z (µ) for a point-like fermion with electric charge Q = Q M and energy E = E. The precise definition of Z as well as its one-loop UV divergences can be found in equations (62) and (63) of [3]. Its anomalous dimension and RG evolution reads with solution Z (µ) = U (µ, µ 0 )Z (µ 0 ). The multiplicative evolution factor 2 is We note that in the QED-generalized factorization theorem for non-leptonic B-meson decays the ratio Φ M (u)/Z naturally appears [3], and we will thus normalize the LCDA accordingly in Sec. 5. The explicit quark charges in the second line of (18) violate isospin symmetry and hence the new local terms render the LCDA asymmetric in u ↔ 1 − u.
Further, an expansion in Gegenbauer polynomials will no longer diagonalize the kernel.

Endpoint behaviour
In this section, we study how the QED contributions to the anomalous dimension change the well-known linear endpoint behaviour of the light-meson LCDAs in QCD. The linear endpoint behavior can be inferred from conformal symmetry arguments, see e.g. [16]. This approach uses a construction of conformal operators at the RG fixed point, where the QCD β-function vanishes, which dictates the form of the anomalous dimension. In QCD×QED there is no related RG fixed point since the β-functions depend on both α s and α em .
Therefore, the anomalous dimension is not restricted by these arguments. The following analysis of the endpoint behaviour is based on the one-loop kernel and may not apply at (uninterestingly) large scales, when the QED coupling becomes strong.
To study solutions of the RGE, it is useful to rewrite the plus-distribution in the anomalous dimension (18) as a distribution in the variable v: . . .
We then split the integral in v in the RGE into different momentum regions in order to construct an asymptotic expansion of Φ M (u; µ) near the endpoints u = 0 and u = 1.
Throughout the rest of this section we focus on the limit u → 0; the behaviour near u → 1 is qualitatively similar and follows from the replacements of charge factors, (17) receives contributions from two regions: the soft region, v ∼ u 1, and a "true" collinear region, v ∼ 1, u 1. In these regions the integral of the plus-distribution part simplifies to Collinear region: Soft region: Based on power counting, the plus distribution has been omitted in the collinear region, while the expansion in the soft region implies that the evolution kernel acts on a function space with support on the whole positive real axis [0, ∞). After these simplifications, the integrals in the collinear and soft region may become divergent for v → 0 and v → ∞, respectively, and need to be regularized. However, in the expansion by regions [17], the integrand determines the power-counting, and we can neglect regions that give a suppressed contribution irrespective of whether the integral converges or not. In order to determine which of the two regions dominates we assume that Φ M (u; µ 0 ) ∼ u b for the small-u behaviour at the initial scale µ 0 . We then distinguish the following cases: (17) is dominated by the collinear region which scales as u 1 . The soft region as well as the local terms in the RGE count as u b u 1 and can be dropped. Hence, an infinitesimal evolution µ 0 → µ 0 + dµ generates a term Φ M (u; µ 0 + dµ) ∼ u 1 which now dominates over u b . It follows that, for b > 1, RG evolution in the collinear region always drives the LCDA immediately back to linear asymptotic behaviour.
• For b = 1 the soft and the collinear region contribute equally and there is no apparent simplification of the evolution kernel. We will argue below that we can nevertheless determine the asymptotic form of Φ M (u; µ) from the soft approximation.
• For b < 1 the soft region dominates and is of the same power u b as the local terms in the RGE (expanded for u → 0). This implies that the endpoint behaviour of Φ M (u; µ 0 + dµ) is now fully determined by the endpoint behaviour of Φ M (u; µ 0 ).
• Lastly, for b ≤ −1 the convolution integral in the RGE (17) is ill-defined. Interestingly, for charged mesons we find that RG evolution to extremely large scales inevitably drives the solution to this scenario due to the local terms and the coefficient of the ERBL kernel in (18).
To make the discussion more transparent it is instructive to first analyze the well-known evolution via the QCD-only one-loop kernel along these lines. Assuming an initial condition at the scale µ 0 with b < 1, the asymptotic soft evolution kernel takes the form Interestingly, up to a constant in the local part, this precisely reproduces the evolution kernel [18] for the B-meson LCDA φ + B (ω). This is intuitive, since for u → 0 the D quark in the M − has only soft fluctuations and the large momentum component n + kū of the anti-u quark becomes frozen. To analyze the asymptotic behaviour we make use of some techniques developed in [19,20]. It turns out to be convenient to study the RG evolution in Mellin space: where c is a real parameter. The definition obviously holds also for the QCD-only LCDA φ M (u; µ). If we assume some µ-dependent exponent Φ M (u; µ) ∼ u bµ for u → 0, the Mellin transform converges for Re(η) < b µ , so that c < b µ must be chosen. Without QED, the Mellin-space RGE is where H n = γ E + ψ(n + 1), with ψ(n) = Γ (n)/Γ(n) the digamma function, is the harmonic number function. This equation is solved by [21] The Mellin variable η of the initial condition on the right-hand side shifted by the evolution variable which is always negative for evolution to higher scales, that is, a < 0 for µ > µ 0 . Here the one-loop QCD β function is where N c = 3 and n f is the number of active quark flavours. Eqs. (25) and (27) require the contour parameter c to lie in the interval −1 − a < c < b µ ≡ b − a(µ, µ 0 ) such that the integral converges. The asymptotic behaviour of φ M (u; µ) for small u is then determined by the location of the left-most pole on the real axis for Re(η) > c of the solution (27) after closing the integration contour in the right half-plane. The initial condition φ(u; µ 0 ) ∼ u b implies that its Mellin transform has a singular point at η = b and thus, the shifted function in (27) at η = b µ . As long as b µ < 1, the analytic structure ofφ M (η + a; µ 0 ), rather than the Γ(1 − η) in the prefactor in (27), determines the asymptotic behaviour of the evolved function to be φ M (u; µ) ∼ u bµ , see Appendix B for more details. The precise functional form, i.e. whether or not this power-law is modified by additional powers of ln u, depends on the initial condition. For the rest of this section we restrict ourselves to pure power-like initial conditions. With increasing µ one will eventually reach the point b µ = 1. (Recall that for exponents greater than one the collinear region dominates and immediately generates a linear term.) Due to the gamma function Γ(1 − η) in the Mellin-space solution (27), one might conclude that φ M (u; µ) ∼ u 1 . This is true, although the asymptotic kernel (24) is not the right object to begin with, since the collinear region is of the same order in power-counting as the soft region. As this collinear region does not contribute additional ln u enhanced terms, the soft region produces the correct asymptotic form of the LCDA. This becomes apparent from the expressions in (22). The integral in the soft region in (23), on the other hand, does contribute such terms due to the logarithmic integral We summarize this discussion in the left plot of Fig. 2 by showing the RG flow of the exponent b µ = b − a(µ, µ 0 ) for a given b at the reference scale µ 0 . In QCD-only, a power-like initial condition φ M (u; µ 0 ) ∼ u b , evolved to higher scales µ > µ 0 , results in an asymptotic expansion φ M (u; µ) ∼ u bµ as long as b µ < 1, and φ M (u; µ) ∼ u 1 for b µ ≥ 1. But once the linear endpoint behaviour is reached, the LCDA remains linear under further upward scale evolution, in agreement with the asymptotic form φ M (u; µ → ∞) → 6uū. Initial conditions with b > 1 become automatically linear after one infinitesimal evolution step. Thus, linear endpoint behaviour can be viewed as a UV fixed point of the RG flow in QCD. For electrically neutral mesons the QED anomalous dimension has the same form as in QCD with the replacement α s C F → α s C F + α em Q 2 q , and the above conclusions also hold (see also (43) and discussion thereafter). However, since QED becomes strongly coupled at (phenomenologically uninterestingly) large scales, the perturbative analysis breaks down and results based on the one-loop kernel are no longer reliable.
We will now follow the same steps to analyze to which extent the additional local logarithmic QED terms in (18) alter this result for electrically charged mesons. As we are only interested in the behaviour for small u, we normalize the LCDA to the point-like limit (19) to get rid of some u-independent local terms in the anomalous dimension. We define The asymptotic RGE in Mellin-space forΦ M (u; µ) is with the Mellin transform of the normalized LCDAΦ M (η; µ). Adapting [20], the equation is solved by the ansatẑ The integrand of the exponent in the second line is proportional to the difference of charge factors on the left-hand side and right-hand side of (31), which arises as a consequence of the new ln u term in the evolution kernel. The generalized evolution variable reads In addition, we have to introduce a second variable which globally multiplies the solution and is thus irrelevant for the functional form in u: It remains to understand the analytic properties of (32) in the complex η-plane. The integral in the exponent in the second line can be solved analytically for some special cases that are instructive to study: i) for QED-only with one-loop running of α em (µ) and ii) in the approximation of scale-independent gauge couplings in QCD×QED. In both cases, we find a simple expression that is very similar to the QCD solution in (27): but now with the combination of Gamma functions raised to the non-integer exponent Let us discuss the QED-only case i) first, for which p reduces The important difference to QCD-only is that the evolution variable is always positive for an upward scale evolution to µ > µ 0 . This is due to the quark charges and the sign of the one-loop QED β-function where n u (n d ) denotes the number of the active up (down) quark flavours and n the active lepton flavours. The endpoint behaviour is then Φ M (u; µ) ∼ u bµ for b < 1, and Φ M (u; µ) ∼ u 1−a(µ,µ 0 ) for b ≥ 1. In particular, the exponent b µ is driven towards smaller values, and, since a increases without bound, will approach the point b µ = −1 for large enough µ. This not only means that the LCDA becomes divergent at the endpoint, it also implies that RG evolution pushes the solution towards a functional form that is no longer compatible with the RGE itself, since for b µ ≤ −1 the convolution integral in the evolution equation is no longer well-defined. Thus, QED evolution inevitably drives the solution outside the validity of its evolution kernel.
In the case ii) of fixed couplings in QCD×QED, the evolution variable is given by The QCD and QED terms enter with different signs due to the quark and meson charges. The behaviour of the solution (35) now depends on whether the overall sign of (38) is positive or negative. For α s C F + α em Q q 1 (Q q 2 − Q M ) < 0, the evolution variable a is positive and b µ decreases with increasing µ. The endpoint behaviour follows the same pattern as in the previously discussed QED-only case i). For the phenomenologically more relevant situation α s C F + α em Q q 1 (Q q 2 − Q M ) > 0, the evolution variable a is negative and b µ increases with µ. The situation is then similar to QCD-only, namely that the solution is pushed towards linear endpoint behaviour with Φ M (u; µ) ∼ u bµ for b µ < 1, eventually reaching b µ = 1 at some finite µ < ∞. However, for b µ = 1, one additional complication arises: the singularities of the gamma functions in (35) now have an associated branch cut.
In particular, the left-most single pole from Γ(1−η) has turned into a cut. After integration along this cut, we obtain Hence, the linear endpoint behaviour is modified by a logarithmic term raised to the noninteger power p. In the limit α em → 0 the exponent p vanishes and we recover the wellknown linear endpoint behaviour in QCD. Expanding in α em to first order generates a term proportional to α em u ln (− ln u). Now, in the general case of QCD×QED with one-loop running coupling the integrand of (33) naturally defines a critical scale by At this scale, the exponent p(µ c ∓ 0) = ±∞ becomes singular and flips sign. This changes the analytic properties of the Mellin-space RG solution. The resulting RG flow of b µ is qualitatively shown in the right plot of Fig. 2. Similar to case ii), it is useful to discuss the general solution (32) for µ 0 > µ c and µ < µ c separately. For µ 0 > µ c , the lefthand side of (40) becomes negative, and the above discussion for fixed couplings with α s C F + α em Q q 1 (Q q 2 − Q M ) < 0 applies, so that the endpoint behaviour is simply given by Φ M (u; µ) ∼ u bµ with decreasing b µ < 1. Independent of µ 0 and the functional form at this scale, the point where evolution becomes inconsistent is always reached at some finite scale µ L > µ 0 . At these extremely high scales, one enters the strong QED coupling regime and therefore this case is not relevant for realistic scenarios. However, in principle, in theories with a different flavour content where the magnitude of β QED 0 is smaller, the evolution of Φ M (u; µ) can break down before entering the strong coupling regime. For µ < µ c , the endpoint behaviour for increasing b µ < 1 is again given by Φ M (u; µ) ∼ u bµ . Based on the observation for fixed couplings, the particular point of interest is b µ = 1. Remarkably, the result (39) also holds in the general case of scale dependent gauge couplings in QCD×QED. We prove this statement explicitly in Appendix B.2. An important consequence is that the inverse moments relevant for (QED-generalized) factorization theorems remain integrable.
To qualitatively visualize the asymmetry of the QCD×QED LCDA generated by evolution, we show in Fig. 3 a numerical solution of the RGE obtained from solving the integro-differential evolution equation (17) by discretization (see Section 6 for more details) for fictitious values of the electromagnetic coupling constant.  Figure 3: Numerical solution to the RGE for the π − LCDA by discretization in u using N = 1001 logarithmically distributed points (see Section 6 for more details). We model the LCDA by the QCD-only initial conditionΦ π − (u; µ lat ) = 6uū(1+a π 2 (µ lat )C (3/2) 2 (2u−1)) at the reference scale µ lat = 2 GeV (black solid curve), using the lattice value for the Gegenbauer moment a π 2 (µ lat ) given in Table 1. The dotted curves show the function evolved to µ = 10 GeV for three fictitiously large values of α em . The QED coupling is fixed, while α s (µ) runs at one-loop. The following qualitative features induced by QED effects can be observed: First, the norm is no longer conserved. Second, the distribution becomes asymmetric and favours the up quark, which carries a larger fraction of the pion momentum due to its larger electromagnetic coupling. Third, the endpoint behaviour gets modified, and above a certain threshold the LCDA starts to diverge.

Gegenbauer coefficients and analytic O(α em ) solution
In the real world, QED effects are expected to be small. In this section we provide an analytic solution to the QCD×QED RGE that treats QED effects at first order while summing ln µ/µ 0 to all orders, i.e. a solution accurate to α k em α n s ln n+k µ/µ 0 with k = 0, 1. For this purpose, we expand the LCDA in Gegenbauer polynomials as usual such that the integro-differential evolution equation (17) becomes an infinite dimensional system of ordinary differential equations of Gegenbauer coefficients a M n (µ). To separate universal from structure-dependent QED effects we again normalize to the point-like limit and consider where C (3/2) n (2u − 1) are Gegenbauer polynomials and a M n (µ) the n-th Gegenbauer coefficients for meson M . The RGE for the coefficients a M n (µ) reads For neutral mesons, Q q 1 = Q q 2 = Q q and Q M = 0, the scale evolution of the a M n (µ) is diagonal and solved by a n (µ) = α s (µ) α s (µ 0 ) in the leading-logarithmic (LL) approximation. The n-dependent anomalous dimension reads [10] γ n = 1 − 2 (n + 1)(n + 2) and grows logarithmically for large n, γ n ≈ 4 ln(n). In QCD, the evolution to larger values of µ suppresses higher Gegenbauer coefficients since α s (µ) < α s (µ 0 ) and the exponent C F γ n /β QCD 0 in (43) is positive and growing with n. The QED factor in (43) additionally suppresses the coefficients as α em (µ) > α em (µ 0 ) and Q 2 q γ n /β QED 0 is negative. This justifies the truncation of the Gegenbauer series at some fixed value n = n 0 for neutral mesons even after QED corrections are included. In QCD only, one typically chooses n 0 = 2, up to which data from lattice QCD is available. Note that γ 0 = 0, so a 0 is not renormalized for neutral mesons. In fact, the normalization condition 1 0 du Φ M (u; µ) = 1 implies a 0 = 1. For charged mesons we observe two important differences. First, the QED contribution enters the first term in (42) with a different sign. However, at physically relevant scales, α em (µ) can be viewed as a small perturbation compared to α s (µ), justifying a truncation also in that case. 3 The solution to (42) can then be obtained by numerically solving the resulting finite-dimensional system of first-order differential equations. Second, for Q M = 0 the local ln u and lnū terms in the anomalous dimensions cause a mixing of Gegenbauer coefficients under RG evolution, with the infinite-dimensional mixing matrix f nm given by Contrary to the triangular structure of the anomalous dimension matrices in QCD, the mixing matrix f nm has no particular structure. This means that in QED the Gegenbauer coefficients a n can also mix into lower coefficients a m with m < n. Since f 0m = 0, there is even a mixing of all Gegenbauer coefficients into the zeroth coefficient a M 0 (µ), such that the standard normalization condition of the LCDA is no longer valid for charged light mesons due to scale evolution, i.e. 1 0 du Φ M (u; µ) = 1. The diagonal terms f nn approach the constant Q M ln 4 and are thus negligible compared to the logarithmically growing first term in (42) for large n. In all other cases, f nm falls off for large n or m, indicating that the mixing of Gegenbauer coefficients with largely different n is strongly suppressed (even for large values of α em (µ)). For example, for n m, the mixing matrix drops like f nm ∼ 1/n 3 , and for m n one has f nm ∼ 1/m 2 . Also for n and m large but of the same order the off-diagonal terms fall off like 1/n.
For scales relevant to hard exclusive processes it is sufficient to consider α em (µ) as a small perturbation. In the following we provide an analytic solution to the RGE for the Gegenbauer coefficients to first order in α em . In other words, we sum large logarithms L to all orders in QCD and count α s × L ∼ O(1), but retain the fixed-order expansion in α em , since α em ×L 1 is still small. This corresponds to the summation of the leading logarithms α k em α n s ln n+k µ/µ 0 with k = 0, 1. However, we always retain the resummed expression for the large QED double-logarithms in (20) associated with the point-like limit. To this end, we expand the Gegenbauer coefficients as a n (µ) = a QCD The QCD moments a QCD n (µ) ∼ O(1) count as order one. If one expands the initial condition at µ 0 in the above form, a (1) n (µ 0 ) ∼ O(1), and the QED initial condition is a correction beyond the LL order. However, as will be seen from the solution below, evolution generates a QED correction a (1) n (µ) ∼ O(ln(µ/µ 0 ), α s (µ 0 )/α s (µ)), which identifies the second term on the right-hand side of (46) as the LL term linear in the QED coupling. The RGE in (42) evidently reduces to the standard QCD evolution equation at zeroth order in O(α em ), and the inhomegenous equation for the first-order coefficients in the electromagnetic coupling α em . The inhomogeneous term reads and depends on all QCD Gegenbauer coefficients at the scale µ. The index M indicating the light meson has been dropped on the a n (µ) for better readability. The solution to (47) recovers the standard LL QCD expression The desired LL solution for the linear QED correction reads As mentioned above, the QED initial condition a (1) n (µ 0 ) is technically beyond the LL accuracy. Also, if the input values at the scale µ 0 are given by a certain model or (future) lattice calculations in QCD×QED, it does not seem very natural to expand the initial condition itself in α em . In both cases, one sets a (1) n (µ 0 ) → 0 in (51), as we do in the numerical analysis in the next section.

Numerical estimates
In this section we provide more details on the discretization of the RGE used to obtain Fig. 3 and give numerical estimates of the QED corrections to the Gegenbauer coefficients and the inverse moments of which are relevant to (QED-generalized) factorization theorems for hard exclusive processes. We recall that the LCDA (but not its evolution) is IR divergent for electrically charged mesons, as it is part of the non-radiative amplitude and only accounts for virtual QED contributions. The initial conditions employed below should be interpreted as a model for the properly IR-subtracted, and thus scheme-dependent, LCDA at the scale µ 0 . Under this assumption, the UV scale evolution can be studied for which we give numerical estimates. We use the input values given in Table 1. For the gauge couplings, we decouple the bottom (charm) quark in α s (µ) at its pole mass µ = m b (µ = m c ). For the electromagnetic coupling α em (µ), we also include the threshold at m τ . To discretize the evolution equation (17), we divide the interval u ∈ [0, 1] into N − 1 sub-intervals by distributing N points u i , with i = 1, . . . , N , between 0 and 1. For large N the integral in v can be approximated by a Riemann sum and the integro-differential equation turns into a coupled system of N ordinary first-order differential equations which Gegenbauer coefficients at µ 0 = 1 GeV Table 1: Numerical inputs. The Gegenbauer coefficients at µ lat = 2 GeV are the lattice QCD results from [22]. The quark masses are to be understood as two-loop pole masses.
we solve numerically. To increase the accuracy, we use the trapezoidal rule, whose error is roughly proportional to the third power (∆u) 3 of the difference ∆u of two points. As the integration kernel is divergent at the endpoints u = 0 and u = 1, we shrink the interval u ∈ [0, 1] to u ∈ [ , 1 − ], with 1. Especially when studying the endpoint behaviour, it is important to choose sufficiently small. As a default value we use = 10 −10 . To increase the accuracy it is useful to distribute points non-uniformly in such a way that the density of points is logarithmically enhanced towards the endpoints. Using, for example, N = 1001 points, we reproduce the known analytic LL solution in QCD with an error of less than 0.04% for the initial condition and scale choice described in the caption of Fig. 3. We emphasize that this error analysis serves as a cross-check of the implementation, but provides only a rough estimate of the expected accuracy in QED, as contributions from the endpoints can be enhanced for large, but unphysical values of α em .
To estimate the size of QED contributions, we compare the scale evolution of the Gegenbauer coefficients in pure QCD at the LL, NLL, and NNLL order with the evolution in QCD×QED using the first-order QED solution (51). 4 For the NLL evolution, we use the NLO anomalous dimension matrix and the two-loop running coupling α s , whereas the NNLO anomalous dimension and three-loop running is used for the NNLL results. In QCD, we use the input values for the Gegenbauer coefficients at µ lat = 2 GeV obtained from lattice QCD [22] given in Table 1, and evolve to a higher scale µ. In QCD×QED, we first run down to the natural scale µ 0 = 1 GeV at LL in pure QCD, resulting in the values given in Table 1 and then evolve to µ including QED. We give numerical results for two different scales: µ = 5.3 GeV ≈ m B and µ = 80.4 GeV ≈ M W . The first is relevant for exclusive B decays, the latter would be relevant e.g. for the very rare W − → π − γ decay [23].
In QED all Gegenbauer coefficients at the low scale µ 0 with n > 2 mix into a 0,1,2 (µ). However, as currently the Gegenbauer coefficients with n > 2 are unknown and since they are expected to be small, we set them to zero at the initial scale µ 0 . In addition, these higher Gegenbauer coefficients will be generated through scale evolution in QED, and also in QCD beyond the LL accuracy. We include this mixing when calculating the inverse moments, as discussed below. At the moment, also the QED corrections to the Gegenbauer coefficients at the reference scale µ 0 = 1 GeV are unknown. Nevertheless, as discussed above, the dominant logarithmically enhanced (LL) contributions are captured by evolution.
For a negatively charged pion π − = (dū), we obtain at µ = 5.3 GeV the values and for the kaon K − = (sū) At µ = 80.4 GeV We refrain from quoting uncertainties as the current uncertainty on the lattice input values at µ lat = 2 GeV is O(15%), which is larger than the QCD NLL contribution in most cases. Our aim is to compare the relative size of QED effects to higher-order evolution in QCD. We find that the QED effects are almost an order of magnitude larger than the NNLL evolution, although consistently below the NLL contribution. For example, the relative contribution of the QED effects in a π − 2 rises from roughly 1% at 5.3 GeV to almost 4% at 80.4 GeV. Compared to the NLL contribution it increases from roughly 12% to 21%. When evolving to higher scales, the QED effects become relatively more important, because the absolute value of the QCD coupling decreases, whereas simultaneously the QED value increases slightly. In fact, in QCD-only, the value of all Gegenbauer coefficients will tend to zero in this upward scale evolution and the LCDA approaches the asymptotic form φ M (u; µ → ∞) = 6uū. This does not hold for QED. In addition, QED modifies coefficients that are fixed to all orders in QCD by isospin symmetry (e.g. a π − 1 = 0) or by the normalization of the LCDA (a M 0 = 1). Besides the Gegenbauer moments themselves, the inverse moments of the LCDA are important, as they often appears at leading order in factorization theorems for observables. They are defined by Here we separated the effect arising from the point-like limit contained in Z (µ), denoted by "point charge", from the structure-dependent contributions in the Gegenbauer coefficients, denoted by "partonic". Again, we keep the resummed form for Z (µ) in (20) but use the fixed-order O(α em ) solution (51) for the Gegenbauer coefficients. 5 We set the energy E that enters the evolution of Z (µ) to E = µ/2, which is the hard scale in two-body decays. We note that the inverse moments in (57) and (58) depend on an infinite sum of Gegenbauer coefficients, which were set to zero at the initial scale µ 0 for n > 2. Scale evolution, both in QED and QCD at NLL, will then generate these higher Gegenbauer moments, which we include up to n max = 100. We find that the sum for QCD NLL converges rather slowly. However, based on a naive convergence analysis, increasing the number of Gegenbauer coefficients may change the last digit of our results by at most one. The convergence in QED is much better and in almost all cases it was sufficient to truncate at a maximum value n max = 10 to obtain the quoted result. We note that we do not include an additional uncertainty from neglecting unknown Gegenbauer coefficients with n > 2 at the reference scale µ lat = 2 GeV. In this sense, our numerical analysis should be again interpreted as an estimate of the relative effects of higher order evolution in QCD versus QED effects, given a model at the low scale. For the inverse moments the relative size of the various effects is different from those observed for the individual Gegenbauer coefficients. We first note that both QCD NLL as well as QED effects are typically of O(1%) and of the same size as the uncertainty from the lattice values given in Table 1. For comparison, we therefore provide the uncertainty from the lattice Gegenbauer moment input for the QCD LL results, but do not give the corresponding errors on the QCD NLL and QED as well as NNLL contributions here, which are only a subleading correction. The structure-dependent QED effects become larger than the NLL QCD correction from evolution above scales of order O(10 GeV). At the high scale 80.4 GeV it also exceeds the uncertainty from the current lattice determination of the input values. This can be understood from the fact that, unlike QCD, QED evolution couples all a M n to the zeroth Gegenbauer coefficient a M 0 , i.e. the norm of the LCDA, which is an order of magnitude larger than the first and second Gegenbauer coefficients. For example, at µ = 80.4 GeV for the π − , QED corrections cause an O(1%) effect on a M 0 = 1, whereas higher-order QCD evolution is a 15% effect on the LL value a π − 2 = 0.0657, which enhances the relative size of QED effects. In comparison, the isospin breaking effects arising from QED due to a difference ū −1 − u −1 π − are at the few permille to 1% level. Lastly, we note that the structure-dependent QED corrections enhance the value for the inverse moments, whereas the contribution from the point charge acts in the opposite direction. The separation of the two effects is nervertheless useful, since the point-charge contribution is typically factored out and by definition not considered as part of the nonradiative amplitude.

Conclusion
In this paper, we studied the evolution of the leading-twist light-cone distribution amplitude (LCDA) for light mesons in QCD×QED. This QED-generalized LCDAs was introduced as part of the generalization of the QCD factorization formula for non-leptonic decays to QCD×QED [3,4]. We solved the (one-loop) RGE numerically and provided analytical expressions for the RGE at O(α em ) which resum the large logarithms in QCD on top of the fixed-order expansion in α em .
For electrically neutral mesons, the RGE kernel is simply the QCD ERBL evolution kernel [10][11][12] with the modified coupling α s C F → α s C F + α em Q 2 q , where Q q denotes the electric charge of the quark in the meson and the evolution is similar to the case of QCD. For electrically charged mesons, however, the solution to the RGE exhibits interesting, qualitatively novel features, which ultimately arise from the coupling of soft photons to the net charge of the meson: • QED effects change the endpoint behavior of the LCDA as shown in Fig. 3. The linear vanishing of the LCDA near the endpoints u = 0, 1 is no longer the UV fixed point of the evolution.
• Due to new local logarithmic terms in the RGE, the evolution of the Gegenbauer coefficients is no longer of triangular form in Gegenbauer moment space. Opposed to QCD, it induces mixing of higher Gegenbauer coefficients into lower ones. Moreover, the normalization condition a M 0 = 1 is not stable under RGE evolution and can no longer be imposed. This is a manifestation of the scale dependence of the analogue of the decay constant for a charged pion in QCD×QED (defined by t = 0 in (10)).
• The QED evolution violates isospin symmetry and the LCDA of charged π mesons becomes antisymmetric. The QED-generalized LCDA favours larger momentum of the up-type quark due to the larger modulus of its charge.
• The QED-generalized LCDA is not boost-invariant and depends on the energy of the meson, measured in a soft reference frame, determined by the process under consideration. The energy dependence appears in universal double-logarithmic terms, which can be factored out and equal the contribution of an electrically charged point particle.
In addition to the investigation of the RGE solution, we presented numerical estimates of the QED corrections from evolution for the Gegenbauer coefficients and the inverse moments, relevant for hard exclusive processes, at the two scales µ = 5.3 and µ = 80.4 GeV. For the Gegenbauer coefficients, we found that the QED effects are almost one order magnitude larger than NNLL evolution. For the inverse moments, both the QED effects and the QCD NLL are at the percent level, and of similar size as the uncertainties from the lattice input values. For these moments, we separated the QED effects in structuredependent and point-charge terms, which is useful in light of the QCD×QED factorization formulas, where the QED-generalized LCDA naturally appears [3,4].
for which we can calculate the endpoint behaviour explicitly. This choice of cut-off has the advantage that the Mellin transform yields only a single simple pole at η = b in the complex η-plane. The Mellin integral converges for Re(η) < b, so that c < b must be chosen for the inverse.

B.1 Inverse Mellin transform in QCD
In QCD-only, the solution to the soft RGE in Mellin space is now given by (27) together with (70): The evolution variable a is defined in (28), and b µ ≡ b − a(µ, µ 0 ). In the complex η-plane, the integrand is analytic up to strings of simple poles from the gamma functions in the numerator and one pole from the Mellin transformed initial condition. The two strings of poles are located at η = −1 − a − n and η = 1 + n, where n = 0, 1, 2, . . . is a non-negative integer, and the single pole from the initial condition is located at η = b µ . In the limit u → 0 we can safely assume u < Λ. The factor (u/Λ) η then exponentially suppresses the integrand for Re(η) → ∞, such that we can close the integration contour in the right half-plane. The asymptotic behaviour of the gamma functions for large real arguments requires us to choose c in the interval −1 − a < c < min(1, b µ ). This implies the constraints b > −1 and a > −2, which excludes an overlap between the two strings of poles. 7 After deformation, the new contour encircles the poles at η = b µ , 1 + n in the mathematically negative direction, and the contour integral is given by the sum of all residues at these poles. This in particular implies that the left-most pole satisfying Re(η) > c gives the dominant contribution in the limit u → 0: The first term is independent of the regulator Λ. This is expected as for b µ < 1 the collinear region is power-suppressed, and the endpoint behaviour at the scale µ is fully determined by the endpoint-behaviour of the initial condition at the scale µ 0 . The regulator dependence in the linear term, however, is cancelled only after including the collinear region in the evolution. We note that for b µ → 1 the sum of both terms is finite. Hence, depending on whether b µ is smaller or larger than one, the endpoint behaviour will be either linear or u bµ . Since we have a < 0 in QCD, we can conclude that the RG solution always flows towards a linear endpoint behaviour. The same analysis holds for the endpoint u = 1, showing consistency with the asymptotic form of φ M (u; µ → ∞) → 6uū.

B.2 Inverse Mellin transform in QCD×QED
We now turn to the general case of QCD×QED and show how additional non-integer powers of logarithms ln u arise from the contour integral. As the evolution above and below the critical scale µ c defined in (40) yields slightly different results, we analyze them separately. We start with the more relevant case µ < µ c , which implies that the evolution variable a(µ, µ 0 ) < 0 is negative for all scales with µ 0 < µ < µ c . In QCD×QED, the analysis of the soft region is more involved due to the additional exponential of the integral over harmonic number functions in (32). With the initial condition (70), (32) readŝ where b µ is defined in QCD×QED by (33). The analytic structure of this exponential is nontrivial, since the integration variable µ appears in the harmonic number function argument. The harmonic number functions have simple poles located at the negative integers, and consequently the integral has branch cuts on the real η-axis, as discussed below. As in QCDonly, the parameter c for the inverse transformation must be chosen in interval −1 − a < c < min(1, b µ ), separating the cuts into left (due to H η+a(µ,µ ) ) and right contributions (due to H −η−a(µ,µ ) ). Asymptotically, the harmonic numbers behave like ln η for large real arguments, such that for u → 0, we can still deform the contour to enclose the discontinuities on the right-hand side. Again, the left-most pole or cut for Re(η) > c then determines the endpoint behaviour. For b µ < 1, we just pick up the residue at η = b µ as in QCD-only and the resulting endpoint behaviour is u bµ . Once reaching b µ = 1, we need to analyze the solution (73) for the initial condition b = 1 as it is a stable point of the evolution. At η = 1, the exponential factor in the second line of (73) has a branch-point due to the left-most pole of the harmonic number function, H −η−a(µ,µ ) = −1/(1 − η − a(µ, µ )) + regular terms and poles/cuts further to the right. The relevant non-analytic term can be extracted from In the integration domain in µ , the function a(µ, µ ) takes values in the interval [a, 0], with a = a(µ, µ 0 ) < 0. For given µ, the integral in (74) thus has a branch-cut of finite length, starting at η = 1 and extending to η = 1 − a. We can conclude in general that the harmonic number function leads to cuts on the intervals [1 + n, 1 − a + n] with n = 0, 1, 2, . . . a nonnegative integer. To determine the contribution from the left-most branch-cut, we recall that the branch-point at η = 1 requires µ → µ, such that we can expand a(µ, µ ) = − α s (µ)C F + α em (µ)Q q 1 (Q q 2 − Q M ) π ln µ µ + . . .
in the denominator of the integrand in (74), and replace α em (µ ) → α em (µ) in the numerator. The integral in µ then essentially reduces to the scale-independent case already discussed in the main text: − µ µ 0 dµ µ α em (µ )Q q 1 Q M π H −η−a(µ,µ ) = p(µ) ln 1 1 − η + regular terms and poles/cuts further to the right, but with µ-dependent p(µ) = α em (µ)Q q 1 Q M /(α s (µ)C F + α em (µ)Q q 1 (Q q 2 − Q M )). The function p(µ) > 0 is always positive for evolution below µ c ; at µ = µ c , it is singular and flips its sign. The result (76) for the exponential needs to be combined with the pole from the Laurent expansion of Γ(1 − η) in (73). The regular and other singular terms are irrelevant for the asymptotic small-u behaviour, and can thus be absorbed, along with other irrelevant factors, into a constant κ. The inverse Mellin transformation now takes the form where the contour C encloses the cut for Re(η) > 1, as discussed above. According to [25], the contour C must be carefully chosen. We decompose it into a small circle C ε of radius ε around η = 1, that is parametrized by η = 1 + εe iϕ for ϕ ∈ (2π, 0), and the contour C cut consisting of straight lines above and below the real axis extending from 1 + ε to ∞. The contributions from the cut and the circle are given bŷ For p(µ) > 0, the first terms n < p of the sum are divergent as ε → 0, but they cancel exactly after adding (78) and (79). In total, we obtain Φ M (u; µ) =Φ Ccut M (u; µ) +Φ Cε M (u; µ) + O(u 2 ) = κ Γ(1 + p(µ)) u(− ln u) p(µ) + O(u 2 ) , (81) which is also dictated by Watson's lemma. This confirms our statement in (39) for scaledependent gauge couplings. The discussion changes when evolution above µ c is considered. For µ > µ 0 > µ c the evolution variable a(µ, µ 0 ) becomes positive and hence the cut from the harmonic number function starts at η = 1 − a. The relevant contribution in (73) is then either given by the point η = b µ or η = 1 − a. In principle, the latter also gives rise to ln u-modified terms from the harmonic numbers. However, this contribution will never be physically relevant since, independent of the value of b, evolution will immediately generate the endpoint behaviour u 1−a(µ 0 +dµ,µ 0 ) with 1 − a(µ 0 + dµ, µ 0 ) < 1. This justifies to assume b < 1 from the beginning, and consequently b µ < 1 − a < 1, such that the dominant contribution from the inverse transformation of (73) is always given by the residue at η = b µ . We conclude that, above the critical scale, the power-like behaviour is always proportional to u bµ in QED. Hence, at some large but finite µ, the evolution will run to solutions with endpoint behaviour u −1 (b µ = −1), for which the RGE is no longer defined.