Final-state QED Multipole Radiation in Antenna Parton Showers

We present a formalism for a fully coherent QED parton shower. The complete multipole structure of photonic radiation is incorporated in a single branching kernel. The regular on-shell 2 to 3 kinematic picture is kept intact by dividing the radiative phase space into sectors, allowing for a definition of the ordering variable that is similar to QCD antenna showers. A modified version of the Sudakov veto algorithm is discussed that increases performance at the cost of the introduction of weighted events. Due to the absence of a soft singularity, the formalism for photon splitting is very similar to the QCD analogon of gluon splitting. However, since no color structure is available to guide the selection of a spectator, a weighted selection procedure from all available spectators is introduced.


Contents
1 Introduction Phenomenology for high energy particle collisions such as those at the LHC is currently strongly reliant on the implementation of theoretical knowledge in Monte Carlo based event generators [1]. In these programs, the simulation of collisions is partitioned into several tasks, of which one is the parton shower. These programs dress hard events with soft and collinear radiation, resumming the associated leading logarithms, and preparing them for hadronization. The Standard Model allows for many types of radiation, of which QCD radiation is the most important in the context of LHC phenomenology. For that reason, QCD parton showers have been developed for decades. Traditionally, parton showers were based on the DGLAP equation [2][3][4] and 1 → 2 kinematics. Soft coherence can then be ensured by some form of angular ordering [5,6]. More recently, showers based on 2 → 3 kinematics have become increasingly popular. These showers also obey the DGLAP equation, but are automatically soflty coherent and have access to exact phase space factorization. QCD parton showers based on 2 → 3 factorization can be divided into two categories. On the one hand are the antenna showers based on the antenna factorization scheme for fixed order calculations [7][8][9][10][11]. These showers include the ARIADNE code [12] which was very successful in the LEP era, the SHERPA-based ANTS shower [13], and more recently the VINCIA shower [14][15][16][17]. Antenna showers treat the partons participating in a branching on the same footing. From the factorization properties of matrix elements, an antenna function is derived which serves as probability measure for branchings. These branching kernels capture both the soft and collinear divergent behaviour of the matrix element associated with the branchings.
On the other hand are the dipole showers based on the Catani-Seymour factorization scheme [18,19]. These types of showers are currently included in the SHERPA [20] and HERWIG [21] event generators. Dipole showers essentially divide the antenna function into two pieces that radiate independently. For both of these kernels, one of the participating partons acts as a recoiler, while the other parton branches.
Finally, the recently developed DIRE [22] parton shower is a DGLAP-based shower that contains all the advantages of dipole and antenna showers. It achieves soft coherence by using modified Altarelli-Parisi kernels and uses 2 → 3 for exact momentum conservation.
Other types of radiation are allowed in the Standard Model. This radiation should be interleaved with the dominant QCD radiation. Most event generators have implementations for QED radiation [23][24][25], which includes the emission of photons from charged fermions, and the splitting of photons into fermion-antifermion pairs. However, photon emission is not formally coherent in any of these implementation. In leading-color QCD DGLAP-style showers, coherence can be achieved by angular ordering, but this does not straightforwardly extend to QED. In antenna or dipole approaches, the branching kernel cannot be partitioned into independently radiating pieces, since there is no color structure distinguishing their post-branching states. On the other hand, YFS exponentiation [26] is used by some event generators [27,28] to add soft photon radiation to particle decays, and by others to simulate process-specific photon radiation for precision physics [29][30][31][32]. This type of radiation is coherent, but collinear logarithms can only be included order-by-order, and it cannot be interleaved with QCD radiation. In this paper, we introduce a parton shower formalism which produces coherent QED radiation and which can be interleaved in a regular QCD shower. Our formalism does not nessecarily fit into either the antenna or dipole categories, but since its kinematical approach is closer to an antenna shower, and to VINCIA in particular, we refer to it as such.
In section 2 the shower formalism for photonic radiation is introduced. A comparison with leading-color QCD showers is used to explain the complications that have to be dealt with in a fully coherent photon shower. The algorithmic implementation of the shower and a method of improving performance are also explained. In section 3 the shower formalism for photon splitting into charged fermion-antifermion pairs is explained. This type of radiation is very similar to gluon splitting, with the exception of the absence of a color structure to dictate a spectator parton. The shower formalism is tested in various ways in section 4, including comparison with exact matrix element calculations, the DGLAP equation and YFS radiation patterns. Finally, section 5 contains a summary and some outlook.

Notation and conventions
For momenta p i and p j and masses m i and m j we will make extensive use of the notations In a massive 2 → 3 branching, we denote the pre-branching system with upper case letters and the post branching-system with lower case letters.

Photon emissions
All antenna and dipole showers function in the leading-color QCD approximation. In this context, it makes sense to partition the total gluon contribution into antennae or dipoles corresponding with differing color ordered states. Because of the leading-color approximation, the number of contributing antennae or dipoles is limited and their interference structure is automatically disentangled. In contrast, for photon emissions, there is no color structure or leading color approximation. This means that every pair of charged (anti)fermions contributes equally, and there is no way to divide the kernel into several disconnected pieces. As a consequence, every charged fermion in principle participates in the emission of a photon. The implementation of such a procedure in a parton shower formalism would be problematic, especially if it is to be interleaved with a regular QCD shower. For that reason, we will employ an approach similar to the sector shower detailed in [17] to cast photon emissions in an antenna shower-like procedure. Below, we first derive the emission kernel, and then proceed to detail the shower implementation by the definition of the ordering variable and the kinematics.

