Towards a fully massive five-flavour scheme

In this work we explore first necessary steps to contruct a fully massive version of a variable flavour number scheme. In particular we focus, as an example, on an extension of the five-flavour scheme, where instead of neglecting explicit initial state quark mass effects, we retain all massive dependence, while keeping the resummation properties of the massless five-flavour scheme. We name this scheme five-flavour-massive (5FMS) scheme. Apart from consistently modified parton distribution functions, we provide all the ingredients that are needed to implement this scheme at MC@NLO accuracy, in a Monte Carlo event generator. As proof of concept we implement this scheme in SHERPA, and perfom a comparison of the new scheme with traditional ones for the simple process of scalar particle production in bottom quark fusion.


Introduction
Processes with heavy quarks (bottom or charm) in the initial state present an interesting challenge for theoretical predictions at the LHC and other hadron collider experiments. First, the finite quark masses introduce another scale to the process, which may or may not prove to be relevant for different observables and different processes. In addition, a decision has to be made in how far heavy quarks can act as incident partons -due to their mass being larger than the QCD scale parameter m Q ≫ Λ QCD one could argue that they are disallowed to have a parton distribution function (PDF), thereby decoupling them from the QCD evolution in the initial state, described by the DGLAP equations. This leads to two complementary solutions. On one hand, heavy quarks Q in the initial and final state may be treated on the same footing as any other light quark, such as the u, d, or s quark, by ignoring their mass in the evaluation of the matrix elements. In such a picture the heavy quark acts as an active quark in the QCD evolution equations, and consequently, possibly large collinear logarithms are resummed to all orders into a Q PDF. On the other hand, for some processes and observables the effects of the finite heavy quark mass, m Q become relevant and in such cases these quarks must be treated as fully massive. Traditionally, this immediately translates into the heavy quarks only appearing as final state particles. This dichotomy is most pronounced for the case of b quarks, due to their mass m b ≈ 4.5 GeV being larger than the charm mass by a factor of about 3. It gives rise to ongoing comparisons of calculations of the same processes and observables in the five-flavour and four-flavour schemes. Here, the former refers to a consistently massless treatment of the b quark, which can therefore be found in both initial and final states, while the latter treats the b quarks as massive and allows them to be in the final state only. For a recent example focusing on the production of Z or Higgs bosons in association with b quarks at the LHC, cf. [1].
There, a slight preference for five-flavour scheme calculations in a multijet merging approach has been found. Broadly speaking, for a wide range of kinematical observables such as the p ⊥ spectrum of jets or gauge or Higgs bosons away from small momenta, this is in agreement with other, similar studies [2,3,4]. A preferable solution, would be to perform a matching between these two scheme, see for examples [5,6,7,8,9,10,11,12,13,14]. However, so far, these schemes have only been worked out for inclusive enough observables, and are not yet suitable for a Monte Carlo implementation. This finding motivates to extend the five-flavour scheme to allow massive particles in the initial state. In this paper, we present the necessary ingredients for next-to-leading order calculations with massive initial state partons, including these mass effects in initial state parton shower. We refer to the extended scheme as five-flavour-massive scheme, or 5FMS in short. This scheme thus has massive b-quarks that contribute both to the running of the coupling constant and to the evolution of PDFs. There are a number of obstacles to this goal: 1. in order to calculate cross sections at next-to leading order accuracy in the strong coupling, a scheme to identify, isolate and subtract infrared divergences is yet to be worked out in full detail. In particular, we follow the logic of the Catani-Seymour subtraction formalism, which was first presented for massless partons in [15], and later extended to massive fermions in QED in [16], to massive final state QCD partons in [17], and to massive initial state quarks for initial-final dipoles in [18]. We work out phasespace mappings and differential and integrated splitting kernels for the emission of a gluon off a massive quark line in the initial state, with an initial state spectator.
The treatment of massive initial state particles in QED has already been discussed in [16]. However, in contrast to the results there, obtained in D = 4 dimensions with a massive photon with m γ as infrared regulator, we consistently work in D = 4 − 2ε dimensions with a massless gluon. Of course, expressions can be mapped onto each other by suitably replacing: and, working in the MS scheme, 2. standard five-flavour PDFs introduce massive quarks purely perturbatively, through gluon splitting within the evolution. In so doing, special care is devoted to the treatment of threshold effects due to the finite masses, resulting in variable-flavour number schemes, such as the ones detailed in [5,6,7,8,10]. However, all these schemes treat mass effects only through thresholds and usually ignore other kinematical effects. We modify standard PDFs through a number of plausible choices detailed below to obtain some handle on the size of such effects. However, a full and comprehensive study of mass effects in PDFs is beyond the scope of this publication. This is also true for more conceptual questions in how far such mass effects must be treated as process-dependent corrections, similar to higher-twist effects. While we acknowledge that these may be important considerations, we leave the detailed study of these effects in hadron-hadron collision for a separate work.
An additional problem, as established in [19,20,21,22,17], is that starting at NNLO there are non-cancelling infrared divergent contributions that are proportional to the initial state quark mass. This renders the scheme presented in this paper only valid up to NLO accuracy. Lastly, approaches that use a finite heavy quark mass in the parton shower have been studied in literature [23,24,25]. Although a comparison is certainly interesting, we leave this to future studies. The outline of this work is the following. In the next section we very briefly summarize the Catani-Seymour subtraction procedure. There, we also report the ingredients needed to extend this method to include massive initial state quarks. In section 3 we present the relevant modifications for the matching of the parton shower and next-to-leading order matrix elements. We provide the discussion of results in section 4 where we show explicit results for the production of a scalar boson, A, in bottom-quark fusion for various combinations of (m A , tan β) both at fixed-order and matched to the parton shower. There we also compare our results with the DIRE parton shower, which includes (collinear) NLO corrections to the DGLAP equation [26,27]. 2 Catani-Seymour subtraction for massive initial states 2.1 Nomenclature The differential leading-order cross section (LO) for a hard scattering process with N particles in the final state is given by where B denotes the Born matrix element squared and the differential phase space element dΦ N implicitly contains the incoming flux of the incident particles, and parton distribution functions, where applicable. Later, where they matter, we will make these factors explicit. The measurement function, F (N ) J , guarantees that the N -jet final state is well defined at Born level and for Born kinematics. In a similar fashion, and suppressing the obvious four-vectors as arguments, the cross section at next-to leading order (NLO) cross section is given by where the additional terms V N and R N signify the virtual and real corrections to the original Born term. They of course relate to final states with N and N + 1 particles, respectively, as indicated by the phase space elements. The measurement function must satisfy to ensure a meaningful cross section definition at NLO accuracy. In later parts of this publication we will assume that this function is implicitly included. The soft and collinear divergences related to the emission of the additional particle in R N are cancelled by similar structures in the virtual part V N , but in order to facilitate this cancellation the poles in both must first be isolated and dealt with. Subtraction methods, such as Catani-Seymour subtraction, make use of the universal property of QCD amplitudes in the soft and collinear limits, where the corresponding divergent poles factorize. This allows the construction of process-independent subtraction terms S N (Φ N ⊗ Φ 1 ), such that the first term on the right hand side of is finite. Assuming that infrared divergences stem from regions in phase space where the momentum p k of a particle k becomes soft or collinear to a particle i with momentum p i , the degree of infrared divergence can be parametrized by a small λ → 0 such that |p k | ∼ λ or p i · p k ∼ λ. In these divergent phase space regions, the difference in the first term on the r.h.s. of the equation above behaves as i.e. all infrared poles have been cancelled. In addition, the functions S(Φ N ⊗ Φ 1 ) are constructed in such a way that their integral over the extra emission phase space -the second term on the right hand side of Eq.(2.4) -can be calculated analytically in D = 4 + 2ε dimensions, with their divergent parts giving rise to poles 1/ε 2 and 1/ε. These poles are ultimately combined with the infrared poles from the loop contributions to cancel exactly. Combining Eqs.(2.2) and (2.4) yields In Catani-Seymour subtraction, the terms S are formulated in terms of dipoles, made from three particles, an emitter, a spectator and the emitted particle k. The subtraction term factorizes into a product of processindependent emission terms and Born-like configurations, possibly with parton flavours that differ from the original Born term. These dipoles D are then classified by the splitter and spectator parton being either in the initial (I) or final (F) state, as II, FI, IF, or FF -the emitted particle k obviously always is in the final state. The overall subtraction term therefore reads where each dipole contribution is given by the sum of all possible emitter-spectator pairs, In the context of this study we primarily focus on II configurations, which can most conveniently be studied in quark-annihilation processes such as bb → H and similar. For these processes, represents one individual dipole contribution, where the emitter is the initial state particle a and the spectator is the other initial state particle b. The matrix element M emerges from the original Born-level matrix element by taking into account that the emission of k off parton a might alter the flavour of the resulting partonã, and it is evaluated at a kinematical configuration, where the modified four-vectorsp a andp i account for four-momentum conservation by absorbing the recoil from emitting p k . The functions 1 2 pi·pj V ij,k (p i , p k , p j ) are generically called splitting kernels, or dipole splitting functions and reduce to the well-known Altarelli-Parisi [28] splitting kernels P ij in the collinear limit. Note that the ⊗ symbol implies possible summation in color and helicity space.

