A new approach to QCD final-state evolution in processes with massive partons

We present an algorithm for massive parton evolution which is based on the differentially accurate simulation of soft-gluon radiation by means of a non-trivial azimuthal angle dependence of the splitting functions. The kinematics mapping is chosen such as to to reflect the symmetry of the final state in soft-gluon radiation and collinear splitting processes. We compute the counterterms needed for a fully differential NLO matching and discuss the analytic structure of the parton shower in the NLL limit. We implement the new algorithm in the numerical code Alaric and present a first comparison to experimental data.


I. INTRODUCTION
The production and evolution of massive partons are an important aspect of collider physics, and they play a particularly prominent role at the Large Hadron Collider at CERN.Key measurements and searches, such as t tH and double Higgs boson production, involve final states with many b-jets.The success of the LHC physics program therefore depends crucially on the modeling of heavy quark processes in the Monte-Carlo event generators used to link theory and experiment.With the high-luminosity phase of the LHC approaching fast, it is important to increase the precision of these tools in simulations involving massive partons.
Heavy quark and heavy-quark associated processes have been investigated in great detail, both from the perspective of fixed-order perturbative QCD and using resummation, see for example [1][2][3][4].Various proposals were made for the fully differential simulation in the context of particle-level Monte Carlo event generators [5][6][7][8].Recently, a new scheme was devised for including the evolution of massive quarks in the initial state of hadron-hadron and leptonhadron collisions [9].In this manuscript, we will introduce an algorithm for the final-state evolution and matching in heavy-quark processes, inspired by the recently proposed parton-shower model Alaric [10].The soft components of the splitting functions are derived from the massive eikonal and are matched to the quasi-collinear limit using a partial fractioning technique.In contrast to the method of [7,11], we partial fraction the complete soft eikonal, leading to strictly positive splitting functions and thus keeping the numerical efficiency of the Monte-Carlo algorithm at a maximum.We also propose to use a kinematic mapping for the collinear splitting of gluons into quarks that treats the outgoing particles democratically.This algorithm can be extended to any purely collinear splitting (i.e., after subtracting any soft enhanced part of the splitting functions) while retaining the NLL precision of the evolution.Our algorithm does not account for the effect of spin correlations.
Multi-jet merging and matching of parton-shower simulations to NLO calculations in the context of heavy-quark production were discussed, for example, in [12][13][14][15].The NLO matching is typically fairly involved, because of the complex structure and partly ambiguous definition of the infrared counterterms.In this publication, we compute the integrated counterterms for our new parton-shower model, making use of recent results for angular integrals in dimensional regularization [16].This calculation provides the remaining counterterms needed for the matching of the Alaric parton-shower model at NLO QCD.We will discuss the extension to initial-state radiation in a future publication.
This paper is structured as follows: Section II briefly reviews the construction of the Alaric parton-shower model and generalizes the discussion to massive particles.Section III introduces the different kinematics mappings.Section IV discusses the general form of the phase-space factorization and provides explicit results for processes with soft radiation and collinear splitting.The computation of integrated infrared counterterms is presented in Sec.V. Section VI discusses the impact of the kinematics mapping on sub-leading logarithms, and Sec.VII provides first numerical predictions for e + e − →hadrons.Section VIII contains an outlook.