Emission kernel
Let |M 1 ({p}, k)| 2 be a squared amplitude of a set of n final-state charged fermions and antifermions with momenta p i and charges Q i , and a photon with momentum k. The factorization properties of such amplitudes are well-known [33]. In the soft limit k → 0, the matrix element factorizes into a photonic piece which is usually called the eikonal factor, and a lower order matrix element that does not include the photon: where α is the fine-structure constant. The soft factorization of the eikonal factor is a general property of QED matrix elements and serves as one of the bases of the shower. We further note that the momentum of every charged particle appears in the eikonal factor, but they remain unchanged since the photon is soft. This property must be reflected in the parton shower. Factorization of the matrix element also occurs in the the quasi-collinear limit, which is the massive generalization of the collinear limit [11,34]. For momentum p a it is defined as Compared with the massless collinear limit, the main difference is that m 2 a and m 2 a should be kept of the same order while m 2 a tends to zero. In this limit the matrix element factorizes into a quasi-collinear piece involving only k and p a , and the remaining photon-less matrix element with p a : is the generalized Altarelli-Parisi splitting kernel. Note that in this limit, p a is the only momentum participating in the emission. Again, the parton shower should reflect this property. Similar factorization theorems apply in QCD. However, in the leading-color approximation, only terms involving partons which are adjacent in color space survive. The radiative factors of soft and collinear gluon emission inbetween two partons a and b can be combined into an antenna function where k is the gluon momentum and C is a color factor. This function captures the soft and the collinear factorization properties of a partial amplitude in the leading-color limit. The total gluon emission amplitude can then be approximated by summing over such an antenna function for every possible pair of partons {ab}. The functional form of a QCD e can be computed from matrix elements [10,11] as where X is some decaying particle such as a Z or γ * . Depending on X, the computation yields different non-singular terms as part of the antenna function. These terms serve as tuning parameters in the parton shower.
There is no leading-color limit in QED, and none of the eikonal factors of eq. (2.1) is of lesser importance. It is therefore not feasible to describe the emission of a photon as a sum of antennae without resorting to approximations. We instead capture all factorization properties in a single emission kernel. Similar to QCD, the singular structure for photon emissions is extracted from a matrix element calculation. The total emission kernel up to nonsingular terms is then just the antenna function summed over all charged pairs.
where we note that additional nonsingular terms may be included depending on the process used to calculate the individual terms. While in eq. (2.1) the indices of the sum run over all charged fermions, here the sum runs over all pairs {ab}. The a = b terms in eq. (2.1) lead to the mass-dependent terms in eq. (2.7) which are redistributed amongst the pairwise sum using charge conservation.

Ordering variable and phase space
We now cast the result in the standard antenna shower formalism [12][13][14] including fermion mass effects. At its core lies a 2 → 3 branching step which yields on-shell massive momenta and obeys momentum conservation. The probability for emissions are proportional to the emission kernel of eq. (2.7), weighted with a Sudakov factor which resums the leading logarithms associated with the soft and collinear singularities. To regulate subsequent branchings, an ordering variable is defined which acts as a resolution scale [35]. As the shower runs, the ordering variable decreases and branchings occur, resolving increasingly soft and collinear radiation. At some point, the ordering variable hits a cutoff value, terminating the parton shower. Here, we first briefly describe the antenna shower formalism for leading-color gluon emissions, before continuing with photon emissions.
We denote the participating momenta as The momentum p k is the emitted gauge boson and has p 2 k = 0.

Leading-color QCD
For leading-color gluon emissions a sensible choice for the ordering variable is the transverse momentum with respect to the pre-branching system. It is defined as where the factor 4 is included to ensure t ≤ m 2 AB . This ordering variable is used in [12,13] up to a factor, and is is one of the available options in [15]. Most importantly, all singular behaviour is contained in the phase space region where t → 0. A cutoff on the scale therefore effectively regularizes the singularities, preventing situations where the emission probability diverges.
Next, we also have to factorize the (n + 1)-body phase space dΦ n+1 to completely separate the radiative piece of a cross section. To this end the antenna factorization dΦ n+1 = dΦ ant ab dΦ n (2.10) is used, where c is the three-body Gram determinant. A brief derivation of this factorization is sketched in Appendix A. We can now write down the shower approximation for a gluonic color ordered matrix element as The crucial point here is that this equality is exact in the singular regions of phase space where the shower aims to correctly reproduce the matrix element, but merely approximate elsewhere. The antenna phase space is then transformed to include the ordering variable t.
To this end, an auxilliary variable z = s ak s ak + s bk (2.13) is introduced. If all particles are massless, the boundaries of z are given by (2.14) However, for massive particles, these expressions are more complicated and are most conveniently given by the positivity of the Gram determinant. We note that the functional form of z has no physical effect and is only chosen for convenience. This particular choice provides a clear physical picture in the sense that for a given value of t, z → 0 and z → 1 correspond with the collinear phase space regions, while the soft region lies inbetween. Transforming the phase space to these variables introduces a Jacobian |J| = The shower proceeds by generating the shower variables t, z and φ and selecting a coloradjacent parton pair ab. The pre-branching system is then transformed to the postbranching system using the generated variables and the kinematics map which we detail in section 2.3. Finally, through the Sudakov veto algorithm [36][37][38], the event passes a rejection step, giving it a probability (2.16) to be accepted. Here, u is the upper boundary of the evolution variable.