Massive II dipoles: initial state splitter with initial state spectator
We now present the relevant modification to the described picture, due to the inclusion of finite masses. In the following we make extensive use of the following kinematical quantities: The only dipoles involving massive partons which exhibit infrared divergences are those corresponding to the emission of a gluon into the final state. Therefore, we only have to find a suitable expression for the term V ak,b (p a , p b , p k ) for the case where a heavy initial quark Q emits a gluon with momentum p k . In analogy to the treatment for massive final state particles in [17] this is given by In any other case -g → QQ and Q → gQ -there are no singular contributions rendering the need for subtraction obsolete. It is clear, however, that the collinear divergences present for the very same splittings in the massless case give rise to logarithmically enhanced terms of the form log m 2 Q /µ 2 , where µ is some large scale related to the dipole kinematics. While one may be tempted to use subtraction terms to smooth these structures and make them more amenable to numerical integration, we have tested explicitly that they do not pose any problem for processes at LHC energies and masses down to 1 GeV. In any case detailed expressions in the QED case in four space-time dimensions can be found in [29]. Coming back to the case of gluon emissions off a heavy quark line, we set the subtraction term to zero for x ab < α, the kinematical lower bound, (2.14) In a next step we need to define the phase space map, connecting the original momenta {p i } of the real emission configuration to the modified momenta {p i } for the reduced matrix element in the subtraction term. This map has to preserve mass-shell conditions, and in particularp 2 a = p 2 a = m 2 a , and it is also customary to keep the spectator momentum fixed. As a consequence of these conditions, all other final state momentap j and their total momentumQ = jp j absorb the recoil in the reduced matrix element. The transformations are given by, where the Lorentz transformation Λ µ ν is given by and applied to all final state particles apart from k, including colourless ones. It is straightforward to check that these relations fulfil the mass-shell conditions, such that p 2 a = m 2 a and Q 2 = Q 2 , and that they possess the right infra-red and collinear asymptotic limits.

