Introduction
Understanding the effect of Quantum Electrodynamics (QED) radiation in particle decays is essential for precise measurements of particle decay properties. This is particularly relevant for the evaluation of reconstruction and selection efficiencies in experimental analyses and for the modelling of kinematic distributions. A notable example is the reconstruction of resonance structures in invariant mass spectra, where QED radiation leads to characteristic low-mass tails arising from “missing” (i.e., unreconstructed) photons. Simulations based on Monte Carlo event generators [1 , 2 , 3 ] play a central role in quantifying these effects and their associated uncertainties, and furthermore enable the embedding of decays into fully realistic collider-event environments. For QED radiation in hadron decays, state-of-the-art generators currently include the Sherpa’s Photons++ module [4 , 5 , 6 ] , the Herwig 7 generator [7 , 8 , 9 ] , and the Photos package [10 , 11 , 12 , 13 , 14 ] . The latter is commonly used in combination with the Pythia [15 , 16 ] and/or Evtgen generator [17 , 18 , 19 ] .
QED radiation generators account for the four-momentum “lost” to QED radiation by modifying the four-momenta of the decay products in decay chains that are initially generated without any final-state radiation (i.e., taking each complete decay chain, including all its particles, as input). The four-momenta are adjusted in such a way that the QED factorisation properties are preserved in the infrared regime, specifically in the soft and/or collinear limits. In practice, this typically requires adopting a formalism based on the soft-photon approximation, as realised in the YFS approach [20 ] employed by Photons++ and Herwig, and implemented as a correction weight in Photos, or on the collinear approximation, as described by DGLAP evolution [21 , 22 , 23 ] , which is used in Pythia and Photos and can be incorporated through higher-order corrections in Photons++.
An alternative and comparatively less explored framework for QED radiation is the antenna formalism [24 , 25 , 26 , 27 , 28 , 29 , 30 , 31 ] . This formalism combines exact phase-space factorisations [27 ] with (multipole) QED antenna functions [31 ] that have YFS-style eikonal factors as their soft limits and DGLAP splitting kernels as their (quasi)collinear limits. From a theoretical standpoint, this offers a unified and accurate description of both soft and collinear limits, without requiring a prioritisation of one over the other. It provides the basis for the treatment of QED radiation in the Vincia parton-shower module [32 ] of Pythia 8 [16 ] , but so far has only been developed for QED radiation from elementary particles such as quarks, leptons, and W W W
In this paper, we present the initial developments towards a new framework for the simulation of QED radiation in generic hadron decays within the Vincia shower module. The algorithm is based on the sector-antenna approximation [28 , 30 , 31 , 32 , 33 , 34 ] , implemented with full mass dependence [32 , 35 , 36 ] , and extends beyond the YFS formalism by incorporating spin-dependent (quasi)collinear single-unresolved limits to all orders. In its present form, the framework includes the relevant spin-dependent contributions for photon emission from spin-0, spin-1/2, and spin-1 particles. In addition, secondary γ → ℓ + ℓ − \gamma\to\ell^{+}\ell^{-} γ → ℓ + ℓ − γ → q q ˉ \gamma\to q\bar{q} γ → q q ˉ 1
An additional important development is the inclusion of a framework that incorporates matrix-element (ME) corrections [33 , 37 , 38 , 39 , 40 , 41 ] directly into the QED parton-shower evolution via an interface to ME libraries. It can be created semi-automatically with FeynRules [42 , 43 ] and MadGraph [44 , 45 ] , enabling modern effective field theory (EFT) and ME-matching methods for QED radiative corrections in hadronic decays. We demonstrate the tool chain (Figure 1 ) for QED corrections involving pseudoscalar- and vector-meson decays into fermion pairs. Extensions to additional classes of decays (for example, further vector-meson and baryonic decay modes) will be the subject of future studies.
Figure 1. Illustration of the tool chain from an effective field theory (EFT) + QED for the hadron sector to matrix elements that can be interfaced from Pythia 8.
We compare the obtained results with those from other well-established algorithms for a representative set of benchmark decay channels. For these comparisons, we simulate decays using the Evtgen generator, which is widely employed in numerous experiments to model the decays of unstable heavy hadrons. Evtgen relies on specialised external packages to incorporate final-state QED radiative corrections into the simulated decays, using Photos in its default configuration. Recently, the Evtgen generator has been upgraded to support both Photons++ and Vincia as alternative options to Photos, enabling comparisons between QED radiation generators within a single framework.
The outline of this paper is as follows. In sec. 2 we discuss the underlying theory approaches of common general-purpose QED-correction generators. In sec. 3 we detail the Vincia QED shower formalism. In sec. 4 we discuss the dead-cone effect in the centre of mass frame. In sec. 5 we look at the effects of changing the infrared cutoff between generators. In sec. 6 we look in more detail at the size of the ME corrections. In sec. 7 we discuss results for benchmark decay modes and compare the results obtained with Vincia with the results obtained using other generators. In sec. 8 we compare the performance of the different generators in both space and time. Finally, in sec. 9, we provide a summary and an outlook.
Approaches for QED radiation in current generators
In this section, we summarise the methodologies currently implemented in the event generators that are used as benchmarks for our comparison with the Vincia simulations. In general, QED radiation generators exploit the factorisation properties of QED splitting kernels such that the differential cross section factorises into the Born-level (leading-order) cross section multiplied by appropriate QED splitting kernel(s) and — beyond tree level — Sudakov factors.
The contributions from real-photon emission are, by themselves, infrared divergent. However, these divergences are cancelled order by order in perturbation theory by virtual corrections of opposite sign. According to the KLN theorem [46 , 47 ] , the combined sum of the virtual and bremsstrahlung contributions is free of infrared divergences. In practical implementations,
Fixed-order generators numerically reconstruct virtual corrections from the real-photon contributions following the slicing approach [48 ] , which partitions the phase space into a resolved and unresolved region based on a cutoff parameter. In the unresolved region, the real-emission matrix element factorises and can be analytically integrated up to the cutoff by exploiting its universal infrared limits. The resolved region, being finite, is computed numerically.
Shower generators, on the other hand, implement exact unitarity by construction. Their probability-conserving Sudakov factors embody the KLN theorem explicitly,
Sing [ V i r t u a l ] = − Sing [ ∫ R e a l ] , (1) \operatorname{Sing}[\mathrm{Virtual}] = -\operatorname{Sing}\left[\int\mathrm{Real}\right], \tag*{(1)} Sing [ Virtual ] = − Sing [ ∫ Real ] , (1) where the notation Sing[X] means the singular part of X (typically cast as explicit poles in a regularisation parameter), order by order in the infinite-order sum generated by the shower. The infrared cutoff of the QED shower (which is typically user-definable, with a low default value ≪ 1 \ll1 ≪ 1
Note that different applications may choose different functional forms for the slicing or cutoff variable. E.g., in one algorithm it may be a photon energy (in the decay rest frame) or energy fraction, while in another it may be a transverse-momentum type variable. This is important to keep in mind when comparing the numerical sizes of slicing or cutoff values to each other, and will be important when comparing results between different QED generators below.
Standard QED factorisation strategies effectively reduce to two main classes depending on the baseline approximation for the QED emission kernels: those that choose the soft–emission limit as their baseline, as formulated in the YFS framework, and subsequently incorporate corrections for the collinear and hard regimes, and those that instead start from the collinear limit, as formulated in the DGLAP framework, and then introduce corrections to account for the soft and hard regimes.
Algorithms for QED corrections also differ in whether they operate as showers or as “N-generators”. Showers use Sudakov-style evolution equations to sweep over phase space in an ordered manner, generating branchings iteratively as they go along, with the first emission always being the hardest (according to the given shower’s evolution variable) and the last being the softest (closest to the infrared cutoff of the shower). Vetoes may be applied during one or more of these steps to implement corrections to the baseline shower kernels. N-generators begin from a prescribed target distribution (e.g., a Poissonian distribution with a specific mean) for the total number of QED branchings, sample a value of N N N N N N
Tab. 1 gives an overview of the type (N-generator or shower) and starting limit (collinear, soft, or both) of each of the current Monte Carlo approaches to QED radiation in hadron decays.
PHOTOS
As indicated in tab. 1, the PHOTOS generator [10 , 11 , 12 , 13 , 14 ] follows an approach based on the collinear approximation: the QED splitting kernel of a charged decay product into itself and a photon, X → X γ X \to X\gamma X → X γ 1 / 2 1/2 1/2
splitting kernel incorporates the sum of Born-level, virtual, and real photon contributions, where the latter are integrated up to a cutoff energy. To control systematic errors, the algorithm covers the full phase space without approximations, imposes four-momentum conservation, and makes use of the explicit analytical form of the approximate one-photon matrix element.
If the exact matrix element for a certain channel is known, Photos can correct for possible deviations with the help of a weight, corresponding to the ratio between the exact matrix element and the approximation. Otherwise, Photos employs generalised weighting factors to account for the soft-photon emission regime, as well as for discrepancies arising from the DGLAP evolution kernels associated with potentially different spins of the decay products. In the case of a (pseudo)scalar meson decaying into two (pseudo)scalar particles, Photos has an implemented calculation of the NLO matrix element [13 ] that must be activated by the user to be effective (we will refer to this configuration as Photos MEC in § 7 ). In the case of decays with multiple charged particles, Photos applies an additional correction weight to account for interference effects, ensuring proper infrared behaviour.
PHOTONS++
The Photons++ module [4 ] (henceforth referred to as Photons), which is embedded in the Sherpa generator [5 , 6 ] , is based on the soft-photon approximation following the YFS formalism [20 ] . The algorithm separates and exponentiates the universal, process-independent components of the QED corrections, while the remaining process-dependent and spin-dependent contributions are incorporated via subtraction terms.
The algorithm supports decays into multi-particle final states and implements a dedicated momentum-mapping procedure to relate the original momenta prior to QED radiation to those after radiation. This is required because the eikonal factors entering the YFS splitting kernel are formulated in terms of the already radiation-corrected final-state momenta, rather than the original pre-radiation kinematics. In this mapping, all final-state momenta are treated symmetrically: their three-momenta are uniformly rescaled by a common factor, the photon momenta are distributed accordingly, and the energy components of all four-vectors are assigned such that overall momentum conservation and all on-shell constraints are simultaneously satisfied.
As an optional feature, the Photons generator can account for the collinear regime by employing dipole splitting functions [49 , 50 ] , generalised to consistently include polarisation degrees of freedom. These dipole splitting functions are used as universal multiplicative correction factors, ensuring the correct reproduction of the collinear limit. Since they smoothly interpolate between the soft and collinear regions, they inherently contain the soft limit. In constructing the correction weights, the soft contributions are explicitly subtracted, as these are already resummed in the original YFS eikonal factors.
To achieve higher precision, the Photons generator incorporates exact first-order, full matrix elements. Substantial deviations from approximate treatments of the matrix elements occur in nonsingular regions characterised by comparably hard, wide-angle emissions. Accordingly, the algorithm employs infrared-subtracted one-photon real-emission matrix elements, which are subsequently evaluated at the generated momentum configurations. Exact first-order matrix elements obtained in this manner are available for a number of relevant processes, such as the decay of a neutral vector boson into a pair of charged fermions.
In the comparative analyses presented in § 7 , the uncorrected configuration employing the YFS splitting kernels will be denoted by Photons 1. The configuration applying approximate corrections derived from dipole splitting functions will be denoted by Photons 20, while the configuration incorporating exact next-to-leading-order (NLO) ME corrections will be indicated by Photons 21.
PYTHIA
The default parton shower in the Pythia generator [16 ] implements a fully interleaved treatment of final-state QCD+QED radiation [51 , 52 ] . The shower evolution is ordered in decreasing transverse-momentum (p ⊥ p_\perp p ⊥ § 3 ), but differs in that its starting point is the collinear approximation, using DGLAP splitting kernels for f → f γ f \to f\gamma f → f γ γ ( ∗ ) → f f ˉ \gamma^{(*)} \to f\bar{f} γ ( ∗ ) → f f ˉ f f f p ⊥ p_\perp p ⊥ X R → x r γ XR \to xr\gamma X R → x r γ X X X R R R p ⊥ p_\perp p ⊥ [51 ]
p ⊥ E 2 = z ( 1 − z ) Q 2 , (2) p_{\perp\mathrm{E}}^2 = z(1-z)Q^2, \tag*{(2)} p ⊥ E 2 = z ( 1 − z ) Q 2 , (2) where the subscript E stands for “evolution”, z = z x ∈ [ 0 , 1 ] z = z_x \in[0,1] z = z x ∈ [ 0 , 1 ]
Q 2 = s x γ ≡ ( p x + p γ ) 2 − m x 2 Q^2 = s_{x\gamma} \equiv(p_x+p_\gamma)^2-m_x^2 Q 2 = s x γ ≡ ( p x + p γ ) 2 − m x 2 represents the virtuality (offshellness) of the branching. In the expression above, p x p_x p x p γ p_\gamma p γ s s s 2
s i j ≡ 2 p i ⋅ p j = ( p i + p j ) 2 − m i 2 − m j 2 , (3) s_{ij} \equiv2p_i \cdot p_j = (p_i+p_j)^2-m_i^2-m_j^2, \tag*{(3)} s ij ≡ 2 p i ⋅ p j = ( p i + p j ) 2 − m i 2 − m j 2 , (3) for particles i i i j j j p i p_i p i p j p_j p j m i m_i m i m j m_j m j [25 , 53 ]
p ⊥ A 2 = s x γ s r γ s x γ + s r γ + s x r → γ ∥ x ( 1 − z ) Q 2 , (4) p_{\perp\mathrm{A}}^2 = \frac{s_{x\gamma}s_{r\gamma}}{s_{x\gamma}+s_{r\gamma}+s_{xr}} \xrightarrow{\gamma\parallel x} (1-z)Q^2, \tag*{(4)} p ⊥ A 2 = s x γ + s r γ + s x r s x γ s r γ γ ∥ x ( 1 − z ) Q 2 , (4) where the subscript A stands for both Ariadne (which used the same evolution variable [54 ] ) and antenna. The rightmost expression corresponds to the γ ∥ x \gamma\parallel x γ ∥ x
It is readily observed that eq. (2) and the γ ∥ x \gamma\parallel x γ ∥ x (4) both vanish in the same way for z → 1 z \to1 z → 1 Q 2 → 0 Q^2 \to0 Q 2 → 0 x γ → 1 ⇒ z → 0 x_\gamma\to1 \Rightarrow z \to0 x γ → 1 ⇒ z → 0 p ⊥ E 2 p_{\perp\mathrm{E}}^2 p ⊥ E 2 Q 2 Q^2 Q 2 z z z z z z p ⊥ A 2 p_{\perp\mathrm{A}}^2 p ⊥ A 2
The shower evolution in Pythia is performed by identifying radiation dipoles associated with the decay process. In configurations with three or more charged particles, a principle of “maximal screening” is used to determine an optimal set of such dipoles. This does not capture the full multipole radiation pattern but should provide a reasonable approximation to it for the typical case of higher multipole moments being suppressed relative to the leading dipole contributions.
For simple ( 2 → 1 ) (2 \to1) ( 2 → 1 ) 1 → 2 1 \to2 1 → 2 V 0 → ℓ + ℓ − V^0 \to\ell^+\ell^- V 0 → ℓ + ℓ − V + → ℓ + ν V^+ \to\ell^+\nu V + → ℓ + ν Z → ℓ + ℓ − Z \to\ell^+\ell^- Z → ℓ + ℓ − W + → ℓ + ν W^+ \to\ell^+\nu W + → ℓ + ν 1 / 2 1/2 1/2
The Pythia parton shower can be configured to evolve exclusively via its QED component. In sec. 7, we use this restricted parton shower to obtain an additional reference sample to study QED radiation in the decay of a vector meson into a pair of charged leptons, for which Pythia’s description should be accurate and includes the NLO ME correction. We note that no form factor for the decaying meson is included in this description, but such form factors are not included in any of the other approaches discussed here.
The VINCIA antenna-shower approach
The Vincia generator [34 ] implements algorithms for QED and QCD parton showers with an interleaved [51 , 55 ] , p ⊥ p_\perp p ⊥ [32 , 53 ] based on sector antennae [31 , 32 , 33 , 34 ] .
In the present study, we restrict our attention to Vincia’s QED antenna showers which are detailed in [31 , 32 ] . Conceptually, these are closely analogous to their QCD counterparts [34 ] , except of course that the photon is uncharged, electric charges are not confined, and there is no “colour flow” or 1 / N C 1/N_C 1/ N C 1 / N C 1/N_C 1/ N C [31 , 32 ] , although a simpler (and somewhat faster) option analogous to Pythia’s maximal-screening principle is also available. Comparisons between these options are included in some of the validation plots presented in sec. 7.
The antenna formalism provides a unified description of both the soft and collinear divergent behaviour of radiative matrix elements. In contrast to shower formalisms based on Catani-Seymour (CS) dipoles [56 , 57 ] , in which the antenna function is effectively partitioned into two separately radiating “dipole ends”, the antenna framework treats the radiation pattern as arising coherently from the full dipole pair3
Compared to approaches based on YFS exponentiation, such as Photos+++, or to algorithms based on collinear logarithms that are included only order by order, such as Photos, and which are used as “afterburners” to generate coherent soft-photon radiation in particle decays, the antenna parton-shower approach not only reproduces coherent QED radiation but can also be systematically interleaved with QCD radiation and with resonance decays [55 ] . In the present work, however, we restrict ourselves to using Vincia as an afterburner and defer a detailed study of its interleaving with hadronic decays to future investigations.
The default QED antenna functions implemented in VINCIA [31 , 32 ] incorporate fully coherent multipole soft-interference contributions in addition to the standard collinear DGLAP structures. Mass effects are systematically taken into account [34 , 35 ] both through mass-dependent corrections to the antenna functions, ensuring that all relevant (quasi-)collinear limits are correctly reproduced, and through the application of exact phase-space factorisations for massive particles. Hadronic form factors are neglected, as are possible effects arising from photon mixing with hadronic vector resonances (for example, through vector meson dominance) and effective photon splitting into pairs of charged particles. Treatment of these aspects is deferred to future work.
The QED multipole shower formalism in VINCIA
The default QED parton shower in VINCIA employs a fully coherent multipole formalism for the simulation of photon radiation. In its original implementation, the algorithm was designed to describe radiation from systems of elementary charged fermions, vector bosons, and scalar particles, as well as photon splitting into pairs of charged fermions [31 , 32 ] .
The principal extension introduced in the present work is a generalisation of the antenna functions to allow for the description of radiation emitted from charged hadrons, incorporating an explicit dependence on the hadron spin. In the current implementation, this includes dedicated antenna functions for spin-0, spin-1/2, and spin-1 hadrons. Hadrons with higher integer spins, if provided, are effectively treated as scalar particles, whereas higher-spin fermions are treated the same as spin-1/2 particles.
To introduce the underlying antenna formalism, we consider I K → i k γ IK \to ik\gamma I K → ik γ 4 I I I K K K p I p_I p I p K p_K p K p i p_i p i p k p_k p k p j p_j p j
I ( p I ) + K ( p K ) ↦ i ( p i ) + k ( p k ) + γ ( p j ) . I(p_I) + K(p_K) \mapsto i(p_i) + k(p_k) + \gamma(p_j). I ( p I ) + K ( p K ) ↦ i ( p i ) + k ( p k ) + γ ( p j ) . If both I I I K K K I I I 5
The emission probability densities are described by (unpolarised) antenna functions, whose leading (eikonal) terms are universal, i.e., independent of the total spins of the emitting particles and the same for both FF and IF antennae. We therefore express the photon-emission dipole antenna functions as
a e m i t i k ( s i j , s j k , s i k , m i 2 , m k 2 ) = a s c a l a r i k + a s p i n i + a s p i n k , a_{\mathrm{emit}}^{ik}(s_{ij},s_{jk},s_{ik},m_i^2,m_k^2) = a_{\mathrm{scalar}}^{ik} + a_{\mathrm{spin}}^i + a_{\mathrm{spin}}^k, a emit ik ( s ij , s j k , s ik , m i 2 , m k 2 ) = a scalar ik + a spin i + a spin k , where the spin-independent eikonal term (presented here without the corresponding electromagnetic charge and coupling factors) is:
a s c a l a r i k = 4 s i k s i j s j k − 4 m i 2 s i j 2 − 4 m k 2 s j k 2 , (5) a_{\mathrm{scalar}}^{ik} = 4\frac{s_{ik}}{s_{ij}s_{jk}} - 4\frac{m_i^2}{s_{ij}^2} - 4\frac{m_k^2}{s_{jk}^2}, \tag*{(5)} a scalar ik = 4 s ij s j k s ik − 4 s ij 2 m i 2 − 4 s j k 2 m k 2 , (5) with the Lorentz-invariant quantity s x y ≡ 2 p x μ p y μ s_{xy} \equiv2p_x^\mu p_{y\mu} s x y ≡ 2 p x μ p y μ (3) , and rest masses m i = m I m_i = m_I m i = m I m k = m K m_k = m_K m k = m K
The hard-collinear terms for spin-1/2 emitters are also independent of whether the antenna is FF or IF:
a s p i n − 1 2 i = 2 s I K s j k s i j (6) a_{\mathrm{spin-\frac{1}{2}}}^{i} = \frac{2}{s_{IK}}\frac{s_{jk}}{s_{ij}} \tag*{(6)} a spin − 2 1 i = s I K 2 s ij s j k (6) a s p i n − 1 2 k = 2 s I K s i j s j k (7) a_{\mathrm{spin-\frac{1}{2}}}^{k} = \frac{2}{s_{IK}}\frac{s_{ij}}{s_{jk}} \tag*{(7)} a spin − 2 1 k = s I K 2 s j k s ij (7) while for spin-1 emitters, the signs of some of the collinear invariants depend on whether the antenna is FF or IF (related by crossing),
a s p i n − 1 F F , i = 8 3 1 s i j ( s j k s I K + s j k s I K − s j k − s j k 2 s I K 2 ) (8) a_{\mathrm{spin-1}}^{\mathrm{FF},i} = \frac{8}{3}\frac{1}{s_{ij}}\left(\frac{s_{jk}}{s_{IK}}+\frac{s_{jk}}{s_{IK}-s_{jk}}-\frac{s_{jk}^{2}}{s_{IK}^{2}}\right) \tag*{(8)} a spin − 1 FF , i = 3 8 s ij 1 ( s I K s j k + s I K − s j k s j k − s I K 2 s j k 2 ) (8) a s p i n − 1 F F , k = 8 3 1 s j k ( s i j s I K + s i j s I K − s i j − s i j 2 s I K 2 ) (9) a_{\mathrm{spin-1}}^{\mathrm{FF},k} = \frac{8}{3}\frac{1}{s_{jk}}\left(\frac{s_{ij}}{s_{IK}}+\frac{s_{ij}}{s_{IK}-s_{ij}}-\frac{s_{ij}^{2}}{s_{IK}^{2}}\right) \tag*{(9)} a spin − 1 FF , k = 3 8 s j k 1 ( s I K s ij + s I K − s ij s ij − s I K 2 s ij 2 ) (9) a s p i n − 1 I F , i = 8 3 1 s i j ( s j k s I K + s j k s I K + s j k + s j k 2 s I K 2 ) (10) a_{\mathrm{spin-1}}^{\mathrm{IF},i} = \frac{8}{3}\frac{1}{s_{ij}}\left(\frac{s_{jk}}{s_{IK}}+\frac{s_{jk}}{s_{IK}+s_{jk}}+\frac{s_{jk}^{2}}{s_{IK}^{2}}\right) \tag*{(10)} a spin − 1 IF , i = 3 8 s ij 1 ( s I K s j k + s I K + s j k s j k + s I K 2 s j k 2 ) (10) a s p i n − 1 I F , k = 8 3 1 s j k ( s i j s i k + s k j + s i j s I K + s j k − s i j 2 ( s I K + s j k ) 2 ) (11) a_{\mathrm{spin-1}}^{\mathrm{IF},k} = \frac{8}{3}\frac{1}{s_{jk}}\left(\frac{s_{ij}}{s_{ik}+s_{kj}}+\frac{s_{ij}}{s_{IK}+s_{jk}}-\frac{s_{ij}^{2}}{(s_{IK}+s_{jk})^{2}}\right) \tag*{(11)} a spin − 1 IF , k = 3 8 s j k 1 ( s ik + s k j s ij + s I K + s j k s ij − ( s I K + s j k ) 2 s ij 2 ) (11) The full dipole antenna function, for a pair of charged particles i i i k k k
a e m i t i k ( s i j , s j k , s i k , m i 2 , m k 2 ) = a s c a l a r i k + δ s p i n − 1 2 i a s p i n − 1 2 i + δ s p i n − 1 2 k a s p i n − 1 2 k + δ s p i n − 1 i a s p i n − 1 i + δ s p i n − 1 k a s p i n − 1 k (12) a_{\mathrm{emit}}^{ik}(s_{ij},s_{jk},s_{ik},m_i^{2},m_k^{2}) = a_{\mathrm{scalar}}^{ik} + \delta_{\mathrm{spin-\frac{1}{2}}}^{i}a_{\mathrm{spin-\frac{1}{2}}}^{i} + \delta_{\mathrm{spin-\frac{1}{2}}}^{k}a_{\mathrm{spin-\frac{1}{2}}}^{k} + \delta_{\mathrm{spin-1}}^{i}a_{\mathrm{spin-1}}^{i} + \delta_{\mathrm{spin-1}}^{k}a_{\mathrm{spin-1}}^{k} \tag*{(12)} a emit ik ( s ij , s j k , s ik , m i 2 , m k 2 ) = a scalar ik + δ spin − 2 1 i a spin − 2 1 i + δ spin − 2 1 k a spin − 2 1 k + δ spin − 1 i a spin − 1 i + δ spin − 1 k a spin − 1 k (12) where δ s p i n − S x = 1 \delta_{\mathrm{spin-S}}^{x}=1 δ spin − S x = 1 x x x S S S
The fully coherent (multipole) antenna function, in which all eikonal contributions are incorporated simultaneously within a single splitting kernel, is given by
a e m i t m u l t i p o l e ( { p } , p j ) = − ∑ { i < k } σ i Q i σ k Q k a e m i t i k ( s i j , s j k , s i k , m i 2 , m k 2 ) (13) a_{\mathrm{emit}}^{\mathrm{multipole}}(\{p\},p_j)=-\sum_{\{i<k\}}\sigma_i Q_i\sigma_k Q_k a_{\mathrm{emit}}^{ik}(s_{ij},s_{jk},s_{ik},m_i^{2},m_k^{2}) \tag*{(13)} a emit multipole ({ p } , p j ) = − { i < k } ∑ σ i Q i σ k Q k a emit ik ( s ij , s j k , s ik , m i 2 , m k 2 ) (13) where the sum over { i < k } \{i<k\} { i < k } 6 Q i Q_i Q i Q k Q_k Q k σ i \sigma_i σ i σ k \sigma_k σ k + 1 +1 + 1 − 1 -1 − 1
This splitting kernel contains all soft multipole contributions and also reproduces the correct collinear limits [31 ] , although its singularity structure is quite intricate. In the Vincia sector shower, double counting of singular configurations is avoided by allowing only a single such splitting kernel to contribute at any given phase-space point. This is achieved by partitioning the phase space into sectors, following the formalism of ref. [33 ] , with the modification that the antenna function itself is employed as the sector-defining variable.
Consequently, the shower approximation to the NLO matrix element can be written as
∣ M n + 1 P S ( { p } , p j ) ∣ 2 = a e m i t m u l t i p o l e ( { p } , p j ) 4 π α ∑ { i < k } Θ j ( a e m i t i k ) ∣ M n ( { p } ) ∣ 2 , (14) \left|M^{\mathrm{PS}}_{n+1}(\{p\},p_j)\right|^2 = a^{\mathrm{multipole}}_{\mathrm{emit}}(\{p\},p_j)\,4\pi\alpha\sum_{\{i<k\}}\Theta_j\left(a^{ik}_{\mathrm{emit}}\right)\left|M_n(\{p\})\right|^2, \tag*{(14)} M n + 1 PS ({ p } , p j ) 2 = a emit multipole ({ p } , p j ) 4 π α { i < k } ∑ Θ j ( a emit ik ) ∣ M n ({ p }) ∣ 2 , (14) where the superscript PS stands for “Parton Shower”, α \alpha α Θ \Theta Θ
Θ j = { 1 if a e m i t i k ( s i j , s j k , s i k , m i 2 , m k 2 ) ≥ a e m i t v w ( s v j , s j w , s v w , m v 2 , m w 2 ) ∀ pairs { v , w } , 0 otherwise , \Theta_j =
\begin{cases}
1 & \text{if } a^{ik}_{\mathrm{emit}}(s_{ij},s_{jk},s_{ik},m_i^2,m_k^2) \ge a^{vw}_{\mathrm{emit}}(s_{vj},s_{jw},s_{vw},m_v^2,m_w^2)\ \forall\ \text{pairs }\{v,w\},\\
0 & \text{otherwise},
\end{cases} Θ j = { 1 0 if a emit ik ( s ij , s j k , s ik , m i 2 , m k 2 ) ≥ a emit v w ( s v j , s j w , s v w , m v 2 , m w 2 ) ∀ pairs { v , w } , otherwise , which effectively selects the sector corresponding to the pair of charged particles i i i k k k { p } \{p\} { p } j j j ( { p } , p j ) (\{p\},p_j) ({ p } , p j ) { i < k } \{i<k\} { i < k } (14) is ever non-zero. A technical description of the implementation of the VINCIA parton shower, as well as its interface within PYTHIA, is provided in app. A.
The simple dipole-shower option
The VINCIA generator also provides a simpler and computationally more efficient alternative in which the complete multipole sum is replaced by a set of individual dipole contributions, constructed according to a principle of maximal screening analogous to that used in PYTHIA’s QED shower (see § 2.3 ).
In this formulation, radiation is generated solely from a collection of dipoles, with each charged particle only appearing in exactly one dipole, and the corresponding dipole shower adopts a faster, unsectorised scheme. The antenna function in eq. (13) is no longer evaluated as a sum over all charged pairs, but instead only over a specific set of pairs which each satisfy σ i Q i σ k Q k < 0 \sigma_i Q_i \sigma_k Q_k < 0 σ i Q i σ k Q k < 0
Following the principle of maximal screening [32 ] , the selection of pairs { i , k } \{i,k\} { i , k } ∑ { i < k } [ ( p i + p k ) 2 − ( m i + m k ) 2 ] \sum_{\{i<k\}} [(p_i+p_k)^2-(m_i+m_k)^2] ∑ { i < k } [( p i + p k ) 2 − ( m i + m k ) 2 ] O ( n 3 ) O(n^3) O ( n 3 ) [58 , 59 , 60 ] .
Although this algorithmic strategy is generally more efficient, it only generates a (summed) dipole approximation to the full multipole structure. In particular, for configurations involving three or more charged particles, the exact soft-photon multipole radiation pattern is not fully captured.
Matrix-element corrections
In the presence of hard, wide-angle photons, the nonsingular contributions in the matrix-element factorisation become phenomenologically significant. These are not guaranteed to be captured by shower kernels (irrespective of whether they are based on DGLAP, dipoles/antennae, or YFS). To incorporate these effects, the antenna function in eq. (13) can be replaced by the exact tree-level matrix-element ratio, up to a given (fixed) tree-level order. In most current descriptions, such ME corrections are only implemented to first order. We have implemented corrections for up to two-photon emissions in VINCIA by default, but it is fairly straightforward to generate and interface corrections for more photons if so desired. Instructions for how to generate and implement further ones will be given in the PYTHIA’s HTML manual [61 ] .
VINCIA currently supports generic QED ME corrections for two-body decays of a scalar into two leptons or two scalars, and for three-body decays of a scalar into all-scalar final states.
We plan to extend these to processes involving generic fermions and spin-1 hadrons in future work. The ME corrections are implemented by employing a combination of FeynRules [42 , 43 ] and MadGraph [44 , 45 ] . Details of the implementation and a summary of the ME corrections included with Vincia in the next release of Pythia (tentatively, Pythia 8.319) are given in app. B.
Currently, matrix elements are evaluated only for one photon at a time. That is, the QED shower proceeds by iterated single emissions, which can include iterated ME corrections but does not incorporate so-called direct multiple-emission contributions such as the direct 2 → 4 2 \to4 2 → 4 [62 , 63 ] . Thus, at each step, the ME-corrected antenna corresponds to the ratio
a e m i t Q E D ( { p } , p j ) = ∣ M n + 1 ( { p } , p j ) ∣ 2 ∣ M n ( { p } / p j ) ∣ 2 , a_{\mathrm{emit}}^{\mathrm{QED}}(\{p\},p_j) = \frac{\left|M_{n+1}(\{p\},p_j)\right|^2}{\left|M_n(\{p\}/p_j)\right|^2}, a emit QED ({ p } , p j ) = ∣ M n ({ p } / p j ) ∣ 2 ∣ M n + 1 ({ p } , p j ) ∣ 2 , where the notation { p } / p j \{p\}/p_j { p } / p j p j p_j p j { p } \{p\} { p }
We note that ME corrections for branchings involving γ → ℓ + ℓ − \gamma\to\ell^{+}\ell^{-} γ → ℓ + ℓ − ℓ + ℓ − \ell^{+}\ell^{-} ℓ + ℓ −
The treatment of virtual corrections proceeds analogously to their implementation in standard parton-shower algorithms, namely through the enforcement of unitarity via Sudakov form factors. In parton showers, one exploits the fact that, at the level of singular contributions, the virtual correction is exactly equal in magnitude and opposite in sign to the phase-space integral of the corresponding real-emission correction, cf. eq. (1) . The nonsingular remainder is neglected, which is consistent with the approximation that the correction being implemented is dominated by the infrared and collinear singular structures.
When employing tree-level matrix elements for real emissions instead of their antenna approximations, the integral over real emissions is modified by a nonsingular correction. Since the corresponding virtual corrections are still incorporated via eq. (1) , this also induces a change in the effective virtual correction. The formal perturbative accuracy of the calculation is neither improved nor degraded by this. To include explicit virtual corrections in a controlled manner, one would have to perform a loop-level matching to fixed-order calculations, such as in MC@NLO [64 ] , Powheg [65 ] , or equivalent strategies.
At the level of inclusive decay rates, it is sufficient to normalise the total decay width to the known NLO prediction (or to experimentally measured total widths) in order to effectively incorporate the NLO virtual corrections. Virtual photon contributions induce relative corrections of O ( α ) \mathcal{O}(\alpha) O ( α )
The dead-cone effect
The dead-cone effect is a well-established phenomenon in both classical and quantum gauge field theories, manifesting itself as a strong suppression of gauge-boson emissions within a polar angle of ∼ m / E = 1 / γ \sim m/E = 1/\gamma ∼ m / E = 1/ γ m m m E E E [66 , 67 , 68 ] for recent demonstrations of this effect for QCD radiation off charm and bottom quarks.
Briefly stated, the differential photon emission probability for bremsstrahlung off a massless charged particle at small opening angles θ \theta θ
d P ∝ d E E d θ γ 2 θ γ 2 . (15) \mathrm{d}P \propto\frac{\mathrm{d}E}{E}\frac{\mathrm{d}\theta_\gamma^2}{\theta_\gamma^2}. \tag*{(15)} d P ∝ E d E θ γ 2 d θ γ 2 . (15) Plotted logarithmically in photon energy and opening angle, the emission density is roughly constant (and proportional to α \alpha α massive particle, the emission probability has the same d θ γ 2 / θ γ 2 \mathrm{d}\theta_\gamma^2/\theta_\gamma^2 d θ γ 2 / θ γ 2 (15) for angles above θ d e a d ∼ 1 / γ = m / E \theta_{\mathrm{dead}} \sim1/\gamma= m/E θ dead ∼ 1/ γ = m / E θ γ 2 d θ γ 2 \theta_\gamma^2\mathrm{d}\theta_\gamma^2 θ γ 2 d θ γ 2 θ γ → 0 \theta_\gamma\to0 θ γ → 0 Figure 2 , where we show the local effective exponent of the soft-limit antenna function defined as p e f f ( θ γ ) = d ln ( a n t ) / d ln θ γ p_{\mathrm{eff}}(\theta_\gamma) = \mathrm{d}\ln(\mathrm{ant})/\mathrm{d}\ln\theta_\gamma p eff ( θ γ ) = d ln ( ant ) / d ln θ γ γ = 2 \gamma= 2 γ = 2 10 10 10 50 50 50 θ = 1 / γ \theta= 1/\gamma θ = 1/ γ γ \gamma γ θ = cos − 1 ( β ) \theta= \cos^{-1}(\beta) θ = cos − 1 ( β )
Figure 2. The (left) local effective exponent and (right) massive-to-massless antenna ratio for soft photon emissions, as functions of emission angle (note the logarithmic axes), for a few representative Lorentz factors γ = 2 \gamma= 2 γ = 2 10 10 10 50 50 50 θ = 1 / γ \theta= 1/\gamma θ = 1/ γ γ \gamma γ θ = cos − 1 ( β ) \theta= \cos^{-1}(\beta) θ = cos − 1 ( β )
Notation and basics
Consider a QED dipole-antenna spanned between a pair of oppositely charged scalar particles (see eq. (5) ). For simplicity, consider the case with only m I = m i ≠ 0 m_I = m_i \ne0 m I = m i = 0 m K = m k = 0 m_K = m_k = 0 m K = m k = 0 p j = p γ p_j = p_\gamma p j = p γ [69 ]
a ( p i , p γ , p k ) ∝ s i k s i γ s k γ − m i 2 s i γ 2 . a(p_i,p_\gamma,p_k) \propto\frac{s_{ik}}{s_{i\gamma}s_{k\gamma}}-\frac{m_i^2}{s_{i\gamma}^2}. a ( p i , p γ , p k ) ∝ s iγ s k γ s ik − s iγ 2 m i 2 . Requiring the antenna function to be non-negative yields
s i γ s i k − m i 2 s k γ ≥ 0. (16) s_{i\gamma}s_{ik}-m_i^2s_{k\gamma} \ge0. \tag*{(16)} s iγ s ik − m i 2 s k γ ≥ 0. (16) This is exactly equal to the expression one obtains from requiring a positive Gram determinant, which defines the accessible phase space for massive kinematics. That is, eq. (16) represents both the physical boundary of the massive phase space as well as the contour along which the (negative) mass correction to the scalar antenna function becomes equal to the (positive) massless eikonal term.
Parametric rapidity and dead-cone angle
Before discussing the dead cone in terms of angular variables, we first note that the constraint in eq. (16) implies hard boundaries on an invariant quantity, which we call the parametric rapidity of the photon,
y γ ∗ ≡ 1 2 ln ( s K γ s i γ ) ≤ 1 2 ln ( m I K 2 ( 1 − x γ ) − m I 2 m I 2 ) (17) y_\gamma^* \equiv\frac{1}{2}\ln\left(\frac{s_{K\gamma}}{s_{i\gamma}}\right) \le\frac{1}{2}\ln\left(\frac{m_{IK}^2(1-x_\gamma)-m_I^2}{m_I^2}\right) \tag*{(17)} y γ ∗ ≡ 2 1 ln ( s iγ s K γ ) ≤ 2 1 ln ( m I 2 m I K 2 ( 1 − x γ ) − m I 2 ) (17) where m I K m_{IK} m I K I K IK I K I I I K K K x γ = 2 E γ / m I K x_\gamma= 2E_\gamma/m_{IK} x γ = 2 E γ / m I K m K ≠ 0 m_K \ne0 m K = 0
y γ ∗ ≤ m K ≠ 0 1 2 ln ( m I K 2 ( 1 − x γ ) − m I 2 − m K 2 m I 2 ) + 1 2 ln ( 1 2 + 1 2 1 − 4 m I 2 m K 2 ( m I K 2 ( 1 − x γ ) − m I 2 − m K 2 ) 2 ) (18) y_\gamma^* \mathrel{\mathop{\leq}^{m_K\ne0}} \frac{1}{2}\ln\left(\frac{m_{IK}^2(1-x_\gamma)-m_I^2-m_K^2}{m_I^2}\right) + \frac{1}{2}\ln\left(\frac{1}{2}+\frac{1}{2}\sqrt{1-\frac{4m_I^2m_K^2}{\left(m_{IK}^2(1-x_\gamma)-m_I^2-m_K^2\right)^2}}\right) \tag*{(18)} y γ ∗ ≤ m K = 0 2 1 ln ( m I 2 m I K 2 ( 1 − x γ ) − m I 2 − m K 2 ) + 2 1 ln ( 2 1 + 2 1 1 − ( m I K 2 ( 1 − x γ ) − m I 2 − m K 2 ) 2 4 m I 2 m K 2 ) (18) We note that this bound is always smaller than that of eq. (17) . In the soft-photon (x γ → 0 x_\gamma\to0 x γ → 0 y γ ∗ y_\gamma^* y γ ∗
− ln ( m I K 2 − m K 2 m K 2 ) + m I 2 2 m I K 2 ≤ y γ ∗ ≤ 1 2 ln ( m I K 2 − m I 2 m I 2 ) − m K 2 2 m I K 2 (19) -\ln\left(\frac{m_{IK}^2-m_K^2}{m_K^2}\right)+\frac{m_I^2}{2m_{IK}^2} \le y_\gamma^* \le\frac{1}{2}\ln\left(\frac{m_{IK}^2-m_I^2}{m_I^2}\right)-\frac{m_K^2}{2m_{IK}^2} \tag*{(19)} − ln ( m K 2 m I K 2 − m K 2 ) + 2 m I K 2 m I 2 ≤ y γ ∗ ≤ 2 1 ln ( m I 2 m I K 2 − m I 2 ) − 2 m I K 2 m K 2 (19) where terms quartic in m / m I K m/m_{IK} m / m I K
The quantity y γ ∗ y_\gamma^* y γ ∗
y γ ∗ = 1 2 ln ( s K γ s i γ ) = 1 2 ln ( E k − ∣ p ⃗ k ∣ cos θ k γ E i − ∣ p ⃗ i ∣ cos θ i γ ) = 1 2 ln ( E k E i ) + 1 2 ln ( 1 − β k cos θ k γ 1 − β i cos θ i γ ) (20) y_\gamma^*=\frac{1}{2}\ln\left(\frac{s_{K\gamma}}{s_{i\gamma}}\right)=\frac{1}{2}\ln\left(\frac{E_k-\lvert\vec{p}_k\rvert\cos\theta_{k\gamma}}{E_i-\lvert\vec{p}_i\rvert\cos\theta_{i\gamma}}\right)=\frac{1}{2}\ln\left(\frac{E_k}{E_i}\right)+\frac{1}{2}\ln\left(\frac{1-\beta_k\cos\theta_{k\gamma}}{1-\beta_i\cos\theta_{i\gamma}}\right) \tag*{(20)} y γ ∗ = 2 1 ln ( s iγ s K γ ) = 2 1 ln ( E i − ∣ p i ∣ cos θ iγ E k − ∣ p k ∣ cos θ k γ ) = 2 1 ln ( E i E k ) + 2 1 ln ( 1 − β i cos θ iγ 1 − β k cos θ k γ ) (20) with E i , k = ∣ p ⃗ i , k ∣ 2 + m i , k 2 E_{i,k}=\sqrt{\lvert\vec{p}_{i,k}\rvert^2+m_{i,k}^2} E i , k = ∣ p i , k ∣ 2 + m i , k 2
y γ ∗ ⟶ x γ → 0 1 2 ln ( m I K 2 + m K 2 − m I 2 m I K 2 + m I 2 − m K 2 ) + 1 2 ln ( 1 − β K cos θ k γ 1 − β I cos θ i γ ) y_\gamma^* \mathrel{\mathop{\longrightarrow}^{x_\gamma\to0}} \frac{1}{2}\ln\left(\frac{m_{IK}^2+m_K^2-m_I^2}{m_{IK}^2+m_I^2-m_K^2}\right)+\frac{1}{2}\ln\left(\frac{1-\beta_K\cos\theta_{k\gamma}}{1-\beta_I\cos\theta_{i\gamma}}\right) y γ ∗ ⟶ x γ → 0 2 1 ln ( m I K 2 + m I 2 − m K 2 m I K 2 + m K 2 − m I 2 ) + 2 1 ln ( 1 − β I cos θ iγ 1 − β K cos θ k γ ) where, in the CM-frame, cos θ k γ ∼ − cos θ i γ \cos\theta_{k\gamma}\sim-\cos\theta_{i\gamma} cos θ k γ ∼ − cos θ iγ y γ ∗ y_\gamma^* y γ ∗
Moving on to the conventional dead-cone angle for soft radiation around a massive emitter, we note that this does not correspond to an actual zero of the massive antenna function but to the ratio of the massive to massless antenna functions becoming small. To derive it, we start from this ratio,
ant ( m I ≠ 0 ) ant ( m I = 0 ) = ( 1 − cos θ i j ) ( 1 − β i cos θ i k ) ( 1 − β i cos θ i j ) ( 1 − cos θ i k ) − ( 1 − cos θ i j ) ( 1 − cos θ j k ) 2 ( 1 − cos θ i k ) ( 1 − β i cos θ i j ) γ i 2 \frac{\operatorname{ant}(m_I\ne0)}{\operatorname{ant}(m_I=0)}=\frac{(1-\cos\theta_{ij})(1-\beta_i\cos\theta_{ik})}{(1-\beta_i\cos\theta_{ij})(1-\cos\theta_{ik})}-\frac{(1-\cos\theta_{ij})(1-\cos\theta_{jk})}{2(1-\cos\theta_{ik})(1-\beta_i\cos\theta_{ij})}\gamma_i^2 ant ( m I = 0 ) ant ( m I = 0 ) = ( 1 − β i cos θ ij ) ( 1 − cos θ ik ) ( 1 − cos θ ij ) ( 1 − β i cos θ ik ) − 2 ( 1 − cos θ ik ) ( 1 − β i cos θ ij ) ( 1 − cos θ ij ) ( 1 − cos θ j k ) γ i 2 where we have not yet made any approximations apart from setting m K = 0 m_K=0 m K = 0 γ i \gamma_i γ i i i i
In the soft-photon limit, we adopt the no-recoil approximation in the center-of-mass (CM) frame, i.e. cos θ i k → − 1 \cos\theta_{ik}\to-1 cos θ ik → − 1 cos θ j k → − cos θ i j \cos\theta_{jk}\to-\cos\theta_{ij} cos θ j k → − cos θ ij 1 / 4 1/4 1/4
cos ( θ d e a d ) = β i . (21) \cos(\theta_{\mathrm{dead}})=\beta_i. \tag*{(21)} cos ( θ dead ) = β i . (21) Performing an expansion about the relativistic limit β i → 1 \beta_i\to1 β i → 1
θ d e a d = 1 γ i + 1 6 γ i 3 + O ( γ i − 5 ) + O ( γ i − 3 γ k − 2 ) , (22) \theta_{\mathrm{dead}}=\frac{1}{\gamma_i}+\frac{1}{6\gamma_i^3}+O(\gamma_i^{-5})+O(\gamma_i^{-3}\gamma_k^{-2}), \tag*{(22)} θ dead = γ i 1 + 6 γ i 3 1 + O ( γ i − 5 ) + O ( γ i − 3 γ k − 2 ) , (22) where the first term corresponds to the standard dead-cone angle definition, γ i − 1 = m i / E i \gamma_i^{-1}=m_i/E_i γ i − 1 = m i / E i m k m_k m k X → 1 + 2 X\to1+2 X → 1 + 2 (22) or the unexpanded form in eq. (21) .
An alternative definition of the dead-cone angle can be derived by identifying the stationary point of the massive antenna function. At this point, the nominal divergence 1 / θ 2 1/\theta^2 1/ θ 2 1 / 4 1/4 1/4 X 0 → 1 + 2 − X^0\to1^+2^- X 0 → 1 + 2 − ( m X , m 1 , m 2 ) = ( 20 , 1 , 2 ) G e V (m_X,m_1,m_2)=(20,1,2)\,\mathrm{GeV} ( m X , m 1 , m 2 ) = ( 20 , 1 , 2 ) GeV 10 − 8 G e V 10^{-8}\,\mathrm{GeV} 1 0 − 8 GeV cos − 1 ( β ) \cos^{-1}(\beta) cos − 1 ( β ) d I / d θ = 0 d\mathcal{I}/d\theta=0 d I / d θ = 0 1 / θ 2 1/\theta^2 1/ θ 2 θ ∈ [ 0 ∘ , 180 ∘ ] \theta\in[0^\circ,180^\circ] θ ∈ [ 0 ∘ , 18 0 ∘ ]
However, when one or both emitters are non-relativistic, the stationary-point definition becomes more complicated and requires a partial-fractioning of the dipole antenna pattern to separate it into two pieces, one for each emitter. An illustration of such a case is shown in fig. 3b corresponding to the B 0 → D + π − B^0\to D^+\pi^- B 0 → D + π − D D D 1 / 4 1/4 1/4 = 1 / 4 =1/4 = 1/4
Figure 3. Angular radiation patterns for massive antenna functions. The dead-cone angles cos − 1 ( β ) \cos^{-1}(\beta) cos − 1 ( β ) D D D
We now examine the full angular dependence of the ratio of the massive to the massless antenna functions in the soft limit (x γ → 0 x_\gamma\to0 x γ → 0 X 0 → 1 + 2 − X^0 \to1^+2^- X 0 → 1 + 2 − Figure 4 , we show polar plots with this ratio on the radial axis (red lines), as a function of the polar emission angle of the photon, relative to the decay axis in the CM frame. Where the red line touches the outermost circle, the ratio is unity, i.e., massive = massless. Where the red line intersects the dashed black circle, the ratio is 1 / 2 1/2 1/2 1 / 4 1/4 1/4 cos − 1 ( β ∗ ) \cos^{-1}(\beta^*) cos − 1 ( β ∗ ) Table 2 . We consider the following four representative scenarios:
Table 2. CM-frame velocities β ∗ \beta^* β ∗ γ ∗ \gamma^* γ ∗ 1 / γ ∗ 1/\gamma^* 1/ γ ∗ cos − 1 β ∗ \cos^{-1}\beta^* cos − 1 β ∗ [2 ] .
Figure 4. Massive-to-massless antenna ratio for soft photon emission vs. photon polar angle relative to the CM-frame decay axis. The ratio is the radius of the red curve. Dead-cone angles, cos − 1 ( β ∗ ) \cos^{-1}(\beta^*) cos − 1 ( β ∗ ) Table 2 are shown as dark grey regions. Circles mark ratios of 1 (outermost), 1/2 (dashed), and 1/4 (solid).
a) J / ψ → μ + μ − J/\psi\to\mu^+\mu^- J / ψ → μ + μ − ∼ ± 4 ∘ \sim\pm4^\circ ∼ ± 4 ∘ Table 2 ;
b) B 0 → D s + D − B^0 \to D_s^+D^- B 0 → D s + D − 1 / 2 1/2 1/2
c) D 0 → K + K − D^0 \to K^+K^- D 0 → K + K − ∼ 30 ∘ \sim30^\circ ∼ 3 0 ∘ Table 2 );
d) D 0 → K − π + D^0 \to K^-\pi^+ D 0 → K − π +
Phase-space distributions for finite photon energies
As emphasised above, the parametric rapidity, y γ ∗ y_\gamma^* y γ ∗ (20) , is kinematically bounded for all photon energies (see eq. (19) ). In contrast, for opening angles θ γ \theta_{\gamma} θ γ η γ \eta_\gamma η γ 1 / θ 2 1/\theta^2 1/ θ 2
This is illustrated explicitly in Figure 5 and Figure 6 . Both figures present the distribution of the first photon emission for a representative decay, D 0 → K − π + D^0 \to K^-\pi^+ D 0 → K − π + p ⊥ 2 p_\perp^2 p ⊥ 2 y γ ∗ y_\gamma^* y γ ∗ Figure 5 and versus ln ( 1 / θ ) \ln(1/\theta) ln ( 1/ θ ) Figure 6 (horizontal axes).
Figure 5. Distribution of p ⊥ 2 p_{\perp}^{2} p ⊥ 2 y ∗ y^{*} y ∗ D 0 → K − π + D^{0}\to K^{-}\pi^{+} D 0 → K − π + y ∗ y^{*} y ∗ θ = 0 \theta=0 θ = 0 y ∗ y^{*} y ∗ θ = cos − 1 ( β ) \theta=\cos^{-1}(\beta) θ = cos − 1 ( β ) θ \theta θ y ∗ y^{*} y ∗
(a) p ⊥ 2 p_{\perp}^{2} p ⊥ 2 ln ( 1 / θ π γ ) \ln(1/\theta_{\pi\gamma}) ln ( 1/ θ π γ )
(b) p ⊥ 2 p_{\perp}^{2} p ⊥ 2 ln ( 1 / θ K γ ) \ln(1/\theta_{K\gamma}) ln ( 1/ θ K γ )
Figure 6. Distributions of p ⊥ 2 p_{\perp}^{2} p ⊥ 2 ln ( 1 / θ ) \ln(1/\theta) ln ( 1/ θ ) θ \theta θ θ = cos − 1 ( β ) \theta= \cos^{-1}(\beta) θ = cos − 1 ( β ) θ ≥ π 2 \theta\ge\frac{\pi}{2} θ ≥ 2 π
In Figure 5 , a pure massless eikonal would give a constant density in these phase-space variables. The thin green line marks where the Gram determinant vanishes, i.e., the physical phase-space boundary for full kinematics. The pink dashed lines show the soft-photon-limit asymptotes of y γ ∗ y_\gamma^* y γ ∗ y ∗ y^* y ∗
In contrast, Figure 6 shows the radiation pattern around a single emitter. The pion (Figure 6a ) radiates at angles smaller than those of the kaon (Figure 6b ) for a given p ⊥ p_\perp p ⊥ ∼ 2 \sim2 ∼ 2 θ → π \theta\to\pi θ → π
Using log ( 1 / θ ) \log(1/\theta) log ( 1/ θ ) θ ≥ π / 2 \theta\ge\pi/2 θ ≥ π /2 Figure 6 ). To show both limits (θ → 0 \theta\to0 θ → 0 θ → π \theta\to\pi θ → π y ∗ y^{*} y ∗
η ^ = 1 2 log ( 1 + cos θ ^ 1 − cos θ ^ ) or η = 1 2 log ( 1 + cos θ 1 − cos θ ) , \hat{\eta}=\frac{1}{2}\log\left(\frac{1+\cos\hat{\theta}}{1-\cos\hat{\theta}}\right) \qquad\text{or} \qquad\eta=\frac{1}{2}\log\left(\frac{1+\cos\theta}{1-\cos\theta}\right), η ^ = 2 1 log ( 1 − cos θ ^ 1 + cos θ ^ ) or η = 2 1 log ( 1 − cos θ 1 + cos θ ) , where θ ^ \hat{\theta} θ ^ θ \theta θ Figure 7 . A comparison of Figure 7a and Figure 7b shows that the distributions for low p ⊥ p_{\perp} p ⊥ p ⊥ p_{\perp} p ⊥
(a) p ⊥ 2 p_{\perp}^{2} p ⊥ 2 η ^ \hat{\eta} η ^
(b) p ⊥ 2 p_{\perp}^{2} p ⊥ 2 η \eta η
Figure 7. Distributions of p ⊥ 2 p_{\perp}^{2} p ⊥ 2 η ^ \hat{\eta} η ^ η \eta η D 0 → K − π + D^{0}\to K^{-}\pi^{+} D 0 → K − π +
Recoil effects have received much attention lately in the context of QCD showers [70 ] , where so-called local recoil schemes (in which the recoil from each emission is distributed among only two partons, which are often referred to as the radiator and the recoiler) have been shown to yield only leading-logarithmic (LL) accuracy in the soft-collinear limit, while global recoil schemes (in which the recoil is distributed among multiple partons) can reach higher accuracies. Among the QED generators considered here, Pythia employs a local dipole-recoil scheme, in which the emitter absorbs the full transverse recoil while the recoiler recoils only longitudinally in the rest frame of the emitting dipole, Vincia uses a local antenna-recoil scheme in which both parents share the transverse recoil (though for resonance-final emissions its recoil strategy does have global properties [36 ] ), and Sherpa and Photos both use variants of global recoil schemes. These differences in recoil strategies imply that the generators will typically differ more from each other in the θ ^ \hat{\theta} θ ^
axis) than in the post-branching opening angle θ \theta θ θ \theta θ θ ^ \hat{\theta} θ ^ § 7 .
To illustrate the differences between casting angular distributions in terms of θ \theta θ θ ^ \hat{\theta} θ ^ Figure 8 and Figure 9 show the angular distributions of the total emitted photon energy in J / ψ → e + e − J/\psi\to e^+e^- J / ψ → e + e − Figure 8 ) and D 0 → K − π + D^0 \to K^-\pi^+ D 0 → K − π + Figure 9 ) for various QED generators7 J / ψ → e + e − J/\psi\to e^+e^- J / ψ → e + e − θ \theta θ θ \theta θ Figure 10 shows that this agreement persists at small angles around the size of the dead cone (using the standard m i / E i m_i/E_i m i / E i
Figure 8. (a) The x-axis is the opening angle between the photon and the post-branching final-state electron. (b) The x-axis is the opening angle between the photon and the pre-branching born-level electron. Angular distribution of photon energy in J / ψ → e + e − J/\psi\to e^+e^- J / ψ → e + e −
Figure 9. (a) The x-axis is the opening angle between the photon and the post-branching final-state kaon. (b) The x-axis is the opening angle between the photon and the pre-branching born-level kaon. Angular distribution of photon energy in B 0 → K + π − B^0 \to K^+\pi^- B 0 → K + π −
Figure 10. (a) J / ψ → e + e − J/\psi\to e^+e^- J / ψ → e + e − B 0 → π + π − B^0 \to\pi^+\pi^- B 0 → π + π − θ i γ / θ d e a d = θ i γ ( E i / m i ) \theta_{i\gamma}/\theta_{\mathrm{dead}} = \theta_{i\gamma}(E_i/m_i) θ iγ / θ dead = θ iγ ( E i / m i ) J / ψ → e + e − J/\psi\to e^+e^- J / ψ → e + e − B 0 → π + π − B^0 \to\pi^+\pi^- B 0 → π + π − θ i γ \theta_{i\gamma} θ iγ
Infrared cutoff dependence
Both parton showers and N N N Table 1 ) impose a lower (“infrared”) cutoff on the phase space of photon emissions. To illustrate why this is needed and how it matters in practice, consider the leading-order photon-emission probability density for a highly relativistic particle as a function of photon energy and emission angle, which has the familiar double-logarithmic form of eq. (15) in the soft small-angle limit, with singularities for E γ → 0 E_\gamma\to0 E γ → 0 θ γ → 0 \theta_\gamma\to0 θ γ → 0
Table 1. Type and starting limit of each of the current Monte Carlo approaches to QED radiation in hadron decays. The default starting limit of the VINCIA event generator is based on the multipole antenna formalism. This configuration can be converted to a soft starting limit based on multipole eikonals (indicated with grey parenthesis) by setting the flag Vincia::useSpinQEDHadDec={off,off,off}, where the three entries are for spin-1/2, spin-1, and spin-3/2 particles respectively.
The concept of a “bare” charged particle is therefore not quantum mechanically meaningful: only the concept of a charged particle plus some amount of unresolved (very soft/collinear) radiation is. We emphasise that this does not, a priori , present a problem for experimental measurements. However, it does become important when comparing experiment to theory. It cannot be emphasised enough that the number of photons that are emitted, in a given theoretical calculation, completely depends on the infrared cutoff it uses. It would not be meaningful at all to compare this number between two calculations that do not use the same infrared cutoffs. This is related not only to the numerical value of the cutoff but also to the specific kinematic variable that defines it, which varies from generator to generator. What can be meaningfully compared is the number of photons that are resolved with respect to some common definition of what a resolved photon means, e.g., in terms of photon energy, transverse momentum, and/or emission angle. This definition of a resolved photon must be infrared safe : it must be insensitive to soft- and collinear photon emissions. For example, the number of photons with energies above some minimum value and opening angles above some minimum emission angle would be well defined. The cumulative energy carried away by the photons above the minimum emission angle is likewise meaningful.
Considering first the soft singularity (E → 0 E \to0 E → 0
The collinear singularity can be more subtle and, especially for highly relativistic emitters like electrons8 θ ∼ m / E \theta\sim m/E θ ∼ m / E d θ 2 / θ 2 \mathrm{d}\theta^2/\theta^2 d θ 2 / θ 2 θ 2 d θ 2 \theta^2\mathrm{d}\theta^2 θ 2 d θ 2 § 4 . However, for highly relativistic emitters, the d θ 2 / θ 2 \mathrm{d}\theta^2/\theta^2 d θ 2 / θ 2 can be resolvable experimentally, significantly affecting observable quantities such as the cumulative energy loss of the emitter.
Experimentally, what happens if an electron emits a hard photon at almost zero emission angle? In the absence of magnetic fields in the detector, the electron and the photon would propagate together and hit the same calorimeter cell, which would only record their combined energy (disregarding possible information from the cluster shape). In such a setup, the effective experimental cutoff on photon emission angle would simply correspond to the size of the calorimeter cells. To reproduce this theoretically, the infrared cutoff must go deep enough into the collinear region to resolve all detectable photons emitted at angles larger than the angular resolution of the calorimeter, while photons at smaller angles should be summed back with their parent emitters, mimicking the sum done by the calorimeter. In practice, however, the magnetic fields in the detector deflect the emitting electron away from any photons it was originally collinear with. To make the experimental electron energy insensitive to photon emission, one must extrapolate the charged-particle trajectory back to its original direction and sum its energy together with that of any reconstructed photons collinear with it. This defines a “dressed” electron: the electron track plus any calorimeter clusters consistent with the chosen definition of a photon that should be added to the electron.9
For hadron decays, the default infrared cutoffs in event generators are typically small enough that they effectively satisfy the conditions above. A “worst-case” scenario is the J / ψ → e + e − J/\psi\to e^{+}e^{-} J / ψ → e + e − ∼ 0.02 ∘ \sim0.02^\circ ∼ 0.0 2 ∘
In Table 3 , we compare the four QED generators used in this study in terms of their cutoff variables and default values and provide the cutoff values used in this study:
Table 3. Infrared cutoff values used for the QED generators under comparison, along with their defaults. Since Photos uses the energy fraction as the cutoff variable, we divide the chosen energy cutoff by M / 2 M/2 M /2 M M M
z γ = 2 E γ M m o t h e r , z_{\gamma} = \frac{2E_{\gamma}}{M_{\mathrm{mother}}}, z γ = M mother 2 E γ , with a default value of z γ min = 10 − 7 z_{\gamma}^{\min}=10^{-7} z γ m i n = 1 0 − 7 z γ min = 2 × 10 − 5 / M m o t h e r z_{\gamma}^{\min}=2\times10^{-5}/M_{\mathrm{mother}} z γ m i n = 2 × 1 0 − 5 / M mother E γ min = 10 − 5 G e V = 10 k e V E_{\gamma}^{\min}=10^{-5}\,\mathrm{GeV}=10\,\mathrm{keV} E γ m i n = 1 0 − 5 GeV = 10 keV
In Sherpa’s Photons module, the cutoff corresponds to the minimum photon energy with default value E γ min = 10 − 3 G e V = 1 M e V E_{\gamma}^{\min}=10^{-3}\,\mathrm{GeV}=1\,\mathrm{MeV} E γ m i n = 1 0 − 3 GeV = 1 MeV 10 11 E γ min = 10 − 5 G e V = 10 k e V E_{\gamma}^{\min}=10^{-5}\,\mathrm{GeV}=10\,\mathrm{keV} E γ m i n = 1 0 − 5 GeV = 10 keV
cutoff in Photons, so we assume it self-regulates the collinear divergence using particle masses.
In Vincia, the infrared cutoff is applied to a Lorentz-invariant transverse-momentum measure, rather than energies or angles. This is the Ariadne p ⊥ p_\perp p ⊥ [25 , 34 , 53 , 54 ] introduced in sec. 2.3. In terms of the photon energy fraction z γ z_\gamma z γ y γ ∗ y_\gamma^* y γ ∗ (17) ),
p ⊥ Λ = z γ e y γ ∗ + e − y γ ∗ M m o t h e r . p_{\perp\Lambda} = \frac{z_\gamma}{e^{y_\gamma^*} + e^{-y_\gamma^*}} M_{\mathrm{mother}} . p ⊥ Λ = e y γ ∗ + e − y γ ∗ z γ M mother . For relativistic emitters in the soft/collinear limit, p ⊥ p_\perp p ⊥
p ⊥ γ → E γ θ γ , p_{\perp\gamma} \to E_\gamma\theta_\gamma, p ⊥ γ → E γ θ γ , so the cutoff effectively constrains a combination of energy and angle and is not directly comparable to an energy cutoff. Since at fixed E γ E_\gamma E γ p ⊥ γ p_{\perp\gamma} p ⊥ γ p ⊥ γ min = 10 − 6 G e V = 1 k e V p_{\perp\gamma}^{\min} = 10^{-6}\,\mathrm{GeV} = 1\,\mathrm{keV} p ⊥ γ m i n = 1 0 − 6 GeV = 1 keV
In Pythia, QED parton showers are cut off in p ⊥ E p_{\perp E} p ⊥ E [51 ] (see eq. (2) ), which has the same soft-collinear limit as the Vincia evolution variable (see discussion in sec. 2.3). We therefore use p ⊥ e v o l min = 10 − 6 G e V = 1 k e V p_{\perp\mathrm{evol}}^{\min} = 10^{-6}\,\mathrm{GeV} = 1\,\mathrm{keV} p ⊥ evol m i n = 1 0 − 6 GeV = 1 keV
We study how the spectrum of total emitted photon energy, ∑ E γ \sum E_\gamma ∑ E γ J / ψ → e + e − J/\psi\to e^+e^- J / ψ → e + e − Figure 11 ). In each case, the cutoff is varied from 100 MeV (solid yellow) down to 100 eV (dotted purple), with ratios shown relative to 10 keV. The photon-energy sum corresponds closely to individual photon energies; deviations arise only from multi-photon emission, which is suppressed by an extra power of α \alpha α Figure 11 b). Only for the largest cutoff, E γ min = 100 M e V E_\gamma^{\min} = 100\,\mathrm{MeV} E γ m i n = 100 MeV ∑ E γ = 100 M e V \sum E_\gamma= 100\,\mathrm{MeV} ∑ E γ = 100 MeV ∑ E γ \sum E_\gamma ∑ E γ ∑ E γ > 1 2 m J / ψ ≈ 1.55 G e V \sum E_\gamma> \frac{1}{2}m_{J/\psi} \approx1.55\,\mathrm{GeV} ∑ E γ > 2 1 m J / ψ ≈ 1.55 GeV ≤ 1 M e V \le1\,\mathrm{MeV} ≤ 1 MeV
Figure 11. Total photon energy radiated in J / ψ → e + e − J/\psi\to e^{+}e^{-} J / ψ → e + e −
We expected a similar behaviour in Photos, but Figure 11 a shows a more complex cutoff dependence than Sherpa. For cutoff values between 10 MeV and 10 keV, the results are stable within a few percent, motivating our default choice of 10 − 5 G e V = 10 k e V 10^{-5}\,\mathrm{GeV} = 10\,\mathrm{keV} 1 0 − 5 GeV = 10 keV
In Vincia and Pythia, the cutoff is defined in terms of p ⊥ p_\perp p ⊥ Figure 11 c and Figure 11 d. If we instead plotted a p ⊥ p_\perp p ⊥
Figure 12 compares p ⊥ p_\perp p ⊥ p ⊥ p_\perp p ⊥ B 0 → K + π − B^0 \to K^+\pi^- B 0 → K + π − Figure 12 a), setting the Vincia p ⊥ p_\perp p ⊥ J / ψ → e + e − J/\psi\to e^+e^- J / ψ → e + e − Figure 12 b), lowering
Figure 12. Base-10 logarithm of the total emitted photon energy in (left) B 0 → K + π − B^0 \to K^+\pi^- B 0 → K + π − J / ψ → e + e − J/\psi\to e^+e^- J / ψ → e + e − x x x y y y
the VINCIA cutoff has a similar effect, bringing it in line with the other generators for another order of magnitude. This further motivates our chosen cutoff values for VINCIA. However, in the extreme case of J / ψ → e + e − J/\psi\to e^+e^- J / ψ → e + e − Figure 12 a it follows the generators down to their energy cutoffs, but then continues well below its expected cutoff. This likely stems from PYTHIA’s handling of emitters whose masses are large relative to the available phase space, producing kinks in the spectral tails when only the kaon radiates. In J / ψ → e + e − J/\psi\to e^+e^- J / ψ → e + e − log 10 ( ∑ E γ ) ∼ − 3.5 \log_{10}(\sum E_\gamma) \sim-3.5 log 10 ( ∑ E γ ) ∼ − 3.5 ∼ − 5.5 \sim-5.5 ∼ − 5.5
In many real-world physics applications, one could likely use cutoffs about ten times larger than those used in these studies, unless resolving very soft and/or collinear photons is required. We instead use low cutoffs to ensure that we resolve all contributions to the observables studied.
Quantifying matrix-element corrections
The first-order real-emission QED corrections for vector-to-dilepton decays such as J / ψ → μ + μ − J/\psi\to\mu^+\mu^- J / ψ → μ + μ − B s 0 → μ + μ − B_s^0 \to\mu^+\mu^- B s 0 → μ + μ − [35 ] :
∣ M ( X 0 → f f ˉ γ ) ∣ 2 ∣ M ( X 0 → f f ˉ ) ∣ 2 = 8 π α Q f 2 M X 2 ( 2 y f f ˉ y f γ y f f ˉ γ − 2 μ f 2 y f γ 2 − 2 μ f 2 y f γ ˉ 2 + y f γ ˉ y f γ + y f γ y f γ ˉ + c N S X ) \frac{\left|M(X^0 \to f\bar{f}\gamma)\right|^2}{\left|M(X^0 \to f\bar{f})\right|^2}
= \frac{8\pi\alpha Q_f^2}{M_X^2}
\left(
\frac{2y_{f\bar{f}}}{y_{f\gamma}y_{f\bar{f}\gamma}}
-\frac{2\mu_f^2}{y_{f\gamma}^2}
-\frac{2\mu_f^2}{y_{f\bar{\gamma}}^2}
+\frac{y_{f\bar{\gamma}}}{y_{f\gamma}}
+\frac{y_{f\gamma}}{y_{f\bar{\gamma}}}
+c_{\mathrm{NS}}^X
\right) M ( X 0 → f f ˉ ) 2 M ( X 0 → f f ˉ γ ) 2 = M X 2 8 π α Q f 2 ( y f γ y f f ˉ γ 2 y f f ˉ − y f γ 2 2 μ f 2 − y f γ ˉ 2 2 μ f 2 + y f γ y f γ ˉ + y f γ ˉ y f γ + c NS X ) where f f f m f 2 m_f^2 m f 2 μ f 2 ≡ m f 2 / m X 2 \mu_f^2 \equiv m_f^2/m_X^2 μ f 2 ≡ m f 2 / m X 2 y i j ≡ s i j / m X 2 = 2 ( p i ⋅ p j ) / m X 2 y_{ij} \equiv s_{ij}/m_X^2 = 2(p_i \cdot p_j)/m_X^2 y ij ≡ s ij / m X 2 = 2 ( p i ⋅ p j ) / m X 2 [35 ] :
for X = V e c t o r X = \mathrm{Vector} X = Vector c N S = 0 c_{\mathrm{NS}} = 0 c NS = 0
for X = S c a l a r X = \mathrm{Scalar} X = Scalar c N S = + 2 c_{\mathrm{NS}} = +2 c NS = + 2
In Vincia, the baseline QED antenna function for fermion-antifermion dipoles [32 ] is chosen to have c N S = 0 c_{\mathrm{NS}}=0 c NS = 0
ant f f ˉ → f f ˉ γ d e f = ∣ M ( X 0 → f f ˉ γ ) ∣ 2 ∣ M ( X 0 → f f ˉ ) ∣ 2 ∣ c N S = 0 . \operatorname{ant}^{\mathrm{def}}_{f\bar{f}\to f\bar{f}\gamma}=\left.\frac{\left|M(X^{0}\to f\bar{f}\gamma)\right|^{2}}{\left|M(X^{0}\to f\bar{f})\right|^{2}}\right|_{c_{\mathrm{NS}}=0}. ant f f ˉ → f f ˉ γ def = M ( X 0 → f f ˉ ) 2 M ( X 0 → f f ˉ γ ) 2 c NS = 0 . Thus, for the vector J / ψ → μ + μ − J/\psi\to\mu^{+}\mu^{-} J / ψ → μ + μ −
R M E P S = Expanded Shower Matrix Element = ant f f ˉ → f f ˉ γ d e f ∣ M ( X 0 → f f ˉ ) ∣ 2 ∣ M ( X 0 → f f ˉ γ ) ∣ 2 (23) R_{\mathrm{ME}}^{\mathrm{PS}}=\frac{\text{Expanded Shower}}{\text{Matrix Element}}=\operatorname{ant}_{f\bar{f}\to f\bar{f}\gamma}^{\mathrm{def}}\frac{\left|M(X^{0}\to f\bar{f})\right|^{2}}{\left|M(X^{0}\to f\bar{f}\gamma)\right|^{2}} \tag*{(23)} R ME PS = Matrix Element Expanded Shower = ant f f ˉ → f f ˉ γ def M ( X 0 → f f ˉ γ ) 2 M ( X 0 → f f ˉ ) 2 (23) to be unity over all phase space. For the pseudoscalar B s 0 → μ + μ − B_s^{0}\to\mu^{+}\mu^{-} B s 0 → μ + μ − c N S = 0 c_{\mathrm{NS}}=0 c NS = 0 c N S = 2 c_{\mathrm{NS}}=2 c NS = 2
P M E C = ( R M E P S ) − 1 , P^{\mathrm{MEC}}=\left(R_{\mathrm{ME}}^{\mathrm{PS}}\right)^{-1}, P MEC = ( R ME PS ) − 1 , generalised so that it can be applied to multiple emissions, not just the first, following the same strategy as for (sector) QCD showers [33 ] . Unless otherwise stated, results labelled as “Vincia MEC” below include tree-level ME corrections for up to two photon emissions.
If we manually add a value of + 2 +2 + 2 R M E P S R_{\mathrm{ME}}^{\mathrm{PS}} R ME PS B s 0 → μ + μ − B_s^{0}\to\mu^{+}\mu^{-} B s 0 → μ + μ − J / ψ → μ + μ − J/\psi\to\mu^{+}\mu^{-} J / ψ → μ + μ − c N S c_{\mathrm{NS}} c NS 12
In Figure 13 , we compare the total emitted photon energy for J / ψ → μ + μ − J/\psi\to\mu^{+}\mu^{-} J / ψ → μ + μ − B s 0 → μ + μ − B_s^{0}\to\mu^{+}\mu^{-} B s 0 → μ + μ − + 2 +2 + 2 − 0.5 -0.5 − 0.5 O ( α 2 ) \mathcal{O}(\alpha^{2}) O ( α 2 ) J / ψ → μ + μ − J/\psi\to\mu^{+}\mu^{-} J / ψ → μ + μ − Figure 13a ), the effect of nonsingular terms and ME corrections in the soft region (low ∑ E γ \sum E_{\gamma} ∑ E γ ∼ m J / ψ / 4 ∼ 800 M e V \sim m_{J/\psi}/4\sim800\ \mathrm{MeV} ∼ m J / ψ /4 ∼ 800 MeV B s 0 → μ + μ − B_s^{0}\to\mu^{+}\mu^{-} B s 0 → μ + μ − Figure 13b ), we find strong agreement in the soft region up to m B s 0 / 4 ∼ 1.3 G e V m_{B_s^{0}}/4\sim1.3\ \mathrm{GeV} m B s 0 /4 ∼ 1.3 GeV + 2 +2 + 2 H → q γ q ˉ H\to q\gamma\bar{q} H → q γ q ˉ [35 ] . The important point is that, regardless of the baseline antenna function chosen, the framework for ME corrections developed in this work imposes the correct tree-level matrix elements in both decays.
Figure 13. Total photon energy radiated in (left) J / ψ → μ + μ − J/\psi\to\mu^{+}\mu^{-} J / ψ → μ + μ − B s 0 → μ + μ − B_s^0 \to\mu^{+}\mu^{-} B s 0 → μ + μ − + 2 +2 + 2 − 0.5 -0.5 − 0.5
We also compare the probability distribution of R M E P S R_{\mathrm{ME}}^{\mathrm{PS}} R ME PS Figure 14 shows this behaviour for the first two photon emissions in J / ψ → μ + μ − J/\psi\to\mu^{+}\mu^{-} J / ψ → μ + μ − B s 0 → μ + μ − B_s^{0}\to\mu^{+}\mu^{-} B s 0 → μ + μ −
Figure 14. Distribution of the base-10 logarithm of R M E P S R_{\mathrm{ME}}^{\mathrm{PS}} R ME PS (23) ), normalised to the total number of generated events and bin width, for the first and second photon emissions in (left) J / ψ → μ + μ − J/\psi\to\mu^{+}\mu^{-} J / ψ → μ + μ − B s 0 → μ + μ − B_s^0 \to\mu^{+}\mu^{-} B s 0 → μ + μ − c N S = 0 c_{\mathrm{NS}} = 0 c NS = 0 c N S = 2 c_{\mathrm{NS}} = 2 c NS = 2 ρ = log 10 ( 1 / 2 ) \rho= \log_{10}(1/2) ρ = log 10 ( 1/2 ) ρ = log 10 ( 1 / 1.5 ) \rho= \log_{10}(1/1.5) ρ = log 10 ( 1/1.5 ) ρ = log 10 ( 1 / 1.1 ) \rho= \log_{10}(1/1.1) ρ = log 10 ( 1/1.1 ) ρ = log 10 ( 1.1 ) \rho= \log_{10}(1.1) ρ = log 10 ( 1.1 ) ρ = log 10 ( 1.5 ) \rho= \log_{10}(1.5) ρ = log 10 ( 1.5 )
For J / ψ → μ + μ − J/\psi\to\mu^{+}\mu^{-} J / ψ → μ + μ − Figure 14a ), with the default c N S = 0 c_{\mathrm{NS}}=0 c NS = 0 x x x R M E P S = 1 R_{\mathrm{ME}}^{\mathrm{PS}}=1 R ME PS = 1 c N S = + 2 c_{\mathrm{NS}}=+2 c NS = + 2 c N S = + 2 c_{\mathrm{NS}}=+2 c NS = + 2 (23) , but broad tails remain where the genuine second-order ME correction is still required for a correct description.
For B s 0 → μ + μ − B_s^0 \to\mu^+\mu^- B s 0 → μ + μ − Figure 14b ), all events with c N S = + 2 c_{\mathrm{NS}}=+2 c NS = + 2 c N S = + 2 c_{\mathrm{NS}}=+2 c NS = + 2 (23) , but the shower still underestimates the matrix element, so the second-order correction is equally necessary to recover the correct expansion at this order. Again, the important point is that, regardless of the baseline antenna function chosen, the ME corrections enforce the correct tree-level matrix elements in both decays.
Results and validations
All results in this section are obtained using the Evtgen generator [17 , 18 , 19 ] to produce events prior to QED radiation. Evtgen is a simulation package specialised in the decays of heavy hadrons, incorporating a variety of models specific to different decay topologies. All of the four QED radiation generators we consider are employed as afterburners, acting on an identical set of input events containing undressed particles, thereby enabling a consistent comparison of the results obtained from the various QED radiation generators within a common framework.
We compare the QED radiation generators using a set of representative decay modes of neutral and charged mesons. Unless otherwise specified, all decays are simulated in Evtgen employing the phase-space model (PHSP). In § 7.1 , we consider two-body decay processes, focussing first on neutral meson decays:
J / ψ → e + e − J/\psi\to e^+e^- J / ψ → e + e − J / ψ → μ + μ − J/\psi\to\mu^+\mu^- J / ψ → μ + μ − → \to →
B s 0 → μ + μ − B_s^0 \to\mu^+\mu^- B s 0 → μ + μ −
B 0 → π + π − B^0 \to\pi^+\pi^- B 0 → π + π − B 0 → K + π − B^0 \to K^+\pi^- B 0 → K + π − B 0 → D s + D − B^0 \to D_s^+D^- B 0 → D s + D −
D 0 → π + π − D^0 \to\pi^+\pi^- D 0 → π + π − D 0 → K − π + D^0 \to K^-\pi^+ D 0 → K − π + D 0 → K + K − D^0 \to K^+K^- D 0 → K + K −
and then focussing on charged meson decays:
In § 7.3 , we consider three-body decay processes:
B 0 → D 0 K + π − B^0 \to D^0K^+\pi^- B 0 → D 0 K + π −
B + → π + π 0 π 0 B^+ \to\pi^+\pi^0\pi^0 B + → π + π 0 π 0 B + → D + K − π + B^+ \to D^+K^-\pi^+ B + → D + K − π +
D + → π + π 0 π 0 D^+ \to\pi^+\pi^0\pi^0 D + → π + π 0 π 0
D + → π + π + π − D^+ \to\pi^+\pi^+\pi^- D + → π + π + π − D + → K + π + π − D^+ \to K^+\pi^+\pi^- D + → K + π + π − [73 , 74 ] .
Finally, in § 7.4 , we present dedicated illustrations of multipole radiation patterns for a specific Born-level phase-space point for D + → π + π + π − D^+ \to\pi^+\pi^+\pi^- D + → π + π + π − B + → D + K + π − B^+ \to D^+K^+\pi^- B + → D + K + π −
To avoid a proliferation of plots, we summarise many results in tables showing the total emitted photon energy (∑ E γ \sum E_\gamma ∑ E γ
Two-body decays
Two-body decays of neutral hadrons
Leptonic decays of neutral vector particles such as the J / ψ J/\psi J / ψ Z Z Z J / ψ → e + e − J/\psi\to e^{+}e^{-} J / ψ → e + e −
Figure 15a shows the total emitted photon energy. For this decay, Vincia without MECs matches Pythia closely. This is expected because Pythia includes an explicit first-order ME correction, while, as discussed in § 6 , Vincia’s baseline spin-1/2 antenna is derived from the Z → ℓ + ℓ − Z \to\ell^{+}\ell^{-} Z → ℓ + ℓ − Figure 14a in § 6 ), which shows that the baseline antenna underestimates the second-emission matrix element, giving a net positive ME correction. In Figure 15b we verify that we can reproduce the plot from ref. [4], where the first-order ME correction mode (Photons 21) sits between the pure scalar QED (Photons 1) and the scalar QED + collinear correction (Photons 20). Comparing the three generators (Photons 21, Vincia, and Photos) in Figure 15c , we find good agreement in the single-emission region (∑ E γ ≤ m J / ψ / 2 ∼ 1.55 \sum E_{\gamma} \le m_{J/\psi}/2 \sim1.55 ∑ E γ ≤ m J / ψ /2 ∼ 1.55 m J / ψ / 2 m_{J/\psi}/2 m J / ψ /2
In Figure 15d , we show the angular distribution of the total emitted photon energy. Here and in all subsequent angular-distribution plots in this section, the opening angle is defined as the angle between the radiated photon and one of the post-branching final-state daughters (see figure captions and discussion in § 4.3 ). Enhancements in the collinear regions (0 ∘ 0^\circ 0 ∘ 180 ∘ 180^\circ 18 0 ∘ 180 ∘ 180^\circ 18 0 ∘ Figure 15c ). The rest are much more similar, with Photos radiating slightly more and Photons 21 slightly less than Vincia at wide angles.
Table 4 shows good agreement between all generators, both on average energy loss and high-energy-loss fractions, with differences greater than 5% only between different Photons configurations. Additionally, the second-order correction in Vincia MEC increases the total energy loss slightly, by about 1.5% relative to Pythia and the baseline Vincia, or by about 2.6% relative to Photons 21.
Table 4. Average amounts of photon energy radiated in J / ψ → e + e − J/\psi\to e^{+}e^{-} J / ψ → e + e − § 7 .
Table 5 shows an equivalent summary table for J / ψ → μ + μ − J/\psi\to\mu^{+}\mu^{-} J / ψ → μ + μ − Figure 14a . However, we note that the total energy loss predicted by Photos is about 5% greater than that generated by Photons 21.
Table 5. Average amounts of photon energy radiated in J / ψ → μ + μ − J/\psi\to\mu^{+}\mu^{-} J / ψ → μ + μ − § 7 .
Moving to B 0 → μ + μ − B^{0} \to\mu^{+}\mu^{-} B 0 → μ + μ − Figure 16a clearly shows the impact of nonsingular terms in the ME corrections in Vincia, as expected from the studies in § 6 . In Figure 16b , Photons 1 (pure YFS) radiates much less in the tail than Photons 20 and 21, which overlap with each other. Comparing the three generators in Figure 16c , we find that, unlike the V 0 → ℓ + ℓ − V^{0} \to\ell^{+}\ell^{-} V 0 → ℓ + ℓ −
(a) Total photon energy radiated in J / ψ → e + e − J/\psi\to e^{+}e^{-} J / ψ → e + e −
(b) Total photon energy radiated in J / ψ → e + e − J/\psi\to e^{+}e^{-} J / ψ → e + e −
(c) Total photon energy radiated in J / ψ → e + e − J/\psi\to e^{+}e^{-} J / ψ → e + e −
(d) Angular distribution of photon energy in J / ψ → e + e − J/\psi\to e^{+}e^{-} J / ψ → e + e −
Figure 15. Total photon energy radiated and angular distribution of photon energy in J / ψ → e + e − J/\psi\to e^{+}e^{-} J / ψ → e + e −
radiates more than the others in the hard tail. Similarly to the V 0 → ℓ + ℓ − V^0 \to\ell^+\ell^- V 0 → ℓ + ℓ − Figure 16d ), we now observe the dead-cone effect around the muon directions at 0 ∘ 0^\circ 0 ∘ 180 ∘ 180^\circ 18 0 ∘ B s 0 → μ + μ − B_s^0 \to\mu^+\mu^- B s 0 → μ + μ − Table 6 , where all other generators predict somewhat smaller total radiation losses. Nevertheless, Vincia MEC and Photons 21 agree well on the fractions of events that have total energy losses greater than 1, 10, and 100 MeV.
Table 6. Average amounts of photon energy radiated in B s 0 → μ + μ − B_s^0 \to\mu^+\mu^- B s 0 → μ + μ −
(a) Total photon energy radiated in B s 0 → μ + μ − B_s^0 \to\mu^+\mu^- B s 0 → μ + μ −
(b) Total photon energy radiated in B s 0 → μ + μ − B_s^0 \to\mu^+\mu^- B s 0 → μ + μ −
(c) Total photon energy radiated in B s 0 → μ + μ − B_s^0 \to\mu^+\mu^- B s 0 → μ + μ −
(d) Angular distribution of energy in B s 0 → μ + μ − B_s^0 \to\mu^+\mu^- B s 0 → μ + μ −
Figure 16. Total photon energy radiated and angular distribution of photon energy in B s 0 → μ + μ − B_s^0 \to\mu^+\mu^- B s 0 → μ + μ −
(a) Total photon energy radiated in B 0 → K + π − B^0 \to K^+\pi^- B 0 → K + π −
(b) Total photon energy radiated in B 0 → K + π − B^0 \to K^+\pi^- B 0 → K + π −
(c) Total photon energy radiated in B 0 → K + π − B^0 \to K^+\pi^- B 0 → K + π −
(d) Angular distribution of photon energy in B 0 → K + π − B^0 \to K^+\pi^- B 0 → K + π −
Figure 17. Total photon energy radiated and angular distribution of photon energy in B 0 → K + π − B^0 \to K^+\pi^- B 0 → K + π −
For purely (pseudo)scalar B 0 → K + π − B^0 \to K^+\pi^- B 0 → K + π − Figure 17a shows that the main Pythia shower still treats the kaon and pion as fermions (since it lacks scalar splitting functions). This is why, in released Pythia versions, the QED shower is only available for purely leptonic hadron decays. As Figure 15a shows, Vincia MEC is around 2% higher than the uncorrected configuration, while having around 0.1% fewer entries in the first bin (0 – 50 MeV). For the three Photons configurations in Figure 17b , we used the same random-number seed in each run to show that, for pure (pseudo)scalar decays, switching between these configurations produces no difference. Comparing the generators in Figure 17c and Table 7 shows excellent agreement, with only Vincia MEC giving slightly higher radiation losses than the others (except Pythia, which, as noted above, is not a realistic model for these decays). This agreement is reflected in the angular distributions (Figure 17d ), where all generators agree very well. The visible asymmetry in the dead-cone effect is caused by the difference between the kaon and pion masses. This pattern repeats overall for the other two-body decays of neutral hadrons considered in this study, summarised in Table 8 –Table 12 .
Table 7. Average amounts of photon energy radiated in B 0 → K + π − B^0 \to K^+\pi^- B 0 → K + π −
Table 8. Average amounts of photon energy radiated in B 0 → D s + D − B^0 \to D_s^+D^- B 0 → D s + D − § 7 .
Table 12. Average amounts of photon energy radiated in D 0 → π + π − D^{0} \to\pi^{+}\pi^{-} D 0 → π + π −
Table 9. Average amounts of photon energy radiated in B 0 → π + π − B^0 \to\pi^+\pi^- B 0 → π + π − § 7 .
Table 10. Average amounts of photon energy radiated in D 0 → K + K − D^0 \to K^+K^- D 0 → K + K − § 7 .
Table 11. Average amounts of photon energy radiated in D 0 → K − π + D^{0} \to K^{-}\pi^{+} D 0 → K − π + § 7 .
Two-body decays of charged hadrons
In Vincia, the description of radiative corrections in decays of charged initial states is based on its unique “resonance-final” antennae, which were initially developed for top-quark decays in [36 ] and generalised to QED in [32 ] . In most other generators we are aware of, radiation in decays of charged initial states is treated by considering the initial hadron to be at rest and hence non-radiating, and then including explicit radiation functions only for the final-state decay products. In singular limits, these treatments should be equivalent.
Taking B + → π + π 0 B^{+} \to\pi^{+}\pi^{0} B + → π + π 0 Figure 18a shows that the two Vincia configurations agree well, while Pythia’s fermion-radiator treatment yields a significantly harder radiation spectrum. As in the case of neutral decays, the ME corrected antennae nearly overlap with the uncorrected antennae. In Figure 18b , we again used the same random-number seeds to show that the Photons distributions remain unchanged in the three configurations. Comparing the three generators in Figure 18c , we find good agreement up to the multi-emission tail near 2 GeV, where Photos MEC and Vincia radiate more than Photons. This is consistent with Table 13 : only Vincia MEC shows a slightly higher total energy loss, about 400 keV more than Photons 21.
Table 13. Average amounts of photon energy radiated in B + → π + π 0 B^{+} \to\pi^{+}\pi^{0} B + → π + π 0 § 7 .
In the angular distribution (Figure 18d ), the radiation pattern is now highly asymmetric since only one of the decay products is charged. The agreement between the generators is very strong in the small-angle region, with Photons 21 emitting less in the wide-angle region than the other generators.
In Table 14 and Table 15 , we summarise the radiation losses for B ± → K ± π 0 B^{\pm} \to K^{\pm}\pi^{0} B ± → K ± π 0 D ± → π ± π 0 D^{\pm} \to\pi^{\pm}\pi^{0} D ± → π ± π 0
Table 14. Average amounts of photon energy radiated in B + → K + π 0 B^{+} \to K^{+}\pi^{0} B + → K + π 0 § 7 .
Table 15. Average amounts of photon energy radiated in D + → π + π 0 D^{+} \to\pi^{+}\pi^{0} D + → π + π 0 § 7 .
Three-body decays
Three-body decays, and many-body decays generally, differ greatly from two-body decays because the configuration of radiating charges is not fixed as a single dipole of fixed invariant mass. Instead, they allow multiple orientations and relative energies of the decay products and possibly multiple charged antennae. We consider three examples of three-body decay topologies:
Neutral particle decays to two charged and one neutral particle, e.g. B 0 → D 0 K + π − B^0 \to D^0 K^+ \pi^- B 0 → D 0 K + π −
Charged particle decays to one charged and two neutral particles, e.g. B + → π + π 0 π 0 B^+ \to\pi^+ \pi^0 \pi^0 B + → π + π 0 π 0 D + → π + π 0 π 0 D^+ \to\pi^+ \pi^0 \pi^0 D + → π + π 0 π 0
Charged particle decays to three charged particles, e.g. B + → D + K − π + B^+ \to D^+ K^- \pi^+ B + → D + K − π + D + → K + π + π − D^+ \to K^+ \pi^+ \pi^- D + → K + π + π − D + → π + π − π + D^+ \to\pi^+ \pi^- \pi^+ D + → π + π − π +
Only the last case should exhibit genuine multipole effects beyond dipole terms.
In this study, QED radiation is handled by treating Born-level decays as point-like and direct, without any special treatment of sequential decays via intermediate resonances (e.g., ρ \rho ρ K ∗ K^* K ∗ D ∗ D^* D ∗
B 0 → D 0 K + π − B^0 \to D^0 K^+ \pi^- B 0 → D 0 K + π − Figure 19 shows the comparisons for B 0 → D 0 K + π − B^0 \to D^0 K^+ \pi^- B 0 → D 0 K + π − Figure 19 a shows that Vincia MEC lies slightly above the uncorrected shower, as in the two-body decays. The dipole approximation matches the uncorrected shower within statistical fluctuations, as expected for a decay with only one antenna. In fig. 19b, all Photons configurations yield again identical plots as in the fully scalar cases. Comparing all generators in fig. 19c shows, unlike in the B 0 → K + π − B^0 \to K^+ \pi^- B 0 → K + π −
Figure 19. (a) Total photon energy radiated in B 0 → D 0 K + π − B^0 \to D^0K^+\pi^- B 0 → D 0 K + π − B 0 → D 0 K + π − B^0 \to D^0K^+\pi^- B 0 → D 0 K + π − B 0 → D 0 K + π − B^0 \to D^0K^+\pi^- B 0 → D 0 K + π − B 0 → D 0 K + π − B^0 \to D^0K^+\pi^- B 0 → D 0 K + π −
Comparing the summary table for B 0 → D 0 K + π − B^0 \to D^0 K^+ \pi^- B 0 → D 0 K + π − B 0 → K + π − B^0 \to K^+ \pi^- B 0 → K + π −
Table 16. Average amounts of photon energy radiated in B 0 → D 0 K + π − B^0 \to D^0 K^+ \pi^- B 0 → D 0 K + π −
B + → π + π 0 π 0 B^{+} \to\pi^{+}\pi^{0}\pi^{0} B + → π + π 0 π 0 D 0 → π + π 0 π 0 D^{0} \to\pi^{+}\pi^{0}\pi^{0} D 0 → π + π 0 π 0 We next consider X + → π + π 0 π 0 X^{+} \to\pi^{+}\pi^{0}\pi^{0} X + → π + π 0 π 0 X = B X = B X = B D D D B + → π + π 0 π 0 B^{+} \to\pi^{+}\pi^{0}\pi^{0} B + → π + π 0 π 0 Figure 20 with that for the two-body decay with one fewer neutral pion (Figure 18 ), we observe qualitatively similar generator behaviour, but slightly different tails: for large energy losses, Vincia sits slightly above Photons, while Photos (as for B 0 → D 0 K + π − B^{0} \to D^{0}K^{+}\pi^{-} B 0 → D 0 K + π − Figure 19 c) decays) sits below it. However, the summary in Table 17 shows that these differences are not particularly significant.
(a) Total photon energy radiated in B + → π + π 0 B^+ \to\pi^+\pi^0 B + → π + π 0
(b) Total photon energy radiated in B + → π + π 0 B^+ \to\pi^+\pi^0 B + → π + π 0
(c) Total photon energy radiated in B + → π + π 0 B^+ \to\pi^+\pi^0 B + → π + π 0
(d) Angular distribution of photon energy in B + → π + π 0 B^+ \to\pi^+\pi^0 B + → π + π 0
Figure 18. Total photon energy radiated and angular distribution of photon energy in B + → π + π 0 B^+ \to\pi^+\pi^0 B + → π + π 0
Table 17. Average amounts of photon energy radiated in B + → π + π 0 π 0 B^{+} \to\pi^{+}\pi^{0}\pi^{0} B + → π + π 0 π 0 § 7 .
(a) Total photon energy radiated in B + → π + π 0 π 0 B^+ \to\pi^+\pi^0\pi^0 B + → π + π 0 π 0
(b) Total photon energy radiated in B + → π + π 0 π 0 B^+ \to\pi^+\pi^0\pi^0 B + → π + π 0 π 0
(c) Total photon energy radiated in B + → π + π 0 π 0 B^+ \to\pi^+\pi^0\pi^0 B + → π + π 0 π 0
Figure 20. Total photon energy radiated in B + → π + π 0 π 0 B^+ \to\pi^+\pi^0\pi^0 B + → π + π 0 π 0
Analogous results for D 0 → π + π 0 π 0 D^{0} \to\pi^{+}\pi^{0}\pi^{0} D 0 → π + π 0 π 0 Table 18 , now with the somewhat lower total energies in Photos showing up in orange shading. However, since the total radiated energy is quite small for this decay, the absolute difference is less than 100 keV on average.
Table 18. Average amounts of photon energy radiated in D + → π + π 0 π 0 D^{+} \to\pi^{+}\pi^{0}\pi^{0} D + → π + π 0 π 0 § 7 .
B + → D + K − π + B^{+} \to D^{+}K^{-}\pi^{+} B + → D + K − π + D + → K + π + π − D^{+} \to K^{+}\pi^{+}\pi^{-} D + → K + π + π − D + → π + π + π − D^{+} \to\pi^{+}\pi^{+}\pi^{-} D + → π + π + π − For fully charged three-body topologies such as B + → D + K − π + B^{+} \to D^{+}K^{-}\pi^{+} B + → D + K − π + § 2 and § 3 , these effects are taken into account in the YFS formalism in the soft limit and in Vincia’s (default) multipole shower mode, while Vincia’s dipole mode approximates them in terms of a direct sum of two (non-interfering) dipoles, identified based on a principle of maximal screening [32 ] .
Figure 21a shows that Vincia with and without MECs exhibit the same behaviour found before, with the ME corrected shower slightly above the uncorrected one. Vincia’s dipole approximation sits consistently above the other Vincia modes, especially in the low ∑ E γ \sum E_{\gamma} ∑ E γ
This behaviour is consistent across the three considered fully-charged three-body decays, summarised in Table 19 , Table 20 , and Table 21 . These tables also show that the relation between Photos and Photons 21 is relatively sporadic: Photos under-radiates (Table 19 ), over-radiates (Table 20 ), or radiates within our 2% bound (Table 21 ). The amount of this radiation seems to correlate negatively with the relative size of the phase space.
Table 19. Average amounts of photon energy radiated in B + → D − K + π + B^{+} \to D^{-}K^{+}\pi^{+} B + → D − K + π + § 7 .
Table 20. Average amounts of photon energy radiated in D + → K + π + π − D^{+} \to K^{+}\pi^{+}\pi^{-} D + → K + π + π − § 7 .
Table 21. Average amounts of photon energy radiated in D + → π + π + π − D^{+} \to\pi^{+}\pi^{+}\pi^{-} D + → π + π + π − § 7 .
(a) Total photon energy radiated in B + → D + K − π + B^{+} \to D^{+}K^{-}\pi^{+} B + → D + K − π +
(b) Total photon energy radiated in B + → D + K − π + B^{+} \to D^{+}K^{-}\pi^{+} B + → D + K − π +
(c) Total photon energy radiated in B + → D + K − π + B^{+} \to D^{+}K^{-}\pi^{+} B + → D + K − π +
Figure 21. Total photon energy radiated in B + → D + K − π + B^{+} \to D^{+}K^{-}\pi^{+} B + → D + K − π +
A fixed three-body multipole configuration
We inspect the angular distributions of photons in decays to more than two bodies by considering a useful fixed phase-space point referred to as “Mercedes” configuration, illustrated in Figure 22 . This configuration serves as direct Born-level input to Pythia, and is then dressed with QED radiation showered with Vincia. This configuration should maximise multipole effects, highlighting differences between Vincia’s multipole and dipole modes.
Figure 22. The Mercedes configuration used for the fixed decay of a charged mother into three charged daughters. The positive-sign (constructive) antennae are labelled in light blue and the negative-sign (destructive) ones in light red. R stands for resonance (initial state), while F stands for final state.
We first consider the case where all three daughter particles have the same mass, using D + → π + π + π − D^{+} \to\pi^{+}\pi^{+}\pi^{-} D + → π + π + π −
To illustrate the sectorisation in Vincia’s multipole mode, fig. 23b shows the relative contributions of each of the three final-final QED sectors to the full emission pattern for Vincia MEC 13 π + π + \pi^{+}\pi^{+} π + π + π + π − \pi^{+}\pi^{-} π + π − π + π − \pi^{+}\pi^{-} π + π −
Another interesting feature of fig. 23a is that it exhibits a change in sign of the multipole corrections. In the angular region in-between the two positive pions (from 300° to 60°), the multipole distribution is below the dipole approximation, reflecting predominantly negative (destructive) multipole corrections. While in the complementary angular region between 60° and 300°, the multipole distribution is generally higher than the dipole approximation, reflecting predominantly positive (constructive) multipole corrections. The lack of the full multipole structure in the dipole case is also visible in the region from 60° to 180° (and, identically, 180 to 300 ∘ 300^\circ 30 0 ∘ 60 ∘ 60^\circ 6 0 ∘ 180 ∘ 180^\circ 18 0 ∘
Turning now to the case of different daughter masses, we consider B + → D + K + π − B^+ \to D^+K^+\pi^- B + → D + K + π −
p B + = ( 5.279 , 0.000 , 0.000 , 0.000 ) T G e V , p_{B^+} = (5.279,\ 0.000,\ 0.000,\ 0.000)^{\mathrm{T}}\ \mathrm{GeV}, p B + = ( 5.279 , 0.000 , 0.000 , 0.000 ) T GeV , p D + = ( 2.348 , 0.710 , − 1.230 , 0.000 ) T G e V , p_{D^+} = (2.348,\ 0.710,\ -1.230,\ 0.000)^{\mathrm{T}}\ \mathrm{GeV}, p D + = ( 2.348 , 0.710 , − 1.230 , 0.000 ) T GeV , p K + = ( 1.504 , 0.710 , 1.230 , 0.000 ) T G e V , p_{K^+} = (1.504,\ 0.710,\ 1.230,\ 0.000)^{\mathrm{T}}\ \mathrm{GeV}, p K + = ( 1.504 , 0.710 , 1.230 , 0.000 ) T GeV , p π − = ( 1.427 , − 1.421 , 0.000 , 0.000 ) T G e V . p_{\pi^-} = (1.427,\ -1.421,\ 0.000,\ 0.000)^{\mathrm{T}}\ \mathrm{GeV}. p π − = ( 1.427 , − 1.421 , 0.000 , 0.000 ) T GeV . Figure 24 shows features similar to the equal-mass case: the asymmetry in the peaks around the mass-suppressed dead-cone regions of the pion (180 ∘ 180^\circ 18 0 ∘ 60 ∘ 60^\circ 6 0 ∘ 300 ∘ 300^\circ 30 0 ∘ D D D Figure 24 b, the region around 300 ∘ 300^\circ 30 0 ∘ 300 ∘ 300^\circ 30 0 ∘ K + D + K^+D^+ K + D + K + π − K^+\pi^- K + π − D + K + D^+K^+ D + K + D + π + D^+\pi^+ D + π + D D D D + K + D^+K^+ D + K + D + π + D^+\pi^+ D + π +
Figure 24. (a) Angular distribution of photon energy for the first emission in Vincia, Vincia MEC and Vincia’s dipole approximation. (b) Contributions from individual antennae in each sector for the first emission in Vincia MEC. Angular distribution of radiated photon energy for the first emission in B + → D + K + π − B^{+} \to D^{+}K^{+}\pi^{-} B + → D + K + π − Figure 22 ). The π − \pi^{-} π − 180 ∘ 180^{\circ} 18 0 ∘ K + K^{+} K + 60 ∘ 60^{\circ} 6 0 ∘ D + D^{+} D + 300 ∘ 300^{\circ} 30 0 ∘ D + K + D^{+}K^{+} D + K + K + π − K^{+}\pi^{-} K + π − D + π − D^{+}\pi^{-} D + π −
The fact that the ME-corrected and -uncorrected distributions agree implies that we do not regard this as in principle problematic; in multibody decays, different sectorisations produce different “histories”; different definitions of which radiative phase-space points are mapped to which Born-level phase-space points. Here, we artificially consider only a single specific Born-level phase-space point. Our sectorisation then defines which hypersurface in the Born+γ \gamma γ
Vector (PDF)
caption: (a) Angular distribution of photon energy for the first emission in Vincia, Vincia MEC and Vincia’s dipole approximation. (b) Contributions from each FF sector to the angular distribution of energy for the first emission in Vincia MEC. fig:23 : Angular distribution of radiated photon energy for the first emission in D + → π + π + π − D^{+} \to\pi^{+}\pi^{+}\pi^{-} D + → π + π + π − Figure 22 ) with the π − \pi^{-} π − 180 ∘ 180^\circ 18 0 ∘ π + \pi^{+} π + 60 ∘ 60^\circ 6 0 ∘ 300 ∘ 300^\circ 30 0 ∘ π + π + \pi^{+}\pi^{+} π + π + π + π − \pi^{+}\pi^{-} π + π − π + \pi^{+} π + 60 ∘ 60^\circ 6 0 ∘ π + π − \pi^{+}\pi^{-} π + π − π + \pi^{+} π + 300 ∘ 300^\circ 30 0 ∘
:::
phase space that point is mapped onto. Sector boundaries correspond to discontinuities in that map; but since the full matrix element is smooth, the ME-corrected distribution considered as a function defined on the full Born+γ \gamma γ
Another notable feature across both figs. 23b and 24b is that the only contributing sectors are final-final ones. This is due to all the resonance-final sectors being vetoed in the sectorisation process in the Mercedes configuration.
In summary, in the dipole approximation, non-singular contributions from other dipole antennae in the event can partly “fill in” the dead cones of the most suppressed massive radiators, while the multipole shower (with or without MECs) captures the full set of soft interference effects as well as (all of) the quasi-collinear limits.
Performance
To evaluate the computational efficiency of the QED radiation generators we have considered in this work, we benchmark both CPU execution time and peak memory consumption through the Evtgen generator framework. The generation of events is conducted using Evtgen release 3.0.0, which provides interfaces to Photos release 3.64, Sherpa release 3.0.1 (with physics consistency verified against release 2.2.15), and Pythia 8.315. The software was compiled within a standardised LHC Computing Grid environment (LCG_107, platform x86_64-el9-clang16-opt) distributed via the Cern Virtual Machine File System (CVMFS) [75 ] .
We evaluate benchmarking metrics for two operational regimes: single-threaded execution to establish baseline single-core performance, and multi-threaded processing to quantify scaling and speed-up behaviour enabled by recent multi-threading adaptations in Evtgen [18 , 19 ] . We perform benchmark tests for four representative decay topologies: J / ψ → e + e − J/\psi\to e^{+}e^{-} J / ψ → e + e − J / ψ → μ + μ − J/\psi\to\mu^{+}\mu^{-} J / ψ → μ + μ − B 0 → K + π − B^{0} \to K^{+}\pi^{-} B 0 → K + π − D + → K + π + π − D^{+} \to K^{+}\pi^{+}\pi^{-} D + → K + π + π − 10 7 10^{7} 1 0 7 10 7 10^{7} 1 0 7
Table 22 reports the elapsed time per event and the corresponding ratio of the runtime with FSR simulation relative to the baseline configuration without FSR, for single-threaded execution. The runtime ratios are also presented in Figure 25 . Table 23 reports the results obtained for a single-threaded execution. The difference in peak RSS between runs performed with and without FSR simulation is presented in fig. 23.
Table 22. Absolute timing (ms per event) and overhead factor t / t n o F S R t/t_{\mathrm{noFSR}} t / t noFSR
Table 23. Maximum resident set size memory (max RSS, in MB) for single-thread execution (10 7 10^{7} 1 0 7
Figure 25. Overhead factor t / t n o F S R t/t_{\mathrm{noFSR}} t / t noFSR
Overall, we find that the FSR simulation using Vincia achieves run times comparable to the fastest available plug-ins. When ME corrections are enabled, the execution time increases by approximately a factor of two. Nevertheless, the resulting generator performance remains competitive. With respect to memory usage, enabling Vincia increases the memory footprint of the Evtgen simulation by only about 10%.
Multi-threaded tests use the C++ std threading library, with one Evtgen instance created per thread. In each test, the 10 7 10^{7} 1 0 7 J / ψ → e + e − J/\psi\to e^{+}e^{-} J / ψ → e + e − Figure 27 shows the speed-up, defined as the ratio of single-threaded to multi-threaded wall-clock time, while Figure 28 shows the corresponding RSS memory consumption normalised to the single-threaded value. The figures corresponding to the other modes are collected in app. § 14 . Although some performance differences are observed among the various decay modes, the overall behaviour is consistent with that found for the J / ψ → e + e − J/\psi\to e^{+}e^{-} J / ψ → e + e −
Figure 27. Speed-up curves for J / ψ → e + e − J/\psi\to e^{+}e^{-} J / ψ → e + e −
(a) Vincia
(b) Photons
(c) FSR without MEC vs. no FSR
(d) FSR with MEC (where available) vs. no FSR
Figure 28. Maximum RSS memory consumption relative to single-thread execution for J / ψ → e + e − J/\psi\to e^{+}e^{-} J / ψ → e + e −
The Photos and Photons libraries are not thread-safe. Therefore, their invocations are serialised via a thread lock inside Evtgen. For Photos, both the event-record conversion and the shower call must be locked, whereas for Photons only the shower call requires serialisation. As a result, multi-threading yields negligible speed-up for these generators, as threads stall waiting for the lock.
When Vincia or Pythia are used, Evtgen creates one Pythia instance per thread. Since these instances are thread-local, no serialisation is required, and a near-ideal linear speed-up is achieved. The cost is a ∼ \sim ∼
In general, Vincia performs competitively, even with ME corrections, while significantly reducing the CPU time of FSR simulation in Evtgen compared to the current default configuration. In multithreaded-threaded runs, Vincia enables near-ideal Evtgen speed-up with a small memory increase, provided it operates in the normally multithreaded regime (no hyper-threading). With hyper-threading, the speed-up is still significant but less linear and more saturated.
Figure 26. Difference in RSS memory between simulation with and without FSR simulated for single-threaded execution.
Summary and outlook
Precise descriptions of hadron decays require a reliable treatment of the QED radiation that accompanies them. We have presented a new framework for QED antenna showers in hadron decays, implemented in the Vincia shower module [34 ] of Pythia 8 [16 ] .
Our antenna functions include full mass dependence and spin-dependent (quasi-)collinear terms for spin-0, spin-1/2, and spin-1 emitters (spin-3/2 ones are so far treated the same as spin-1/2 ones). For soft photon emissions, our antenna functions reproduce the full coherent multipole eikonal structure [31 , 32 ] with the same fidelity as the YFS formalism [20 ] . Unlike YFS, however, they also reproduce the correct DGLAP kernels in collinear limits. By default, the full soft multipole pattern is included via phase-space sectorisation [32 , 33 ] ; a faster dipole option based on maximal screening is also available.
To improve the accuracy of the shower, in particular for semi-hard emissions (which can be important since they can be associated with significant photon energy losses), we have developed a semi-automated tool chain for iterated tree-level ME corrections (MECs). Starting from an effective Lagrangian for generic hadronic degrees of freedom, implemented in FeynRules, we use MadGraph to produce a library of plugins that Vincia calls at run time. By default, this library currently includes MECs for up to two photon emissions, for decays involving only (pseudo)scalars and for decays of (pseudo)scalars and vectors to lepton pairs, see app. § 12 . Since the ME correction factor is the ratio of the full matrix element to the shower kernel, the corrected result does not depend on the choice of baseline antenna function. We illustrated this by varying the nonsingular term of the antenna function for J / ψ → μ + μ − J/\psi\to\mu^{+}\mu^{-} J / ψ → μ + μ − B s 0 → μ + μ − B_s^0 \to\mu^{+}\mu^{-} B s 0 → μ + μ −
We also advocated cos θ d e a d = β \cos\theta_{\mathrm{dead}} = \beta cos θ dead = β 1 / γ 1/\gamma 1/ γ D D D B → D K B \to DK B → D K
In § 5 , we emphasised that the number of emitted photons is not an infrared safe quantity; it depends entirely on the infrared cutoff and on the variable it is imposed on, and both differ between generators. Only infrared-safe quantities, defined with respect to a common notion of a resolved photon, can be meaningfully compared, with cutoffs chosen sufficiently low to resolve the full phase space of such photons. In this context, we studied the effect of varying the cutoff values in various generators and observed an unexplained suppression of the total emitted photon energy in Photos [10 , 11 , 12 , 13 , 14 ] for J / ψ → e + e − J/\psi\to e^{+}e^{-} J / ψ → e + e −
Using EvtGen [17 , 18 , 19 ] to provide identical Born-level events, we compared Vincia to Photos and Sherpa’s Photons++ module [4 ] for a range of two- and three-body decays. Pythia’s internal QED shower was also included for some of these comparisons, although, by construction, it treats all radiators as spin-1/2 particles. For the total radiated energy, and for the fractions of events losing more than 1, 10, and 100 MeV, Vincia typically agrees with Photos with first-order MECs to within a few percent. The second-order MECs appear to consistently increase the average energy loss, by about 1–3%.
For decays to three charged particles, the dipole approximation overestimates the average energy loss of the multipole shower by 5–20%. This indicates that the multipole corrections are predominantly destructive (at least relative to the maximal-screening dipole approximation). In a fixed Mercedes configuration (see § 7.4 ), the dipole approximation also partly fills in the dead cones, since the “other” dipoles are allowed to radiate into them.
Regarding computational performance, Vincia achieves single-core run times comparable to the currently available final-state radiation (FSR) plug-ins (Sherpa’s Photons module and Photos), with ME corrections increasing the execution time by approximately a factor of two while remaining competitive overall. The memory footprint of EvtGen increases only by ∼10% when Vincia is enabled. In multi-threaded operation, the Photos and Photons libraries require serialisation of their calls via a thread lock, yielding negligible speed-up with increasing thread count. In contrast, Vincia and Pythia can be used as thread-local instances, achieving a near-ideal linear speed-up at the cost of a ∼ 10 \sim10 ∼ 10 30 % 30\% 30%
Several extensions are open for future work. These include radiation from spin-3 / 2 3/2 3/2 [55 ] to intermediate hadronic resonances (such as intermediate ρ \rho ρ K ∗ K^* K ∗
The present framework, for pseudoscalars, will be included in the next Pythia release (tentatively 8.319).
Acknowledgements
We are grateful to M. Kreps, T. Latham, and J. Back for their valuable assistance in reviewing the Evtgen-Vincia interface code and for their insightful discussions. We also gratefully acknowledge fruitful discussions related to this work with U. Egede, G. Valencia, M. Schoenherr, and P. Richardson. Regarding the technical aspects of multithreading, we also acknowledge valuable discussions and technical support provided by the Research Software Engineers at the University of Warwick, Heather Ratcliffe and Christopher Brady. FA is supported by UKRI FLF grant MR/W009048/1. PS is supported by ARC grant DP230103014. GM is supported by the Commonwealth through an Australian Government Research Training Program Scholarship [DOI: https://doi.org/10.82133/C42F-K220]. This work was also supported by the Monash-Warwick Alliance in Particle Physics (MWAPP).
In formulating the summary and outlook section (and only that section), we made partial use of suggestions from Claude AI. All statements were explicitly edited by the authors, and we take full responsibility for the statements made in that section.
Implementation of the QED shower algorithm
A description of the antenna formalism underlying Vincia’s QED shower is given in sec. 3. The shower is embedded in Pythia and can be accessed via the generic shower models of Pythia (base classes ShowerModel, TimeShower, and SpaceShower) in three main ways:
As part of the standard interleaved QCD+QED perturbative evolution of a partonic hard-scattering event, via Vincia’s main prepare(), pTnext(), and branch() methods;
As a continuation of the evolution of hadronic collisions, after adding beam remnants but before hadronisation and particle decays, via showerQEDafterRemnants();
As a plugin for QED showering in hadron decays via the showerQED() method, callable within Pythia at hadron decay generation or directly from external programmes such as Evtgen.
The primary QED radiation algorithm for hadron decays is Pythia’s TimeShower-inherited showerQED() method, with the signature:
showerQED(int iBeg, int iEnd, Event& event, double pTmax, int nBranchMax);Figure 29 shows an overview of the functions in the showerQED algorithm. As a first step, the procedure calls a preparation routine called VinciaQED::prepare, which discriminates between hadronic decays and configurations involving leptons and quarks. This distinction is required because, in the case of bare partonic constituents, additional interactions (such as QCD) may need to be interleaved. Prepare also determines whether the initial system provided has MECs (MECs::prepare and MECs::meAvailable) and builds the internal systems (buildSystem). The buildSystem14 15
Figure 29. A diagram showing the key function calls in showerQED for emitting a photon with matrix-element corrections. Arrows coming out of the bottom of a block show a function called by the block.
In the dipole configuration, the buildSystem function aggregates all charges into a single vector of vectors, grouped according to their absolute value and sign. If the charges correspond to the initial state, their signs are inverted (reflecting their effective role in the antenna picture, cf. Eq. (13) ). Subsequently, the Hungarian algorithm is applied, which guarantees a matching between charges of opposite sign and equal magnitude. These pairs are then stored in a vector that represents the complete collection of dipoles. In the multipole case16 buildSystem function constructs a matrix of dipoles that enumerates all possible charge combinations.
As a next step, the shower evolution (showerQED) is initiated at a scale set by the mass of the parent particle. The algorithm now enters a loop which starts with VinciaQED::q2Next choosing a scale for the next emission17 [58 , 59 , 60 ] . In the multipole case, the loop instead runs over the elements of the matrix of dipoles, which is not endowed with any particular ordering.
For each antenna, the algorithm generates a trial emission and selects as the physical emission the trial with the largest transverse momentum p ⊥ p_\perp p ⊥ QEDemitSystem::hadDecTrial step, which checks the availability of MECs for the current system if requested18 QEDemitSystem::acceptTrial accepts or rejects the trial emission based on the antenna function value and, in the multipole case, a sector check that avoids double counting. When the ME correction mode is enabled, the splitting kernel is constructed using the matrix elements provided by the MadGraph plugins. Otherwise, the standard antenna function is employed. In the ME correction mode, the uncorrected antenna function is used for the phase-space sectorisation, and the subsequent accept/reject step for photon emission is performed on the basis of the exact matrix elements. If the event is accepted, a photon is generated and the event record is updated accordingly (QEDemitSystem::updatePartonSystems and VinciaQED::update).
In the next iteration, the p ⊥ 2 p_\perp^2 p ⊥ 2 p ⊥ 2 p_\perp^2 p ⊥ 2
Matrix elements for corrections
To implement the ME corrections (see § 3.3 ), we specify the complete set of particles and encode their interactions in terms of Lagrangian densities in Mathematica syntax [76 ] . In this way, the particle content and the structure of the interaction terms are fully defined. This model definition is then passed to FeynRules [42 , 43 ] . Within FeynRules we formulate an effective field theory (EFT) in terms of generic hadronic degrees of freedom, characterised primarily by their electric charges and spins. FeynRules expands the EFT Lagrangian densities, systematically enumerates all admissible vertex topologies implied by the interaction terms, and writes them into a universal UFO file [42 ] .
The resulting UFO model is imported into, together with a set of run cards that specify the physical processes to be studied (including, for example, the decay chain structure and QED perturbative order). MadGraph [44 , 45 ] is used to construct all combinations of interaction vertices consistent with the specified process, generating all contributing Feynman diagrams. For each diagram, MadGraph produces numerical routines corresponding to the given amplitude, which can be accessed from within Vincia via a dedicated run-time plugin interface to MadGraph libraries.
For each specified process, MadGraph generates a dedicated C++ plugin. This plugin contains routines that iterate over all individual amplitudes, compute them, and assemble the total squared matrix element. The plugin output is used as input to Vincia, where both the pre- and post-branching matrix elements are evaluated. Their ratio is the ME-corrected antenna function.
In our setup, the antenna functions are computed using tree-level matrix elements obtained from the combined FeynRules + MadGraph framework. These matrix elements are then employed to define the antenna radiation pattern, in place of traditional dipole or multipole approximations. The full suite of generic matrix elements currently included for corrections in Vincia is listed below. For all modes, tree-level matrix elements for up to two-photon emissions are included. We list (pseudo)scalar particles as S S S V V V B + → D + K − π + B^+ \to D^+ K^- \pi^+ B + → D + K − π + D D D S 2 ± S_2^{\pm} S 2 ± D + → π + π + π − D^+ \to\pi^+\pi^+\pi^- D + → π + π + π − D D D S 1 ± S_1^{\pm} S 1 ± B 0 → D + π − B^0 \to D^+\pi^- B 0 → D + π − B 0 → K + π − B^0 \to K^+\pi^- B 0 → K + π −
Neutral vector to leptons:
Neutral scalar to leptons:
Charged scalar to leptons and (anti-)neutrinos:
S 1 + → e + ν e S_1^+ \to e^+\nu_e S 1 + → e + ν e
S 1 + → μ + ν μ S_1^+ \to\mu^+\nu_\mu S 1 + → μ + ν μ
S 1 − → e − ν ˉ e S_1^- \to e^-\bar{\nu}_e S 1 − → e − ν ˉ e
S 1 − → μ − ν ˉ μ S_1^- \to\mu^-\bar{\nu}_\mu S 1 − → μ − ν ˉ μ
Neutral scalar to two charged scalars:
S 1 0 → S 1 + S 1 − S_1^0 \to S_1^+S_1^- S 1 0 → S 1 + S 1 − B 0 → K + K − B^0 \to K^+K^- B 0 → K + K − D 0 → π + π − D^0 \to\pi^+\pi^- D 0 → π + π − K S 0 → π + π − K_S^0 \to\pi^+\pi^- K S 0 → π + π −
S 1 0 → S 1 + S 2 − S_1^0 \to S_1^+S_2^- S 1 0 → S 1 + S 2 − D ˉ 0 → K + π − \bar{D}^0 \to K^+\pi^- D ˉ 0 → K + π − B ˉ 0 → D + K − \bar{B}^0 \to D^+K^- B ˉ 0 → D + K −
Neutral scalar to two charged scalars and a neutral scalar:
S 1 0 → S 2 0 S 1 + S 1 − S_1^0 \to S_2^0 S_1^+ S_1^- S 1 0 → S 2 0 S 1 + S 1 − B 0 → D ‾ 0 K + π − B^0 \to\overline{D}^0 K^+\pi^- B 0 → D 0 K + π − D ‾ 0 → K ‾ 0 K + π − \overline{D}^0 \to\overline{K}^0 K^+\pi^- D 0 → K 0 K + π −
S 1 0 → S 2 0 S 1 + S 1 − S_1^0 \to S_2^0 S_1^+ S_1^- S 1 0 → S 2 0 S 1 + S 1 − B 0 → D ‾ 0 π + π − B^0 \to\overline{D}^0\pi^+\pi^- B 0 → D 0 π + π − B 0 → π 0 π + π − B^0 \to\pi^0\pi^+\pi^- B 0 → π 0 π + π − K L 0 → π 0 π + π − K_L^0 \to\pi^0\pi^+\pi^- K L 0 → π 0 π + π −
Charged scalar to a charged scalar and a neutral scalar:
Charged scalar to three charged scalars:
S 1 ± → S 2 ± S 2 ± S 2 ∓ S_1^\pm\to S_2^\pm S_2^\pm S_2^\mp S 1 ± → S 2 ± S 2 ± S 2 ∓ D ± → π ± π ± π ∓ D^\pm\to\pi^\pm\pi^\pm\pi^\mp D ± → π ± π ± π ∓
S 1 ± → S 3 ± S 2 ± S 2 ∓ S_1^\pm\to S_3^\pm S_2^\pm S_2^\mp S 1 ± → S 3 ± S 2 ± S 2 ∓ D ± → K ± π ± π ∓ D^\pm\to K^\pm\pi^\pm\pi^\mp D ± → K ± π ± π ∓
S 1 ± → S 4 ∓ S 3 ± S 2 ± S_1^\pm\to S_4^\mp S_3^\pm S_2^\pm S 1 ± → S 4 ∓ S 3 ± S 2 ± B ± → D ∓ K ± π ± B^\pm\to D^\mp K^\pm\pi^\pm B ± → D ∓ K ± π ±
Charged scalar to one charged scalar and two neutral scalars:
Deadcone derivations
We provide complete derivations of the expressions used in sec. 4: upper-bound on y ∗ y^* y ∗ (17) and eq. eq. (18) . We start with the definition of the massive antenna function and set it to zero:
s i k s i γ s γ k − m I 2 s i γ 2 − m K 2 s γ k 2 = 0. \frac{s_{ik}}{s_{i\gamma}s_{\gamma k}}-\frac{m_I^2}{s_{i\gamma}^2}-\frac{m_K^2}{s_{\gamma k}^2}=0 . s iγ s γ k s ik − s iγ 2 m I 2 − s γ k 2 m K 2 = 0. We define r = s i k s i γ r=\frac{s_{ik}}{s_{i\gamma}} r = s iγ s ik
s i k s i γ s γ k − m I 2 s i γ 2 − m K 2 s γ k 2 = 0 , \frac{s_{ik}}{s_{i\gamma}s_{\gamma k}}-\frac{m_I^2}{s_{i\gamma}^2}-\frac{m_K^2}{s_{\gamma k}^2}=0 , s iγ s γ k s ik − s iγ 2 m I 2 − s γ k 2 m K 2 = 0 , s i k s γ k s i γ − m I 2 s γ k 2 s i γ 2 − m K 2 = 0 , \frac{s_{ik}s_{\gamma k}}{s_{i\gamma}}-\frac{m_I^2s_{\gamma k}^2}{s_{i\gamma}^2}-m_K^2=0 , s iγ s ik s γ k − s iγ 2 m I 2 s γ k 2 − m K 2 = 0 , s i k r − r 2 m I 2 − m K 2 = 0. s_{ik}r-r^2m_I^2-m_K^2=0 . s ik r − r 2 m I 2 − m K 2 = 0. Here, we can solve for r r r
r = s i k ± s i k 2 − 4 m I 2 m K 2 2 m I 2 . r=\frac{s_{ik}\pm\sqrt{s_{ik}^2-4m_I^2m_K^2}}{2m_I^2}. r = 2 m I 2 s ik ± s ik 2 − 4 m I 2 m K 2 . Using the relationship s I K = s i γ + s γ k + s i k s_{IK}=s_{i\gamma}+s_{\gamma k}+s_{ik} s I K = s iγ + s γ k + s ik
r = s I K − s i γ − s γ k ± ( s I K − s i γ − s γ k ) 2 − 4 m I 2 m K 2 2 m I 2 . r=\frac{s_{IK}-s_{i\gamma}-s_{\gamma k}\pm\sqrt{(s_{IK}-s_{i\gamma}-s_{\gamma k})^2-4m_I^2m_K^2}}{2m_I^2}. r = 2 m I 2 s I K − s iγ − s γ k ± ( s I K − s iγ − s γ k ) 2 − 4 m I 2 m K 2 . Using the relationships E j = s i γ + s γ k + 2 m γ 2 2 m I K E_j=\frac{s_{i\gamma}+s_{\gamma k}+2m_\gamma^2}{2m_{IK}} E j = 2 m I K s iγ + s γ k + 2 m γ 2 m I K 2 = s I K + m I 2 + m K 2 m_{IK}^2=s_{IK}+m_I^2+m_K^2 m I K 2 = s I K + m I 2 + m K 2
r = s I K − s i γ − s γ k ± ( s I K − s i γ − s γ k ) 2 − 4 m I 2 m K 2 2 m I 2 , = s I K − 2 E γ m I K ± ( s I K − 2 E γ m I K ) 2 − 4 m I 2 m K 2 2 m I 2 , = m I K 2 ( 1 − x γ ) − m I 2 − m K 2 ± ( m I K 2 ( 1 − x γ ) − m I 2 − m K 2 ) 2 − 4 m I 2 m K 2 2 m I 2 . \begin{aligned}
r&=\frac{s_{IK}-s_{i\gamma}-s_{\gamma k}\pm\sqrt{(s_{IK}-s_{i\gamma}-s_{\gamma k})^2-4m_I^2m_K^2}}{2m_I^2},\\
&=\frac{s_{IK}-2E_\gamma m_{IK}\pm\sqrt{(s_{IK}-2E_\gamma m_{IK})^2-4m_I^2m_K^2}}{2m_I^2},\\
&=\frac{m_{IK}^2(1-x_\gamma)-m_I^2-m_K^2\pm\sqrt{\left(m_{IK}^2(1-x_\gamma)-m_I^2-m_K^2\right)^2-4m_I^2m_K^2}}{2m_I^2}.
\end{aligned} r = 2 m I 2 s I K − s iγ − s γ k ± ( s I K − s iγ − s γ k ) 2 − 4 m I 2 m K 2 , = 2 m I 2 s I K − 2 E γ m I K ± ( s I K − 2 E γ m I K ) 2 − 4 m I 2 m K 2 , = 2 m I 2 m I K 2 ( 1 − x γ ) − m I 2 − m K 2 ± ( m I K 2 ( 1 − x γ ) − m I 2 − m K 2 ) 2 − 4 m I 2 m K 2 . For the sake of readability we will define α = m I K 2 ( 1 − x γ ) − m I 2 − m K 2 \alpha= m_{IK}^{2}(1-x_{\gamma}) - m_{I}^{2} - m_{K}^{2} α = m I K 2 ( 1 − x γ ) − m I 2 − m K 2
r = α ± α 2 − 4 m I 2 m K 2 2 m I 2 , = α ± α 1 − 4 m I 2 m K 2 α 2 2 m I 2 . \begin{aligned}
r = \frac{\alpha\pm\sqrt{\alpha^{2} - 4m_{I}^{2}m_{K}^{2}}}{2m_{I}^{2}}, \\
&= \frac{\alpha\pm\alpha\sqrt{1-\frac{4m_{I}^{2}m_{K}^{2}}{\alpha^{2}}}}{2m_{I}^{2}}.
\end{aligned} r = 2 m I 2 α ± α 2 − 4 m I 2 m K 2 , = 2 m I 2 α ± α 1 − α 2 4 m I 2 m K 2 . We proceed taking the the logarithm of both sides and separating out common factors,
ln ( r ) = ln ( α ± α 1 − 4 m I 2 m K 2 α 2 2 m I 2 ) , = ln ( α m I 2 ) + ln ( 1 ± 1 − 4 m I 2 m K 2 α 2 ) − ln ( 2 ) . \begin{aligned}
\ln(r) &= \ln\left(\frac{\alpha\pm\alpha\sqrt{1-\frac{4m_{I}^{2}m_{K}^{2}}{\alpha^{2}}}}{2m_{I}^{2}}\right),\\
&= \ln\left(\frac{\alpha}{m_{I}^{2}}\right) + \ln\left(1\pm\sqrt{1-\frac{4m_{I}^{2}m_{K}^{2}}{\alpha^{2}}}\right) - \ln(2).
\end{aligned} ln ( r ) = ln 2 m I 2 α ± α 1 − α 2 4 m I 2 m K 2 , = ln ( m I 2 α ) + ln ( 1 ± 1 − α 2 4 m I 2 m K 2 ) − ln ( 2 ) . Now we substitute α \alpha α (18) .
y ∗ ≤ 1 2 ln ( m I K 2 ( 1 − x γ ) − m I 2 − m K 2 m I 2 ) + 1 2 ln ( 1 + 1 − 4 m I 2 m K 2 ( m I K 2 ( 1 − x γ ) − m I 2 − m K 2 ) 2 ) − 1 2 ln ( 2 ) . \begin{aligned}
y^{*} &\leq\frac{1}{2}\ln\left(\frac{m_{IK}^{2}(1-x_{\gamma})-m_{I}^{2}-m_{K}^{2}}{m_{I}^{2}}\right) +\\
&\quad\frac{1}{2}\ln\left(1+\sqrt{1-\frac{4m_{I}^{2}m_{K}^{2}}{\left(m_{IK}^{2}(1-x_{\gamma})-m_{I}^{2}-m_{K}^{2}\right)^{2}}}\right)-\frac{1}{2}\ln(2).
\end{aligned} y ∗ ≤ 2 1 ln ( m I 2 m I K 2 ( 1 − x γ ) − m I 2 − m K 2 ) + 2 1 ln ( 1 + 1 − ( m I K 2 ( 1 − x γ ) − m I 2 − m K 2 ) 2 4 m I 2 m K 2 ) − 2 1 ln ( 2 ) . From this eq. (17) follows clearly from taking m K = 0 m_{K}=0 m K = 0
y ∗ ≤ 1 2 ln ( m I K 2 ( 1 − x γ ) − m I 2 m I 2 ) . y^{*} \leq\frac{1}{2}\ln\left(\frac{m_{IK}^{2}(1-x_{\gamma})-m_{I}^{2}}{m_{I}^{2}}\right). y ∗ ≤ 2 1 ln ( m I 2 m I K 2 ( 1 − x γ ) − m I 2 ) . The definition of y ∗ y^{*} y ∗ (20) is fully laid out as
y ∗ = 1 2 ln ( s j k s i j ) , = 1 2 ln ( 2 E j E k − E j ∣ p k ∣ cos θ j k 2 E j E i − E j ∣ p i ∣ cos θ i j ) , = 1 2 ln ( E k ( 1 − β k cos θ j k ) E i ( 1 − β i cos θ i j ) ) , = 1 2 ln ( E k E i ) + 1 2 ln ( 1 − β k cos θ j k 1 − β i cos θ i j ) . \begin{aligned}
y^{*} &= \frac{1}{2}\ln\left(\frac{s_{jk}}{s_{ij}}\right),\\
&= \frac{1}{2}\ln\left(\frac{2E_{j}E_{k}-E_{j}|\boldsymbol{p}_{k}|\cos\theta_{jk}}{2E_{j}E_{i}-E_{j}|\boldsymbol{p}_{i}|\cos\theta_{ij}}\right),\\
&= \frac{1}{2}\ln\left(\frac{E_{k}(1-\beta_{k}\cos\theta_{jk})}{E_{i}(1-\beta_{i}\cos\theta_{ij})}\right),\\
&= \frac{1}{2}\ln\left(\frac{E_{k}}{E_{i}}\right)+\frac{1}{2}\ln\left(\frac{1-\beta_{k}\cos\theta_{jk}}{1-\beta_{i}\cos\theta_{ij}}\right).
\end{aligned} y ∗ = 2 1 ln ( s ij s j k ) , = 2 1 ln ( 2 E j E i − E j ∣ p i ∣ cos θ ij 2 E j E k − E j ∣ p k ∣ cos θ j k ) , = 2 1 ln ( E i ( 1 − β i cos θ ij ) E k ( 1 − β k cos θ j k ) ) , = 2 1 ln ( E i E k ) + 2 1 ln ( 1 − β i cos θ ij 1 − β k cos θ j k ) . Additional studies on performance
Figs. 30–26d show the speed-up and relative RSS memory increase for the J / ψ → μ + μ − J/\psi\to\mu^{+}\mu^{-} J / ψ → μ + μ − B 0 → K + π − B^{0} \to K^{+}\pi^{-} B 0 → K + π − D + → K + π + π − D^{+} \to K^{+}\pi^{+}\pi^{-} D + → K + π + π − § 8 ).
Figure 30. Speed-up curves for J / ψ → μ + μ − J/\psi\to\mu^{+}\mu^{-} J / ψ → μ + μ −
Figure 31. Speed-up curves for B 0 → K + π − B^0 \to K^+\pi^- B 0 → K + π −
Figure 32. Speed-up curves for D + → K + π + π − D^{+} \to K^{+}\pi^{+}\pi^{-} D + → K + π + π −
Figure 33. Maximum RSS memory consumption relative to single-thread execution for J / ψ → μ + μ − J/\psi\to\mu^{+}\mu^{-} J / ψ → μ + μ −
Figure 34. Maximum RSS memory consumption relative to single-thread execution for B 0 → K + π − B^0 \to K^+\pi^- B 0 → K + π −
Figure 35. Maximum RSS memory consumption relative to single-thread execution for D + → K + π + π − D^+ \to K^+\pi^+\pi^- D + → K + π + π −