QED
From the parton shower perspective the main difference between leading-color gluon radiation and photon radiation is the separation of the emission probability for different color structures. These probabilities can be partitioned in antennae as in eq. (2.16) because they correspond with different color orderings, and therefore different final states. In QED, there is no color structure and a partitioning is not possible without discarding some of the contributions to eq. (2.1). However, to maintain the parton shower picture and to allow for interleaving with QCD radiation, we would prefer to use the antenna phase space factorization of eq. (2.10), even though every charged fermion should participate in the emission. Our strategy will therefore be to divide phase space into sectors, similar to [17]. There, this approach is instead used to simplify the matching and merging procedure of the VINCIA shower.
The sector method is most easily understood by considering the QED sector shower approximation and we denote the momenta of the branching process by {p} → ({p}, k) ab if fermions a and b participate. At every phase space point, only a single term of eq. (2.17) is active. In terms of the parton shower, the photon is emitted by the charged pair it has the lowest transverse momentum with. This ensures that, as long as the 2 → 3 kinematics are infrared and collinear safe, the emission kernel never encounters a divergence. The corresponding ordering variable is This ordering variable has the required property of ensuring that all soft and collinear regions are contained in the limit t → 0, while still allowing for regular antenna shower kinematics. It can be implemented with an additional step in the Sudakov veto algorithm, which we discuss in section 2.4. It will produce emissions distributed according to We emphasize that eq. (2.17) contains the correct soft and collinear behaviour at the current scale of emission. After the shower is continued, more charged particles may be produced through photon splitting as described in section 3, or QCD radiation. In the leading logarithmic approximation, a high scale photon is blind to these lower scale effects. A similar situation appears in QCD gluon emissions, although the dynamic effects of low scale branchings are partly hidden behind the leading color approximation. Low scale branchings can still kinematically influence high scale branchings, but their effect is deemed negligible due to the ordering condition.

Kinematics
Here, we describe how the post-branching momenta are constructed from the pre-branching momenta and the shower variables. To agree with the factorization properties eq. (2.3) and eq. (2.1), this mapping must obey the following rules: Additionally, p a and p b should be treated on equal footing in the antenna picture. We use the massive generalization of the Kosower map [7] which is also used by VINCIA [15]. The pre-branching and post-branching momenta are related by The parameters x a and x b can be fixed by setting p 2 A = m 2 a and p 2 B = m 2 b , leading to where we defined Next, a choice for the functional form of the parameter r has to be made in such a way that the infrared and collinear safety conditions are satisfied. We follow VINCIA and use This choice has the property that interchanging a ↔ b corresponds with r → 1 − r, and it reduces to the massless Kosower map where The parton shower requires the inverse of the map given by eq. (2.21). Since the shower variables and m 2 AB fix the overal azimuthal angle of the system and the invariants s ak , s bk and s ab , which in turn fix the energies and relative angles between the post-branching momenta, the remaining choice for r appears only in the angle between the pre-braching and post-branching systems, at which point the kinematics are completely fixed. One might for instance compute the angle between p A and p a as which has an explicit dependence on r. The momenta are constructed in the center of mass frame of the pre-branching system using the shower variables and this angle, and are then boosted back to the original frame. For more details on the kinematic map for emissions, we refer to VINCIA [15].

Sudakov veto algorithm
The Sudakov veto algorithm is a very important component of parton shower programs [36,37,39]. It can produce samples from probability distributions like eq. (2.16) and eq. (2.20) without evaluation of the integral in the exponent or explicit inversion of the cumulative density. Here, we present the veto algorithm specialized to produce the density of eq. (2.20). We refer to [36] for a more general explanation. We first introduce some definitions to shorten the notation. As was shown in appendix A, the boundaries of the emission phase space are given by the positivity of the Gram determinant G abk . In terms of the shower variables, this tanslates to p 2 ⊥ -dependent boundaries on z. The Sudakov veto algorithm makes use of overestimates of these boundaries which are valid for any value of t. We use where t cut is the lower cutoff on the evolution variable, regulating the infrared divergences. These boundaries are based on the phase space for massless particles, which contains the massive phase space. Furthermore, we define and channel-specific weights which contain most of the phase space factors of eq. (2.15). Finally, we need an overestimate function g(t) for the emission kernel, further discussed in section 2.5, and independent and identically distributed random numbers ρ i . The Sudakov veto algorithm for photon emissions is given in algorithm 1. From the p 2 ⊥ veto step it is immediately clear that the ordering variable is t = min (p 2 ⊥ ) ab and that it is not possible for the photon to be too collinear with any of the other charged fermions. We also note that this algorithm is quite dissimilar from the standard competition algorithms mentioned in [37,39]. However, algorithm 1 is essentially a competition algorithm as well, but one that facilitates competition between sectors instead of between branching kernels. It was shown in [36] that competition can be handled in different ways, and algorithm 1 is just a particularly suitable version chosen because in this case all competing channels have the same emission kernel.
We check if this algorithm produces the density given by eq. (2.20) by writing out the probabilities step-by-step. For the sake of readability, we leave out some details such as the Compute s ak , s bk , and s ab from t, z, and m 2 dependence of the post-branching momenta on the generated shower variables.
Note that the delta function in the last line of eq. (2.30) is written rather symbolically. From the algorithmic standpoint, it just means that the newly generated momenta are accepted. Using (2.15) and taking the derivative with resect to u, we find the following differential which is solved by eq. (2.20). It is however not a unique solution. As is shown in [36], the general class of solutions include a term that corresponds with a cutoff on t. Here, we leave it out for brevity. Algorithm 1 has multiple subsequent veto steps. From the above analysis, the probabilities of these veto steps are multiplicative. In the rest of this work, some additional veto steps will be introduced, of which one is the inclusion of a running coupling. For QED, we fix the value of the coupling at the electron mass and use the ordering variable of eq. (2.19) as renormalization scale. The scale-dependent coupling is then given by where n f (t) is the number of active fermion flavors at scale t, weighted with the appropriate factors of charge and N c . In the veto algorithm, running can be incorporated by using α(u) during generation of the shower variables and vetoing with probability This factor could of course also be included in the emission kernel veto, but since that step is computationally the most expensive, some efficiency is gained by separating them and performing the cheaper veto steps first.