Phase space
The phase space for the real emission correction factorises into a Born-level part and a one-particle phase space integral, where x-dependent momenta can be obtained from p a and Q upon replacing Q 2 → s ab x + m 2 a + m 2 b . The extra particle phase space reads where, for convenience, we define v = y a / (1 − x), and The integrated splitting function V N a,b is given by 20) and can be decomposed according into an end-point contribution V a,b N (ε) containing the 1/ε pole and a finite part The individual pieces read where Note that due to the absence of a collinear divergence -which is shielded by the finite quark mass -there are no terms ∝ 1/ε 2 . The contribution to the partonic differential cross section is thus given by where we make explicit the dependence on the is the initial state flux, φ.
To embed this into the calculation of cross section at hadron colliders, the PDFs must be added. Parametrizing the incoming hadron and parton momenta as yields the allowed intervals for the light-cone momentum fractions Making explicit the flux, φ(s ab ) = 4 λ(s, m 2 a , m 2 b ), the integral over the incoming light-cone momenta, and the parton distribution functions, the integrated splitting function, corresponding to the purely partonic expression in Eq. (2.24), reads (2.27) In this form I is not very useful for direct implementation, because the first term in the second line implies that the parton-level Born cross section must be integrated over all values of x in the interval [α, 1]. To remedy this, we need to rewrite this term, to disentangle the Born cross section from the x-integration. To fix this we define the following variable transformation, which defines a Jacobean, J(η ′ 1 (x), η 2 ), and Note that the Jacobean reduces to the usual 1/x factor in the massless limit. After reversing the integration order, and performing the change of variable, we find is the old variable expressed in terms of the new ones, and where the integration boundary is now given bȳ Renaming η 1 ↔ η ′ 1 finally yields This disentangles the evaluation of the Born cross section from the x-integral such that the whole curly bracket in Eq. (2.33) acts as a local K-factor on top of the partonic cross section.