II. SOFT-COLLINEAR MATCHING
We start the discussion by revisiting the singularity structure of n-parton QCD amplitudes in the infrared limits.If two partons, i and j, become quasi-collinear, the squared amplitude factorizes as n ⟨1, . . ., n|1, . . ., arXiv:2307.00728v2 [hep-ph] 3 Oct 2024 where the notation i \ indicates that parton i is removed from the original amplitude, and where (ij) is the progenitor of partons i and j.The functions P λλ ′ ab (z) are the spin-dependent, massive DGLAP splitting functions, which depend on the momentum fraction z of parton i with respect to the mother parton, (ij), and on the helicities λ [11,[17][18][19][20][21].For the remainder of this work, we will use only the spin-averaged splitting functions.
In the limit that gluon j becomes soft, the squared amplitude factorizes as [22] n ⟨1, . . ., n|1, . . ., n⟩ n = −8πα s i,k̸ =i where T i and T k are the color insertion operators defined in [22,23].In the remainder of this section, we will focus on the eikonal factor, w ik,j , for massive partons and how it can be re-written in a suitable form to match the spin-averaged splitting functions in the soft-collinear limit.The eikonal factor is given by Following Refs.[10,24], we split Eq. ( 3) into an angular radiator function, and the gluon energy , where The parton velocity is defined as v = 1 − m 2 /E 2 and is frame dependent.When we label velocities by particle indices in the following, it is implicit that they are computed in the frame of a time-like momentum n µ .In this reference frame they reduce to the relative velocity v i = v pin , where For convenience, we also introduce the velocity four-vector When this vector is labeled by a single particle index only, it again refers to the four-velocity of that particle in the frame of n µ , for example l µ i = l µ pi,n .When partial fractioning Eq. ( 4), we aim at a positive definite result in order to maintain the interpretation of a probability density and, with it, the high efficiency of the unweighted event generation in a parton shower.Following the approach of Ref. [10], we obtain The partial fractioned angular radiator function can be written in a more convenient form using the velocity fourvectors.We find an expression that makes the matching to the ij-collinear sector manifest In the quasi-collinear limit for partons i and j, we can write the eikonal factor in Eq. ( 3) as [11,21] w ik,j i||j This can be identified with the leading term (in 1 − z) of the DGLAP splitting functions P aa (z, ε), where1 To match the soft to the collinear splitting functions, we therefore replace where the two contributions to the gluon splitting function are treated as two different radiators [25].As in the massless case, this substitution introduces a dependence on a color spectator, k, whose momentum defines a direction independent of the direction of the collinear splitting [10].In general, this implies that the splitting functions which were formerly dependent only on a momentum fraction along this direction, now acquire a dependence on the remaining two phase-space variables of the new parton.It is convenient to define the purely collinear remainder functions They can be implemented in the parton shower using collinear kinematics, while maintaining the overall NLL precision of the simulation.This will be discussed further in Sec.VI.

III. MOMENTUM MAPPING
Every parton shower algorithm requires a method to map the momenta of the Born process to a kinematical configuration after additional parton emission.This mapping is linked to the factorization of the differential phasespace element for a multi-parton configuration.Collinear safety is a basic requirement for every momentum mapping.In addition, a mapping is NLL safe if it preserves the topological features of radiation simulated previously [26,27].We will make this statement more precise in Sec.VI.
This section provides a generic momentum mapping for massive partons that is both collinear and NLL-safe.It is constructed for both initial-and final-state radiation, inspired by the identified particle dipole subtraction algorithm of [23].In the massless limit, we reproduce this algorithm and thus agree with the existing parton-shower model Alaric [10].

A. Radiation kinematics
We will first describe the kinematics mapping needed for soft evolution.This is sketched in Fig. 1.The momenta K and pij are mapped to the momenta K, p i and p j , while pk acts as a spectator that defines the azimuthal direction.
We assume that any of the particles i, j and k can be massive.The momentum K can be chosen in a suitable way to reflect the dynamics of the process [10] and absorbs the recoil in the splitting.We define the variables and the analogous final-state variables µ 2 i/j = m 2 i/j /(2p ij K) and μ2 i/j = 2µ 2 i/j /(1 + v pij K ).The scalar invariants after the splitting are defined in terms of the energy fraction, z [23] The momentum of the radiator after the emission is where 1. Sketch of the radiation kinematics in final-state evolution.See the main text for details.Note that p k is unaltered by the mapping and only acts as a reference for the azimuthal angle ϕ (cf.Eqs. ( 21) and ( 22)).
We define the variable and where v = (p i p j )/(p i n) is defined as in the massless case [10,23].In terms of these quantities, the gluon momentum is given by where In order to determine a reference direction for the azimuthal angle ϕ = arctan(k y /k x ), we note that the soft radiation pattern must be correctly generated.To achieve this, we compose the transverse momentum as where the reference axes n ⊥ and l ⊥ are given by the transverse projections 2 and Since the differential emission phase-space element is a Lorentz-invariant quantity, the azimuthal angle ϕ is Lorentz invariant [10].This allows us to write the emission phase-space in a frame-independent way.Upon construction of the momenta p i , p j and K, the momenta {p l } which are used to define K are subjected to a Lorentz transformation, 2 In kinematical configurations where p µ k is a linear combination of p µ i and n µ , n ⊥ in the definition of Eq. ( 20) vanishes.It can then be computed using n ⊥ = ε µj νρ p ν i n ρ , where j ∈ {1, 2, 3} may be any index that yields a nonzero result.Note that in this case, the matrix element cannot depend on the azimuthal angle.
2. Sketch of the splitting kinematics in final-state evolution.See the main text for details.Note that the unlabeled momentum is unaltered by the mapping, and only acts as a reference direction to define the azimuthal angle ϕ.

B. Splitting kinematics
In this section, we describe the kinematics for the implementation of the purely collinear components of the splitting functions in Eq. ( 12).This is sketched in Fig. 2. We make use of some of the notation in [11], in particular and we define the scaled masses We also introduce the following variables In terms of the additional variables the momenta after the splitting are given by The transverse momentum squared is given by The construction of the transverse momentum vector proceeds as described in Sec.III A. While the choice of the reference vector defining n µ ⊥ is in principle arbitrary, it can be made conveniently, e.g., to aid the implementation of spin correlations in collinear gluon splittings.

IV. PHASE SPACE FACTORIZATION
In this section, we discuss the factorization of the differential n + 1 particle phase-space element into a differential n particle phase-space element and the radiative phase space.We start from the generic four-dimensional expression Replacing particles i and j with a combined "mother" particle ıȷ, we can write the expression for the underlying Born phase space as dΦ n−1 (p a , pb ; p1 , . . ., pij , . . ., \ pj , . . ., pn ) = The ratio of differential phase-space elements after and before the mapping is defined as the differential phase-space element for the one-particle emission process This expression naturally generalizes to D dimensions.It can be computed using the lowest final-state multiplicity, i.e., n = 3.In order to do so, we first derive a suitable expression for the D-dimensional two-particle phase space in the frame of an arbitrary, time-like momentum n.It reads where dΩ p,n is the integral over the D − 1 dimensional solid angle of p µ in the frame of n µ .All energies and momenta are computed in this frame, and we have defined P = p + q.We can re-write the momentum-dependent denominator factor on the right-hand side of Eq. ( 33) as Inserting this into Eq.( 33), and using D = 4 − 2ε, we obtain a manifestly covariant form of the phase-space element, dΦ 2 (P ; p, q) = (4π

A. Radiation kinematics
To derive the emission phase space for final-state splittings in the radiation kinematics of Sec.III A, we make use of the standard factorization formula where q = k̸ =i,j p k is the sum of all final-state momenta except p i and p j .We can use Eq. ( 35) to relate the phase-space factor dΦ 2 (−n; p i , q) to the underlying Born phase space as follows where we have defined The differential two-body phase-space element for the production of n µ and p µ j is given by Finally, we combine Eqs. ( 37) and ( 39) to obtain the single-emission phase space element in Eq. (32) dΦ +1 (p a , pb ; . . ., pij , . . .
We demonstrate in App.A 1 that this expression agrees with the result derived in Ref. [23].

B. Splitting kinematics
To derive the emission phase space for final-state radiation in the splitting kinematics of Sec.III B, the recoiler is chosen to be the sum of the remaining final-state partons.We make use of the standard factorization formula dΦ 3 (q; p i , p j , K) = dΦ 2 (q; K, p ij ) dp where q = k̸ =i,j p k is the sum of final-state momenta except p i and p j .Working in the frame of q = K + pij , we can use Eq. ( 35) to relate the phase-space factor dΦ 2 (q; K, p ij ) to the underlying Born differential phase space element dΦ 2 (q; K, pij ) as follows The decay of the intermediate off-shell parton is associated with the differential phase-space element Finally, we combine Eqs. ( 43) and ( 44) to obtain dΦ +1 (p a , pb ; p1 , . . ., pij , . . ., pn ; In the massless limit, this simplifies to dΦ +1 (p a , pb ; p1 , . . ., pij , . . ., pn ; In App.A 2, we will show the equivalence of Eq. ( 45) to the single-emission differential phase-space element of [11].

V. INFRARED SUBTRACTION TERMS AT NEXT-TO-LEADING ORDER
The matching of parton showers to fixed-order NLO calculations in dimensional regularization based on the MC@NLO algorithm [28] requires the knowledge of integrated splitting functions in D = 4 − 2ε dimensions.Since our technique for massive parton evolution is modeled on the Catani-Seymour identified particle subtraction and the Alaric parton-shower, we can use the methods developed in [29,30].We will limit the discussion to the main changes needed to implement the algorithm for massive partons.The results of this section provide an extension of the subtraction method for identified hadrons first introduced in [23].

A. Soft angular integrals
In this section, we compute the angular integrals of the partial fractioned soft eikonal.While we focus on pure final-state radiation, the results of this integration are generic and apply to the case of final-state and initial-state emission of massless vector bosons.The integrand is given by Eq. ( 8), We can write the angular phase space for the emission of the massless particle j as where Ω(n) = 2π n/2 /Γ(n/2) is the d-dimensional line element, in particular Ω(1 We finally find the following expression for the cases with massive emitter For massless emitter, we obtain (see also [10,23]) The angular integrals I 1,1 and I 2 have been computed in [16,[31][32][33][34].They read [16] where the velocities are defined in the notation of Ref. [16], and where Note that spurious collinear poles appear on the right-hand side of Eq. ( 50), which cancel between I (1) 1,1 (l i l ik /2, l 2 ik /4)/2 and I (1) 1,1 (l i l k , l 2 k ).In addition, in the fully massless case 1,1 is related to the massless two-denominator integral and gives the simple result [10]

B. Soft energy integrals
In this section, we introduce the basic techniques to solve the energy integrals in Eq. ( 40).After performing the angular integrals, we are left with the additional z-dependence induced by the energy denominator in Eq. ( 4).We focus on the cases relevant for QCD and QED soft radiation, where µ ij = µ i and µ j = 0.The differential emission probability per dipole is then given by dΦ +1 (p a , pb ; . . ., pij , . . .
In order to carry out the integration over the energy fraction, z, we expand the integrand into a Laurent series.The differential soft subtraction counterterm, summed over all dipoles, is given by dσ S S = − 8πα s µ 2ε ıȷ,k dΦ +1 (p a , pb ; . . ., pij , . . . where and where the mass correction factor reads and where the massless limit is given by J ik (z, 0, κ) = 1.All 1/ε poles have been extracted, and we can (after expanding the ε-dependent prefactors) compute the final result as an integral over the delta functions and plus distributions in z.In general, some of the terms must be computed numerically, as n implicitly depends on z, see Eq. ( 14).In order to apply Eq. ( 56) to processes with resolved hadrons, we can make use of the formalism derived in [23,29,30].The 1/ε pole proportional to 2z/[1 − z] + can then be combined with the soft enhanced part of the collinear mass factorization counterterms.The extension to initial-state radiation requires a repetition of the derivation in Sec.IV A for initial-state kinematics.We will discuss this in a future publication.

C. Collinear integrals
To compute the purely collinear counterterms, we use the splitting kinematics and make use of the results in App.A 2. Note that we also use the corresponding definitions of the scaled masses.We define the purely collinear anomalous dimension in terms of the collinear remainders in Eq. ( 12) where z ± are given by Eq. (A9).This leads to The complete counterterm can now be obtained by Laurent series expansion, using Eq.(A14) There are three variants of Eq. ( 59) relevant for QCD.The first describes the collinear splitting of a massive quark into a quark and a gluon, and is characterized by μi = μij > 0 and μj = 0. We obtain dΦ +1 (q; pij , K; Note that the result is infrared finite.This is expected because the only poles that occur in the radiation off massive quarks have their origin in soft gluon emission, which is fully accounted for by the eikonal integral in Eq. ( 54).This also shows that in our subtraction scheme, there is a hierarchy between the soft and the collinear enhancements, as required for a proper classification of leading and sub-leading logarithms.
The second case is the splitting of a gluon into two massive quarks, which is characterized by μi = μj > 0 and μij = 0.The result is finite and reads dΦ +1 (q; pij , K; The final case describes the collinear decay of a massless parton into two massless partons, and is characterized by μi = μj = μij = 0.This term is collinearly divergent, and we obtain dΦ +1 (q; pij , K; The treatment of initial-state singularities will be discussed in a future publication.

VI. KINEMATICAL EFFECTS ON SUB-LEADING LOGARITHMIC CORRECTIONS
We now demonstrate that our kinematics mapping satisfies the criteria for NLL accuracy laid out in [26,27], extending the discussion of the massless case in [10].We refer the reader to [10] for more details on the general method of the proof.The analysis is based on the technique for general final-state resummation introduced in [35].Following the notation of Sec. 2 of [35], the observable to be resummed is denoted by v, while the hard partons and soft emission have momenta p 1 , . . ., p n and k, respectively.The observable is a function of these momenta, v = V ({p}, {k}), and for any radiating color dipole formed by the hard momenta, p i and p j , the momentum of a single emission can be parametrized as where with rapidity η ij = 1/2 ln(z i,j /z j,i ).An observable can be expressed as with k T,l = k T,lj and η l = η lj for j ∦ l in the collinear limit.A natural extension of Eq. ( 63) to the case of massive emitters would be In the quasi-collinear limit, this Sudakov decomposition agrees with the one given in [21,36] if the auxiliary (lightlike) vector defining the anti-collinear direction is chosen as the direction of the color partner of the QCD dipole.In particular, for constant z i,j , z j,i and small k T , the gluon momentum behaves as in complete agreement with Eq. ( 63).In the quasi-collinear limit, the value of the observable V (k) will therefore be unchanged from the case of massless evolution.However, both the Sudakov radiator and the F function will change, due to the modified splitting functions, Eq. ( 10), and the effects of masses on the integration boundaries.The Lund plane for gluon emission off a dipole containing a massive quark with mass m and energy E will have a smooth upper rapidity bound at η = ln(E/m), consistent with the well known dead-cone effect [37][38][39].When the similarity transformation introduced in Sec.2.2.3 of [35] is generalized to the quasi-collinear limit, and applied to a process with massive quarks, the location of this boundary in relative rapidity is unchanged, because the quasicollinear limit requires m ∝ k T .The sub-leading logarithmic terms in the resummed cross section can then be extracted similarly to the massless case, but a change to the CMW scheme is required.We will not discuss the complete structure of the result, but address only the effects related to momentum mapping, which have been found to spoil NLL precision [26,27].
In order to prove that the momentum mapping of Sec.III A satisfies the criteria laid out in [26,27], we need to show that it has the same topological features as in the massless case [10].This amounts to showing that hard, (quasi-)collinear emissions, i.e. highly energetic emissions in the direction of the hard partons, do not generate Lorentz transformations that change existing momenta by more than O((k 2 T /K 2 ) ρ), where ρ is a positive exponent.To show this, we analyze the behavior of Eq. ( 23) in the (quasi-)collinear region.We can split K µ into its components along the original recoil momentum, Kµ , the emitter, pµ i , and the emission, p µ j .
For the effect of the Lorentz transformation that determines the event topology after radiation, only the parametric form of X µ is of interest.Using Eqs. ( 15) and ( 18), the small momentum X µ for gluon radiation can be written as In the quasi-collinear limit, this scales as O(k ⊥ ), which is sufficient to make the effect of the Lorentz transformation vanish even for hard quasi-collinear splittings.Finally, we note that the precise treatment of transverse recoil in the kinematics mapping of Sec.III B does not affect the resummed result at NLL precision, because the mapping is applied solely to the purely collinear parts of the splitting functions, Eq. ( 12).This can be understood with the help of Eq. (2.46) in Ref. [35], which will be structurally similar in the case of massive partons.The sub-leading logarithmic terms, which lead to a violation of NLL accuracy in many dipole showers, arise from accidental correlations between multiple soft-enhanced emissions.This means in particular, that the parton-shower equivalent of the integrals in Eq. (2.46) does not factorize in the strongly ordered limit.The differential probability for the emissions is given by the derivative of the radiator function in Eq. (2.21) of Ref. [35] and does not receive a contribution from the purely collinear remainders in Eq. ( 12), such that a change in the kinematics mapping cannot generate this particular type of NLL violation.

VII. COMPARISON TO EXPERIMENTAL DATA
In this section we present first numerical results obtained with the extension of the Alaric final-state parton shower to massive parton evolution.The implementation is part of the event generation framework Sherpa [46][47][48].We do not perform NLO matching or multi-jet merging, and we set C F = (N 2 c − 1)/(2N c ) = 4/3 and C A = 3.The quark masses are set to m c = 1.42 GeV and m b = 4.75 GeV and the same flavor thresholds are used for the evaluation of the strong coupling.The running coupling is evaluated at two loop accuracy, and we set α s (m z ) = 0.118.Following standard practice to improve the logarithmic accuracy of the parton shower, we employ the CMW scheme [49].In this scheme, the soft eikonal contribution to the flavor conserving splitting functions is rescaled by 1 + α s (t)/(2π)K, where K = (67/18 − π 2 /6) C A − 10/9 T R n f , and where n f is scale dependent with the same flavor thresholds as listed above.Velocity dependent corrections to K should be included in principle, but the massless result provides an acceptable approximation at very large velocities [50,51], and can therefore be used for b-quark production at LEP, where v ≈ 0.99.Our results include the simulation of hadronization using the Lund string fragmentation implemented in Pythia 6.4 [52].We use the default hadronization parameters, apart from the following: PARJ(21)=0.3,PARJ(41)=0.4,PARJ(42)=0.45,PARJ(46)=0.5 and compare our predictions with those from the Dire parton shower [25].All analyses are performed with Rivet [53].
Figure 3 displays predictions from the Alaric parton shower for differential jet rates in the Durham scheme compared to experimental results from the JADE and OPAL collaborations [40].The b-quark mass corresponds to y ∼ 2.8 • 10 −3 , and hadronization effects dominate the predictions below ∼ 10 −4 .We observe good agreement with the experimental data.Overall, the quality of the results is similar to Ref. [10], where heavy quark effects were modeled by thresholds.Figure 4 shows a comparison for event shapes measured by the ALEPH collaboration [41].It can be expected that the numerical predictions will improve upon including matrix-element corrections, or when merging the parton shower with higher-multiplicity calculations.Again, the overall quality of the prediction is similar to Ref. [10].
Finally, we show predictions for the b-quark fragmentation function as measured by the ALEPH [42], DELPHI [43], OPAL [44] and SLD [45] collaborations.Figure 5 shows a fair agreement of both the Alaric and Dire predictions with experimental data.We note that both parton shower implementations use the same hadronization tune.

VIII. CONCLUSIONS
We have introduced an extension of the recently proposed Alaric parton-shower model to the case of massive QCD evolution.An essential aspect of the new algorithm is the form of the collinearly matched massive eikonal, which is obtained by partial fractioning the angular component of the eikonal of a complete dipole.This technique preserves the positivity of the splitting function, thus leading to an excellent efficiency of the Monte-Carlo simulation.Inspired by the symmetry of the partonic final state in purely collinear splittings, we also introduced a dedicated kinematics mapping for this scenario and showed that it preserves the NLL precision of the overall simulation.We computed the infrared counterterms needed for the matching to fixed-order calculations at NLO accuracy, and discussed the logarithmic structure of the resummation in the case of heavy-quark evolution.
Several improvements of this algorithm are needed before it can be considered on par with the parton shower simulations used by past and current experiments.Clearly, spin correlations and dominant sub-leading color effects should be included.This can be achieved with the help of the techniques from [54][55][56][57][58][59][60].An extension to initial-state evolution is needed for LHC phenomenology.It will need to account for the non-cancellation of certain types of singularities in processes with two massive initial states [61].Finally, the algorithm should be extended to higher orders based on the techniques developed in [59,62,63].[42][43][44][45].
In this context, we note that the all-orders (in ε) expressions from Sec. V, in conjunction with higher-order expressions for the angular integrals in Sec.V A that can be obtained from [16], can be used to compute the factorizable integrals at NNLO, thus providing a significant part of the components needed for an MC@NNLO matching [64].The computation of the remaining non-factorizable integrals is a further development needed in order to reach the precision targets of the high-luminosity LHC. and finally dΦ 2 (p ij ; p i , p j ) = (4π We also have The last missing component of the emission phase space is the ratio of the production to the underlying Born phasespace element.It is most easily computed in the frame of q = pij + K and results in dΦ 2 (q; p ij , K) where the Källen function is given by λ(a, b, c) = (a − b − c) 2 − 4bc.Combining Eqs.(A11)-(A13) gives the final result dΦ +1 (q; pij , K; for the initial-state momenta p µ a and p µ b and final-state momenta, {p µ 1 , . . ., p µ i , . . ., p µ j , . . ., p µ n }, dΦ n (p a , p b ; p 1 , . . ., p i , . . ., p j , . . ., p n dΦ +1 (p a , pb ; p1 , . . ., pij , . . ., pn ; p i , p j ) = dΦ n (p a , p b ; p 1 , . . ., p i , . . ., p j , . . ., p n ) dΦ n−1 (p a , pb ; p1 , . . ., pij , . . ., \ pj , . . ., pn ) .
[41]ic and Dire predictions in comparison to LEP data from[40].Alaric and Dire predictions in comparison to LEP data from[41].Alaric and Dire predictions in comparison to LEP data from