Determining the overestimate
Algorithm 1 makes use of an overestimate function g(t) ≥ a QED e ({p}, k). Due to the complex structure of the branching kernel of eq. (2.7) and the definition of the ordering variable eq. (2.9), the determination of this function is not easy. Since the singular structure of the branching kernel is regulated by the ordering variable, the overestimate should behave like where c is a constant. Less singular terms can be added to affect behaviour for high values of t. A simple choice for c which ensures that g(t) overestimates the branching kernel is This value can be found by discarding the contribution of the same-sign terms in eq. (2.7) and realising that the largest opposite-sign contribution can be overestimated by −16Q a Q b /t. The branching kernel is exactly equal to this overestimate for events where all same-sign fermions are collinear with each other, and anticollinear to all opposite-sign fermions. However, eq. (2.36) overestimates the branching kernel in the majority of cases by a very large margin. This issue is particularly pronounced for high multiplicities, and significantly impacts the computation time. We therefore offer an alternative that should increase performance at the cost of introducing small fluctuations in weights for the events.
To improve on eq. (2.36), properties of the pre-showering event have to be used. As it is difficult to find an upper limit of the branching kernel as a function of the pre-branching event, we resort to using a value for c such that eq. (2.35) is an incomplete overestimate. This means that it overestimates the branching kernel in the majority of phase space, but falls short in some small regions. The resulting discrepancy is corrected by introducing small deviations in the weights of the events.
We make use of a modified version of the Sudakov veto algorithm that has previously been used in [40][41][42] to estimate shower uncertainties. It is given in algorithm 2 for a general single-variable case to avoid notational cluttering. The function g(t) is an incomplete overestimate of the branching kernel f (t), which means that their ratio r(t) = f (t)/g(t) cannot used as the veto probability. Instead, a veto probability p(r(t)) is introduced and corrected for in the weights. Algorithm 2 can easily be shown to produce the desired distribution using the methods shown in the previous section and [36]. A sensible choice for should provide better performance while maintaining small fluctuations in weights. This can be achieved by choosing a function which behaves as closely as possible to r(t) for ratios between zero and one, where the overestimate g(t) truely overestimates f (t), while moving close to one as r(t) > 1. A possible choice is p(r(t)) = tanh(r(t)). (2.37) Note that the weights of the events should remain very close to unity with this choice, but there is a tiny probability to produce negatively weighted events when a scale is generated in a region where r(t) > 1, but gets rejected. Next, we search for a relation between the pre-branching event and the upper limit of the branching kernel. One variable that is correlated with this upper limit is where θ ab is the angle between the momenta p a and p b . The variable R resembles the eikonal factors without the photon energy, if they were evaluated with pre-branching fermion momenta. The correlation with the upper limit of the branching kernel is shown in figure 3 for some final states. The values c up are the approximate upper limits of the branching kernels of individual events. The incomplete overestimate where n is the number of charged fermions, is chosen to coincide with c over at r max = c over /8. In section 4.6, the performance of algorithm 2 with eq. (2.39) is tested.

Ordering
An ordering strategy is a necessary ingedient for a parton shower to produce properly resummed event samples. It serves as a means to restrict the available Monte Carlo pathways to a phase space point, such that overcounting is prevented and the leading contributions are kept. For QCD antenna showers it has been shown in [14,15] that the absence of ordering produces a substantial overcounting in the hard, wide-angle region. We perform a comparable analysis for QED radiation in section 4.1. The VINCIA parton shower implements two 1 ordering strategies that are closely connected to the method of merging the shower with exact matrix element calculations. These ordering strategies can be shown to produce the same leading-log behaviour [16]. Strong ordering is the traditional way of ordering branchings in a parton shower. An additional veto is applied in the Sudakov veto algorithm, where t prev is the scale of the previous shower branching. The application of this step function is equivalent to restarting the shower from t prev for subsequent branchings. The inclusion of this veto adds an explicitly non-Markovian element to the parton shower. In addition, unpopulated zones appear in the multi-branching phase space when there is no ordered path to reach them. In section 4.1 these zones are shown for several final state configurations. These properties of strong ordering have consequences for the implementation of the unitary matching and merging method employed by VINCIA [44]. An alternative is given by smooth ordering. Subsequent emissions are instead allowed to cover the entire available phase space. An ordering criterion is introduced by means of an additional veto in the Sudakov veto algorithm with probability wheret is the 'current' scale of the event. It is defined as the minimum of the scales that correspond with all available shower clustering of the pre-branching event. This determination of the scale accounts for all possible ways at which the shower could have reached the pre-branching state, and is therefore entirely Markovian. Smooth ordering offers the distinct advantage of covering the entire phase space with every emission. Unitary matching and merging is therefore more straightforward [14]. However, it was noted in [45] that the Sudakov factors in the unordered regions are probably not correct.