Dipole formulae for initial-final configurations
A detailed derivation of dipole formulae in the initial-final and final-initial cases can be found in [18], although in a slightly different notation compared to the one presented in this work. In principle one could also extract all relevant formulae from [16] with the modifications described in the introduction, following the steps presented in the previous section: We consider the splitting: Q a → g k Q with spectator i in the final state. To make the reading of this section more transparent, we also report some useful kinematical quantities used throughout: x ai = p a · p i + p a · p k − p i · p k p a · p i + p a · p k ; y a = p a · p i p a · p i + p a · p k ; (2.34) The subtraction term in this case is given by D ai k (p 1 , . . . , p i , . . . , p k , . . . , p N +1 ; p a , p b ) = − where the only divergent dipole contribution in the massive case reads The mapped momenta can be expressed in terms of the original momenta using The integral of the extra emission phase-space can be split into two contributions, as done in Eq.(2.21), V a i,N (ε) and [K a i (x)] + . They are given by and . (2.41) The rest of the derivation follows exactly as in the previous section.

MC@NLO matching
Having successfully built fixed-order NLO matrix elements in the 5FMS, we now proceed to the matching to the parton shower along the lines of the well-established MC@NLO technique [30] as implemented with small variations in the SHERPA Monte Carlo [31], and referred to as S-MC@NLO [32,33,34]. Note that our implementation closely follows that of [32,33], which we refer to, for further details. We start, by constructing the NLO-weighted Born cross section, can be split into an unresolved, divergent part, and a hard, resolved part, and redefined Eq. (2.7 )S ≡ D (S) to be consistent with the notation commonly used in this context. It is worth mentioning at this stage that D (A),(S) have the same formal structure and they only differ by finite terms * . As a consequence they can both be written using the structure of Eqs. (2.7,2.8), The last ingredient needed is the MC@NLO Sudakov form factor. This is built starting from D (A) (Φ N +1 ), In particular, Eq. (2.7) implies that the Sudakov form factor can be decomposed as The inclusion of mass effects in the initial state only modifies the i = II, IF contributions to ∆ (A) with respect to their original definitions, which is what we focus on in the rest of this section. Finally, the MC@NLO matched fully differential cross section can be written in terms of the previous ingredients as

Sudakov form factor
We now describe explicitly how the ∆ i contributions are constructed in our implementation. As already noted, only the i = II, IF are changed with respect to their original implementation, so in the following we restrict our discussion to them. Most of the ingredients relevant to the matching can be obtained as the four-dimensional limit -ε → 0 -of the equations presented in Sec.2. One comment is in order here. The initial state evolution is partially driven by ratios of PDF factors at different scales, and possibly for different flavours for transitions of quarks to gluons. This may lead to a situation where such a factor reads f g /f Q , i.e. the ratio of a gluon and a heavy quark PDF. Ignoring effects of intrinsic charm and beauty, the quark PDF has no support below its mass threshold, and this ratio becomes ill defined. Different solutions have been constructed in various parton shower algorithms, most of which effectively enforce a splitting such that the heavy quark is replaced by a gluon at threshold. * In their implementation in SHERPA, in particular, they are equal up to phase-space, so that

Initial-Initial configurations
First, we need to express the transverse momentum of the emission -the ordering variable in SHERPA's parton showering -in terms of the variables used to construct the subtraction, x and y: Further, we need the relevant Jacobian factor for the one particle phase-space integration in Eq.(2.18), The fully massive Sudakov form factor for initial-initial dipole configurations is thus given by and V ak,b can be taken from Eq.(2.13).

Initial-Final configurations
Similarly to the previous case we get: This implies that the Jacobian factor becomes 12) which in turn gives the fully massive initial-final contribution to the Sudakov form factor

Physical kinematics
In the practical implementation of the parton shower procedure, the extra emission is attached to an underlying N particle phase space. This corresponds to the reversed procedure used to construct the reduced matrix elements in Sec.2. The new momenta are then obtained from the old ones by inverting Eqs.(2.37),

Initial-Initial configurations
For initial-initial configurations we get: (3.14)

Initial-Final configurations
Similarly, for initial-final configurations, All other configurations are left unchanged by our scheme.