Photon splitting
From the parton shower perspective, photon splitting is much simpler than photon emission. No soft singularity is present in the photon splitting kernels, so their treatment is very similar to the QCD counterpart. As we will explain, the only significant difference is the absence of a color structure to aid in the selection of a secondary participating particle for the photon splitting. We again first define the splitting kernel and proceed with the shower implementation.

Splitting kernel
Let |M 1 ({p}, p a , p b )| 2 be the squared amplitude of some process, where p a and p b are the momenta of a charged fermion-antifermion pair. This matrix element factorizes in the quasi-collinear limit for p a and p b as defined by eq. (2.2).
where m f and Q f are the mass and charge of the fermion-antifermion pair and Since there are no additional soft properties, the factorization is unrelated to a third particle. However, an additional spectator particle is still required for the kinematics. We therefore derive the splitting kernels from matrix elements including these spectators. For a photon spectator, we make use of an effective Hγγ coupling, while for a fermionic spectator we make use of an effective µēγ coupling. As is to be expected, the splitting kernels turn out to be the same up to nonsingular terms. They are where p c is the spectator momentum.

Ordering variable and phase space
We denote the participating momenta as Since the only singularity is of collinear nature, it is sufficient to use the invariant mass of the produced fermion-antifermion pair as the ordering variable. We follow VINCIA and use the shower variables The massless boundaries on z now are (3.6) Transforming the antenna phase space to these shower variables leads to To determine the shower approximation, a method of selecting the spectator particle is required. At first glace, this choice does not seem very significant and one might be tempted to select a spectator at random. However, as will be shown in section 4.1, this can lead to a significant overcounting of the matrix element. Instead, we will generalize a method used by ARIADNE and VINCIA for gluon splitting. In leading-color QCD, the available spectators are the two partons which are color-adjacent to the splitting gluon. Labelling them I and J and the gluon K, the probability to select parton I as spectator is given by Figure 4: A shower history where a gluon splitting follows a gluon emission.
To see why, consider the shower history sketched in figure 4. Both the emission and the splitting are performed by the parton shower. This means that the gluon is on-shell after the emission. If the gluon is collinear with one of the fermions, say p I , then the invariant mass m 2 IK is small. If p I is selected as the spectator for the splitting of p K into p a and p b , the invariant mass m 2 IK remains unchanged, but if p K is selected, m 2 IK can become large. Therefore, the small value that was used in the emission kernel was incorrect by a significant margin. The selection probability of eq. (3.8) gives preference to spectators which have low invariant mass with p K , suppressing this effect.
For QED, eq. (3.8) needs to be expanded to be able to account for more than two candidate spectators. A suitable choice is where J now runs over all available spectators. It is easy to see that this definition satisfies the requirements explained above, and that it reduces to eq. (3.8) for two spectator candidates. The shower approximation is and the corresponding targeted shower probability distribution should be

Kinematics
We employ a simplified kinematic strategy for photon splitting as compared with emissions, since the map only needs to account for collinear safety. The pre-branching and postbranching momenta are related by The momenta p a and p b are multiplied with the same parameter such that the collinear limit corresponds with x → 1 and z → 0. Solving p 2 k = 0 and p 2 A = m 2 a fixes x and z. They are given by where r = (s ac + s bc ) 2 − 4m 2 c m 2 ab (3.14) Note that in the massless limit x → 1 and p C only recoils longitudinally. Then, in the collinear limit s ab → 0, z → 0 making the mapping collinear safe. Computing the angle between p C and p c , we find which fixes the kinematics entirely.

Sudakov veto algorithm
Since for photon splitting, the splitting kernel can be partitioned into independently radiating antenna functions corresponding with the choice of the photon and the specator, the Sudakov veto algorithm for photon splitting is much closer to the standard algorithms used in QCD. We again first define the overestimates for the boundaries of z and channel-specific weights The required overestimate for the antenna function is easily derived. Since eq. (3.3) is only singular in the invariant of the produced fermion-antifermion pair and the terms inside the brackets are easily overestimated with constants, the most sensible choice is where we made use of Q 2 f < 1 for all Standard Model fermions. Note that this overestimate is independent of the selected photon or spectator. The Sudakov veto algorithm for photon splitting is given by algorithm 3. As written, the algorithm only allows for a single flavor of the fermion-antifermion pair. Due to the independence of the overestimate on the selected channel, including additional flavors is straightforward. For n f flavors, the weights should be adjusted according to w KC → n f w KC (3.19) and the algorithm should be modified to select a flavor at random. Even for massive flavors, the probabilities are adjusted accordingly by the veto step.

Algorithm 3
The Sudakov veto algorithm for photon splitting t ← u loop Choose a photon K and a spectator C with probability We again analyze algorithm 3 to show that it produces the density given by eq. (3.11) by explicitly writing out the probabilities.
where again the delta function is written symbolically. The newly generated momenta are primed to distinguish them from the arguments of the probability distribution. Taking the u-derivative leads to which is solved by eq. (3.11).

Validation
In this section, we compare the QED shower method described in the previous sections to theory results. We first compare the shower approximation to fixed order matrix elements to inverstigate how it performs outside the singular regions. Next, we verify the validity of the sector strategy for photon emissions by comparing the shower implementation to a numerical solution of the DGLAP equation. Then, we compare to the YFS formalism used in [27,28] to validate the soft behaviour of the shower. The potential issue of discontinuities in the emission phase space is discussed, and the performance of the method described in section 2.5 is tested.