Results
To study the impact of the inclusion of finite-mass effects, we compare results of the five-flavour massive scheme (5FMS) with the vanilla five flavour scheme, where b-quarks are massless. In particular, just as an example and with no intentions of making any statement about any beyond the Standard Model model, we focus on the production of a scalar particle, A, coupling to b quarks through a Yukawa coupling. As a proxy to test the impact of the inclusion of mass effects we vary the mass of this scalar particle in the range 20 − 500 GeV. Further, we let the coupling of the b-quarks to this particle vary too by varying the parameter tan β, mimicking a two Higgs doublet model. We study this process both at fixed-order next-to-leading order accuracy, and at MC@NLO accuracy. Plots for the former are collected in Fig. (1), while the latter set-up is shown in Fig. (2). All our results are obtained within the SHERPA event generator [31]. Leading-order matrix elements, including those of real radiation processes, are calculated using the AMEGIC++ [35] matrix element generator. The differential subtraction follows closely the algorithms of [36], extended with the ingredients reported in the previous Section, Sec. 1. The integrated subtraction terms are implemented in SHERPA and will be made publicly available in a future SHERPA release. Virtual corrections have been obtained from the OPENLOOPS generator [37]. Finally, for the MC@NLO results, we make use of the modifications described in Sec. 3 to SHERPA standard parton shower, the CSS shower [38]. We apply no cuts at generation level. However, in the following we define as b-jet any jet with p ⊥ ≥ 25 GeV that has at least one b-flavoured parton in it. We further require any particle in the final state to have |η| ≤ 2.5. Jets are clustered using the anti-k T [39] algorithm as implemented in FASTJET [40] while the event analysis is performed using the RIVET package [41,42]. We set the renormalization, factorisation and shower starting scales to be µ R = µ F = µ Q = m A /3. For fixed-order predictions, in Fig. (1), we use the -SHERPA default -NNPDF30 NNLO set with α s (m Z ) = 0.118 [43] for the 5FS and the 5FMS. Results presented in Fig. (2), however, where we further compare the two five-flavour scheme predictions with a LO matrix element matched with the NLL DIRE parton shower [26,27], we employ the CT14 NNLO PDF set [44], to satisfy DIRE's requirement of positive definite PDFs, for all predictions. Note that the DIRE prediction is also with five massless flavour, like the standard 5FS.
We start our discussion from the fixed-order results, Fig.1. In this case, we have two mass scales, p ⊥ and m A . Further, as we let them vary, the mass corrections differ in size in different regions of phase-space. In any case, we always expect the 5FS and the 5FMS to converge to each other in the region p ⊥ ∼ m A , which is indeed evident in Fig.1. Note the different p ⊥ range in Fig.1 for m A ≥ 300 GeV to show this expected behaviour. As m A increases we expect a reduction in the absolute size of mass effects, with differences between the two schemes remaining in the region p ⊥ ∼ m b . Looking at the standard (m H = 125, tan β = 1) point, we expect mass effects to play a marginal role, of order ∼ 1 − 5%, at the level of total cross sections. Furthermore, since they are power suppressed, we expect them to be less important at large p ⊥ , while having the largest impact in the lower bins of the distribution. This is due to the fact that the difference in the mass treatment between the two schemes is only in the hard matrix elements and phase-space. This is confirmed in our plot We now turn to the parton shower-matched results, Fig.2. Leaving for a second the discussion of the DIRE prediction aside, in general, we expect a shower prediction to fall on the fixed order result for values of p ⊥ µ Q , where the discussion of fixed order results apply. In this case, too, we expect less important mass effects in the region where both m A and p ⊥ are large, but not for large m A and small p ⊥ . These general considerations are all confirmed in Fig.2. The additional sample is obtained using DIRE. As the implementation of this shower stops the emission for scales ∼ 2 m b , no transverse momentum is generated for the scalar particle when m A = 20 GeV, with µ Q = m A /3 < 2 m b . As for all other configurations the the 5FMS and DIRE predictions are in good agreement, in the region of the Sudakov peak, p ⊥ µ Q .
In this paper we presented all the ingredients that are necessary to construct an extension of VFNS, like the 5FS, in which heavy quarks are treated as massless partons, to allow for massive partons in the initial state, for processes at hadron colliders, at NLO accuracy. In particular we extended the successful Catani-Seymour scheme for subtraction of infrared divergences to the case of massive initial states. In variance to an earlier treatment by Dittmaier, we do not use a finite photon mass as regulator but consistently work in D = 4 + 2ε dimensions. We also re-parametrize the result for the integrated subtraction term in such a way that the residual integral over the light-cone fraction of the emitted particle decouples from the evaluation of the Born cross section, rendering our result more useful for direct implementation. Further, we used these massive dipoles to extend the S-MC@NLO matching as well as the shower generation. We investigated the effect of finite quark masses at fixed-order accuracy in the process bb → A. Mass effects for this process are generally quite small, both at the inclusive and differential level, of the order of a few percent. This however might not hold true for other processes involving heavy quarks. A five flavour massive scheme will provide insight by producing fully differential results including mass effects in a consistent way at matrix-element level. As an additional example we presented simulations for the production of a scalar particle A in bottom-quark fusion, for various configurations of m A and tan β.