Comparison with matrix elements
In this section, we test the shower approximations given by eq. (2.17) and eq. (3.10). The shower approximations are only exact in their corresponding singular regions, where most radiation is produced. However, the shower populates a much larger part of the available phase space where equations (2.17) and (3.10) are only approximately valid. To gain some insight in the quality of the shower approximations we compare them to fixed-order matrix elements, varying the types and number of branchings, the type of ordering and the inclusion of the selection probability given by eq. (3.9). Similar to VINCIA [14,15], we compare to fixed order calculations by selecting some process with an n-particle final state and sampling the n-body phase space uniformly using the RAMBO algorithm [43]. The fixed-order matrix element |M n | 2 is computed using MadGraph5 [46]. The parton shower approximation is applied multiple times until an mparticle final state matching matrix element |M m | 2 is reached. To achieve this, the sampled momenta are clustered through eq. (2.21) and eq. (3.12), inverting the normal parton shower process. In most cases, the parton shower can reach a phase space point through multiple paths which all contribute to the Monte Carlo probability. In this comparison, all possible shower histories are therefore summed over. The ratio between the parton shower and the fixed order calculation is then computed as where a QED i is the branching kernel for the i-th clustering, which can be both emissions and splittings. The term P n contains additional factors including eq. (3.9) for photon splittings and an ordering factor, which is a step function in case of strong ordering or eq. (2.41) in case of smooth ordering. We keep α constant everywhere such that it drops out in eq. (4.1).

Photon emission
The sector approach to photon emissions is only different from the normal dipole or antenna shower strategy if more than two charged fermions are involved. We therefore consider decays to final states of 4 and 6 charged fermions. In the Standard Model, decays like H → 4l correspond with highly structured matrix elements which should not be probed in this comparison. Instead, we add a scalar φ to the Standard Model which directly couples to either 4 or 6 charged leptons. This causes the comparison to mostly probe the parton shower approximation only. However, we stress that the following results do not offer a direct representation of the accuracy of the shower. Not only do the results depend on the underlying process, but the parton shower does not sample phase space uniformly as is done in this comparison, but rather it prefers the singular regions where it performs best. Instead, these comparisons offer insight into the impact of algorithmic choices such as the type of ordering or the spectator selection probability given by eq. (3.9).
(PS/ME)  We first plot the ratio given by eq. (4.1) in figure 5 for a scalar φ decaying to four leptons and two photons. The matching matrix element is the decay of φ to four leptons, so the showering component only consists of two emissions. Figure 5 illustrates the impact of the types of ordering discussed in section 2.6. From the left-hand plot, it is clear that the matrix element is significantly overestimated by the shower approximation. For the process at hand, the only two available shower histories are defined by the order in which the photons are clustered. When no ordering condition is imposed, both paths contribute to the shower approximation, overestimating the matrix element. When strong ordering is imposed, one of these paths will in most cases not contribute. However, occasionally either both paths or no paths will contribute. This is caused by the changing fermion momenta after the primary clustering, changing the ordering scale of the secondary clustering. In the middle graph of figure 5, all phase space points where neither path contributes are contained in the bars on the left side. These events constitute the dead zone in the parton shower phase space. For smooth ordering, the paths are instead weighted by a continuous ordering probability, preventing the occurence of a dead zone.
In figure 6 the comparison to matrix elements is shown for one, two and three emissions from φ → 6τ . The shower approximation appears to perform better for more massive leptons, but this is caused by the decreased phase space available for the photon, pushing it to be softer more often. Finally, the graphs are somewhat shifted in the higher multiplicity Smooth ordering Figure 6: Comparison of the shower approximation to matrix elements for one, two and three emissions from φ → 6τ using strong and smooth ordering.
cases. This is a direct consequence of the absence of mass-dependent non-singular terms in the emission kernel, which become more relevant for higher lepton masses.

Photon splitting
We now incorporate photon splitting to massless and massive lepton-antilepton pairs. As discussed in section 3.2, the parton shower approximation can be improved by using a weighted selection of a spectator for a photon splitting, instead of choosing uniformly. To check this, we compare processes involving emissions and splittings from the decay Z → 2τ . In figure 7 the shower approximation is compared to matrix elements for a combination of emissions and splittings. In the top row, the spectator is selected uniformly, while in the bottom row it is selected with probability given by eq. (3.9). The weighted selection strategy improves the smoothly ordered results more than the strongly ordered result. We have checked this to be true for multiple matrix elements. In [15] it was shown that the agreement between the matrix element and the parton shower approximation depends on the combination of the choice of ordering variable and the method of selecting a spectator. Until the QED shower is interleaved with a QCD shower, the effect of the choice of ordering variable remains ambiguous and we prefer to maintain the choice made in VINCIA due to the similarity between gluon and photon splitting.
We also note that more significant shifting occurs for the higher multiplicity processes as compared with the case of pure emissions. The photon splitting kernel is singular only in the (quasi-)collinear limit, and even there only single poles occur. Photon emissions are instead associated with double poles. It is therefore to be expected that the influence of the process-specific non-singular terms increases significantly for photon splittings, worsening the quality of the parton shower approximation.  Figure 7: Comparison of the shower approximation to matrix elements for a combination of emissions and splittings from Z → 2τ . In the top row, the spectator is selected uniformly, while in the bottom row, the selection probability given by eq. (3.9) is used.
In figure 8, the impact of non-singular terms is illustrated. In the middle and right-hand plot, the splitting kernel is modified to where m Z only serves as a means to fix the dimensionality of the non-singular terms. The coefficients are loosely chosen to show that the peaks in the middle graph of figure 8 can be aligned and centralized. The same non-singular terms are used in the right-hand pane, showing that they are not a universal improvement.
Smooth ordering Weighted selection Non-singular terms Figure 8: Illustration of the impact of non-singular terms on the comparison between the parton shower approximation and matrix element calculations for two emissions and two splittings. In the left-hand graph, the Z-boson is replaced with a scalar that decays to a fermion-antifermion pair. In the middle and righthand graph, the same comparison is repeated for the decay from both a Z and a φ, but the splitting antenna function now includes non-singular terms.

Comparison with analytic resummation
For many exclusive observables, large logarithmic corrections of the form α n ln n Q 2 /m 2 or α n ln n Q 2 /t cut appear in cross section calculations at perturbative order n. These logarithms are a direct consequence of the singular regions the parton shower aims to cover correctly. Here, Q 2 is the hard squared scale of the process, and the singular behaviour leading to these logarithms is regulated by either t cut for massless particles, or m 2 for particles of mass m. Since these logarithmic contributions can be sizeable at every order in perturbation theory, they have to be resummed to all orders to obtain reliable results. Resummation can be achieved analytically or numerically using a parton shower. In this section, we compare the sector approach to photon emissions with results from analytical resummation, validating that it produces the correct collinear logarthms. For simplicity, we restrict the comparison to massless leptons.
In QCD, the evolution from high to low scales of the final state parton energy is given by the partonic fragmentation functions. While in QCD these functions are related to hadronization, a similar concept can be introduced for leptons. Naming these functions L(x, t), they describe the distribution for a lepton to retain a fraction x of its original energy at an energy squared scale t which is lower than the hard scale Q 2 . Note that these distributions are not sensitive to soft wide-angle radiation and should therefore also be reproduced by incoherent showers. The function L(x, t) is completely analogous to the QCD fragmentation function, and thus also obeys the DGLAP equation where P ll (z) is the regularized Altarelli-Parisi splitting function Usually, a transformation to Mellin space is used to turn eq. (4.3) into an ordinary linear differential equation. We instead opt to solve it numerically using the methods described in [47].
To compare with the shower approach, we start from a RAMBO-generated four-lepton system. The definition of x is where E Q 2 and E tcut are the energies of one of the leptons at respectively the hard scale and the cutoff scale. The result of the comparison is shown in figure 9. The hard scale Q 2 is set to the minimum of the invariant masses of all pairs of leptons for both the numerical DGLAP solution and the parton shower. This is the highest scale such that all sectors are able to radiate. The center of mass energy of the RAMBO event is set to 10 4 GeV. In the left and middle panel, the shower cutoff t cut is set to Λ 2 QCD ≈ 1 GeV 2 and t cut = 10 −12 GeV 2 , which is the default Pythia cutoff for photonic radiation. The coupling α is fixed to the default Pythia value at the electron mass α(m 2 e ) = 0.00729735. In the right panel, α is allowed to run from this value according to eq. (2.33). To enhance the effects of the running of α, we set n f to 35. For both the Pythia shower and the sector approach, events appear with x > 1. This is not possible in the analytical approach, but it is allowed in the shower since the 2 → 3 kinematics will sometimes raise the energy of a participating lepton. Outside this region, we observe strong agreement between both showers and the analytical approach when α is kept fixed. The agreement worsens when α is allowed to run. This is caused by the difference in scales used as argument for α in all three approaches. However, neither of the showers performs significantly better than the other. We note that, as no equivalent to the CWM scheme [48] exists for QED, there is no a priori preference for any scale, which is reflected in this result.

Effects of coherence
In the currently available parton shower approaches to photon radiation such as those of PYTHIA or PHOTOS, not all eikonal factors are included. Instead, independent dipoles are constructed such that every radiating particle is assigned a single kinematic partner, usually of opposite sign to allow for a simple probabilistic interpretation. The correct collinear limits can be achieved through the normal 1 → 2 or 2 → 3 shower schemes, as well as the eikonal factor for these particle pairs only. However, all other eikonal factors are not included. Here, we compare the sector approach with such an incoherent strategy in the antenna formalism. Since the methods are equivalent for just two final state radiators, we consider the LHCrelevant decay process H → ZZ → e − e + µ − µ + . The left-hand plot in figure 10 shows the invariant mass distribution of the leptons after application of the showers. The difference between the algorithms is minor, only appearing at the very end of the mass spectrum. The coherent branching kernel can vary from being a factor of 2 larger than the incoherent branching kernel, to completely vanishing due to destructive interference. In case of the invariant lepton mass, the differences are largely averaged out. On the right-hand side, the distribution for the angle between the electron and muon is shown for events with exactly one photon. This variable separates phase space points where the difference between the coherent and incoherent branching kernels are most pronounced.

Comparison with YFS simulation
In this section, we perform a brief comparison with the implementation of the YFS formalism as implemented by [27,28]. The YFS formalism incorporates all soft logarithms, but collinear logarithms have to be included order-by-order, similar to matrix element corrections in parton showers. The sector approach includes both the soft and the collinear logarithms without resorting to corrections. To confirm that the soft behaviour of the sector approach is consistent with the YFS method, we display the photon radiation profiles for Z → 2τ in figure 11. These radiation profiles are also shown in figure 1 in both [27] and [28] and we observe good agreement. In all cases t cut is set to 0.01 GeV 2 and strong ordering is used. In the left-hand graph, the collinear single-pole terms of eq. (2.7) are turned off, revealing their influence on the preference for hard photon production. The graphs drop off sharply at E γ = m Z /2 due to kinematic constraints. Higher values of E γ can only be reached more than one photon is emitted, which is rare in this decay.
At the particle level, the YFS method can also be used to simulate photonic radiation off W decay by treating emissions off the W as initial state radiation. This is not yet possible in our approach, and we reserve this for later work in a full electroweak parton shower. In such a shower, it makes sense to treat W and Z decay as part of the shower similar to photon splitting. If these decays are just components of the shower, the W is allowed to radiate before it eventually decays, and the decay product radiate afterwards.

Phase space discontinuities
One concern with eq. (2.7) and the sector shower approach is the presence of discontinuities in the radiative phase space on the boundaries between sectors. These discontinuities do not affect the formal accuracy of the shower since the collinear regions are far away from the boundary and in the soft region the fermion momenta are hardly changed. However, there may be an effect for high-scale photon emissions which is relevant for a potental implementation of matching and merging where the entire phase space has to be covered.
To test for artifacts of these discontinuities, we compare the shower to events generated directly according to the branching kernel. The parton shower is run from the kinematic limit on a RAMBO-generated four charged lepton event with E CM = 10 4 GeV and t cut = 1 GeV 2 . The shower is terminated after a single emission, and only the events with an emission are kept. To remove the Sudakov suppression, a CKKW-L-like [49,50] procedure is used where events are rejected with a probability that is generated using trial emissions from the scale of the actual emission. A directly generated event sample was compared with the unweighted parton shower sample, both with O(10 9 ) events, in the emission scale, the photon energy and the various leptonic invariant masses. The samples match up extremely well for all variables, giving no cause for concern for an implementation of matching or merging to matrix element calculations at a later stage.

Performance testing
We perform a brief performace comparison between the regular veto algorithm using the overestimate given by eq. (2.36) and algorithm 2 using the overestimate give by eq. (2.39). In figure 12 their relative performance for final states with an increasing number of charged leptons and a typical distribution of shower weights are shown. All events are produced with E CM = 10 4 GeV and the cutoff scale is set to t cut = 10 −6 GeV 2 . The increase in performance is substantial and the weight distribution peaks strongly at one. Negative weights can occur when a trial scale is rejected in a region where the branching kernel is larger than the overestimate. The acceptance probability eq. (2.37) is close to unity in these regions, so the probability for the appearance of negative weights is strongly suppressed. Note that the incomplete overestimate of eq.(2.39) is by no means the only possible choice. If the weights are found to fluctuate too much, the overestimate can be raised at the cost of some performance.  Figure 12: Performance comparison of the regular veto algorithm using overestimate c over and algorithm 2 using overestimate c linear . On the left, the computation time to shower 10 4 RAMBO-generated events of N electron-positron pairs is compared. On the right, the distribution of weights that results from the application of the shower using algorithm 2 on events with three electron-positron pairs is shown.

Summary and Conclusion
We described a formalism for coherent QED radiation in parton showers that is closely related to QCD antenna showers like ARIADNE and VINCIA. For photon radiation, all soft and collinear singularities are captured in a single branching kernel that is active over all of phase space. The phase space itself is divided into sectors such that branching can be regulated by an ordering parameter that remains similar to the standard antenna shower choice and the usual 2 → 3 kinematics can be used. A modified version of the Sudakov veto algorithm is presented to improve performance at the cost of introducing weighted events. For photon splitting, the methodology is much closer to the QCD analogon of gluon splitting with the exception of the presence of a color structure that can be used to dictate which spectator is used. A solution is provided by a generalized version of the so-called ARIADNE factor.
For validation, we presented several comparisons with exact matrix element calculations, the DGLAP equation and the YFS method. When performing a full phase space scan, the shower approximation shows agreement with matrix elements at a similar level as the QCD counterpart in VINCIA. Good agreement is observed with both the DGLAP equation and the YFS radiation pattern. We also compared the sector approach with an incoherent shower similar to those implemented in PYTHIA and PHOTOS for a final state with four radiators at LHC energies. The differences are currently minor, but they should become increasingly relevant at future colliders as multiple relevant radiators appear in final states, especially once initial state radiation is also included.
In the near future, the QED shower formalism described in this paper will be implemented in the VINCIA parton shower. This implementation should also include initial state radiation, which is a relatively straightforward extension of the work presented in this paper. As a consequence, all relevant interference between initial and final state radiation will be included by construction. Initial state radiation and its interference with final state radiation has already been shown to be relevant for precision measurements at the LHC [51,52] and for future colliders [53,54]. In the future, helicity dependence and an extension to a full electroweak formalism will also be included.

A Antenna phase space factorization
Here we show how phase space can be factorized as indicated by eq. (2.10). We first note that the 2-body phase space is where λ is the Källén function. We now start from the (n + 1)-body phase space where we explicitly factorize three momenta p a , p b and p c where we denoted Q = P − i =a,b,c p i . By a straightforward change of variables, this can be written as where the pre-branching particles are generally labelled with u and v, but which are related to a and b in a way dependent on the branching process.