Higgs-boson production in top-quark fragmentation

We compute the fragmentation functions for the production of a Higgs boson at $\mathcal{O}(y_t^2\alpha_s)$. As part of this calculation, the relevant splitting functions are also derived at the same perturbative order. Our results can be used to compute differential cross sections with arbitrary top-quark and Higgs-boson masses from massless calculations. They can also be used to resum logarithms of the form $\ln(p_T/m)$ at large transverse momentum $p_T$ to next-to-leading-logarithmic accuracy by solving the DGLAP equations.


Introduction
The associated production of a Higgs boson and a top-quark pair has been studied extensively in theory and by the experimental collaborations at the Large Hadron Collider LHC. This process is not only of interest due to the important role the Higgs boson and the top quark play in searches for models of physics beyond the Standard Model, but also because it provides a means of measuring the top-Yukawa coupling y t directly. Indeed, while currently y t is inferred from measurements by assuming the validity of the Standard Model, see e.g. Refs. [1][2][3], the recent observation of the ttH final state by both the CMS [4] and ATLAS [5] collaborations suggests that a direct measurement of y t using this channel may soon be feasible.
Much effort has gone into improving the precision of the theoretical predictions for top-Higgs associated production: the next-to-leading order (NLO) QCD corrections have been known for a long time [6][7][8][9][10][11] and have since been extended to incorporate the decay of the top quarks, including off-shell effects [12]. Soft-gluon resummation has been performed, first at next-to-leading-logarithmic accuracy (NLL) [13] and then also at NNLL [14][15][16], and NNLO corrections have been computed in the soft-gluon limit [17]. A NLO calculation of ttH+jet is available [18] and NLO electroweak corrections have been calculated as well [19][20][21][22].
Potentially large logarithms of the form ln(p T /m), where m denotes the Higgs-boson or top-quark mass, appear at higher orders in perturbation theory. While measurements at the LHC do not currently probe the kinematic regime where such logarithms would spoil the convergence of the perturbative expansion, such regimes could become relevant for the LHC in the future and will be important at future colliders with higher partonic centre-of-mass energies. The resummation of mass logarithms can be performed using the perturbative fragmentation function formalism [23]. The LO and NLO top-Higgs fragmentation functions have been presented in Ref. [24] in the limit m 2 H m 2 t and based on a soft-gluon approximation. A complete LO calculation of the top-Higgs fragmentation function based on QCD factorisation and including the full mass dependence has been performed in Ref. [25].
A further important benefit of the fragmentation function formalism is that in order to obtain the fully massive cross section in the high-p T regime, the partonic cross sections need only be known in the massless limit m t = m H = 0, considerably simplifying the calculation. While this is only of academic value at NLO, as the NLO corrections are known for general values of the masses, it may be of significant value when attempting to compute the NNLO corrections to associated ttH production. As such, the calculation within the fragmentation function formalism at NLO may be a stepping stone towards the full NNLO corrections at high p T . Of course, in order to perform such computations one needs to know the relevant fragmentation functions (and time-like splitting functions) at the same order in perturbation theory.
In this paper we present the first complete calculation of the O(y 2 t α s ) corrections to the perturbative fragmentation functions and to the splitting functions relevant for associated top-Higgs production, namely those describing the final-state transitions t → H and g → H. The paper is organised as follows. The fragmentation function formalism is briefly presented in Section 2. In Section 3, we provide the details of the calculation. Section 4 describes the results, and we conclude in Section 5.

The perturbative fragmentation function formalism
Fragmentation functions were initially introduced to describe the production of hadrons at high p T [26]. In analogy to the initial state factorisation theorem, which states that hadron-collision cross sections can be described, up to power corrections, by a convolution of non-perturbative parton distribution functions and perturbative hard scattering matrix elements, the production of hadrons can be factorised into non-perturbative fragmentation functions D i→h , which describe the transition from the parton i to the hadron h, and perturbative hard scattering matrix elements: The presence of collinear singularities implies that fragmentation functions need to be renormalised. This collinear renormalisation is of the form (see also Section 3.2) and leads to a set of renormalisation group equations for fragmentation functions, called the (time-like) Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) evolution equations [27][28][29]: where P T ij are the time-like splitting functions and µ F r is the final-state factorisation (or fragmentation) scale.
In Ref. [23], it was pointed out that the same formalism can be used to describe the production of a massive quark starting from a massless calculation. The difference is that the production of a heavy quark can be described perturbatively, and the corresponding fragmentation functions can thus be calculated in perturbation theory. The main benefit of factorising the production of the massive quark into fragmentation functions is that this allows one to use the DGLAP equations to resum mass logarithms: whereas the fully massive calculation contains logarithms of the type ln Q 2 /m 2 , the perturbative fragmentation function (PFF) formalism splits the calculation into fragmentation functions, which contain the logarithm ln µ 2 F r /m 2 , and massless cross sections, which contain the logarithm ln Q 2 /µ 2 F r . It is thus possible to evaluate the fragmentation functions at a scale µ F r ∼ m and then use the DGLAP equations to obtain them at a scale µ F r ∼ Q without introducing any large logarithms. The accuracy of the resummation of mass logarithms is given by the order of the splitting functions used, i.e. leading-logarithmic accuracy (LL) for LO splitting functions, NLL for NLO splitting functions, etc.
While the original derivation of the PFFs at NLO in Ref. [23] involved comparing the massive and massless cross sections, the derivation of the NNLO heavy-quark FFs in Refs. [30,31] used a definition equivalent to the field-theoretic definition given in Ref. [32]. This formalism can be used for the production of any heavy particle and was first applied to Higgs-boson production in Ref. [25]. There, the LO FF and splitting function for the transition t → H, and of several other transitions, was derived by directly employing the definition of FFs given in Ref. [32]. The fully-differential factorisation formula for this case reads where the sum is over all QCD partons (including the top quark) and the Higgs boson, assuming only QCD corrections are considered.

Details of the calculation
In this section we present our calculation in considerable detail. We set the notation, discuss the renormalisation procedure and, in particular, describe the techniques used to perform the calculation of the Feynman integrals.

Notation and definitions
If n µ is an arbitrary four-vector satisfying n 2 = 0 and n · n = 1, where n µ = n µ , then any four-vector p µ can be decomposed as follows: The bare fragmentation function D B q→H (z) for a quark q to fragment to a Higgs boson H with momentum p µ H with p + H = zp + q and p µ H,T = 0 is defined as follows [32]: where P (P) is the (reverse) path-ordering operator. Here, 'bare' refers to the need to perform collinear renormalisation, as explained in Section 3.2. The bare gluon fragmentation function D B g→H (z) is defined as [32] Using these definitions, fragmentation functions can be calculated using standard higher-order techniques for the computation of cross sections. We use dimensional regularisation with d = 4 − 2 dimensions and the MS definition of the renormalisation scale: At NLO, the top-quark mass (and the Yukawa coupling) as well as the wave-function are renormalised on-shell. We work in the linear gauge with an arbitrary gauge parameter.
The computation of D B t→H at O(y 2 t ) using eq. (3.2) is simple enough to reproduce here explicitly. There is only a single diagram, shown in fig. A1, which yields (cf. Ref. [25]) (3.5) Note that the terms O( ) of this LO result are required for the different renormalisation contributions to the NLO result.

Collinear renormalisation
The parton-to-Higgs FFs are renormalised in full analogy to the more common heavy-quark case, presented through NNLO in e.g. Ref. [33]. Nonetheless, the collinear renormalisation matrix is reproduced here for the case of Higgs fragmentation. We perform collinear renormalisation in the MS scheme. Setting µ F r = µ R for convenience, the NLO QCD MS collinear renormalisation of D t→H is given by Here has been used, P (0)T qq is the standard LO time-like quark-to-quark splitting function, and which implies β (0) where γ 0 is the leading-order quark-mass anomalous dimension. The NLO QCD MS collinear renormalisation of the gluon-to-Higgs FF is given by The factor 2 in the second line of eq. (3.13) accounts for both top and anti-top quarks, and is the usual LO time-like gluon-to-quark splitting function. P (0)T tH (z) was first calculated in Ref. [25]. It can also be obtained from the results above: since the poles in of the bare FF must be removed by the collinear renormalisation and the only counterterm at LO is proportional to P This result, together with the well-known QCD splitting functions and the top-Yukawa β-function, allows us to predict all the poles of the bare Higgs fragmentation functions at NLO, except for the terms proportional to P (1)T tH and P (1)T gH , since both of these functions were previously unknown. However, similarly to what was done for the LO case above, both of these splitting functions can be readily derived from the corresponding NLO FFs and have thus been obtained as a by-product of the computation presented here.

Outline of the calculation and computational techniques
The general outline of the calculation is identical for both FFs. We start by determining all relevant Feynman diagrams, as shown in Appendix A. All subsequent symbolic manipulations were performed using Mathematica [34]. The Mathematica package FeynCalc [35] was used to perform the contraction of Lorentz indices and Dirac traces. On-shell conditions and the momentum fraction δ-distribution coming from the definition of the fragmentation function were replaced with cut propagators using reverse unitarity [36][37][38][39][40]: For each part of the calculation (D g→H , the real corrections to D t→H and the virtual corrections to D t→H ), all required Feynman integrals can be written as the product of 7 propagators raised to integer powers. These integrals are then reduced to a smaller set of master integrals using IBP relations [41,42]. This reduction was performed using FIRE6 [43]. The master integrals were calculated explicitly in the m t → ∞ limit (or equivalently m H → 0) up to the required order in , making use of HypExp [44][45][46][47] to expand hypergeometric functions in where necessary. Subsequently, the method of differential equations [48][49][50][51][52][53][54] was used to obtain the master integrals for finite values of m t . The set of master integrals was constructed to be in canonical form, significantly simplifying this step [55]. The details of the canonical basis are described in Section 3.4. The list of master integrals for the top-to-Higgs FF, as well as their -expansions through the order required for the present calculation, can be found in Appendix D. The solution of the differential equations was performed using tools provided by PolyLogTools [56]. The results, initially involving multiple polylogarithms (MPLs) of weight three, were simplified using the strategy of symbols [57][58][59] outlined in detail in Ref. [58] and summarised for the present case in Section 3.5, again making use of tools provided by PolyLogTools. The final results are expressed in terms of a comparatively small number of classical polylogarithms.
The calculation was performed under the assumption that m H < √ 2m t . As such, the results presented below are only valid for values of the masses satisfying this condition, which however includes the values relevant for the Standard Model. Of course, the NLO splitting functions are valid for any values of the masses.

The canonical basis
As mentioned in the previous section, the differential equations method was used to derive the top mass dependence of the master integrals. The system of first-order linear differential equations with respect to m 2 t reads where g(m 2 t , ) is the vector of master integrals and A(m 2 t , ) is a N × N matrix, with N being the dimension of g(m 2 t , ). Note that for simplicity of notation the dependence on kinematic invariants other than m 2 t is not explicit in eq. (3.20). This convention will be used henceforth.
A modern approach for the solution of eq. (3.20) consists in the application of the canonical basis method. For the set of master integrals arising from the reduction of both the real and the virtual contributions to D (1) t→H , a change of basis was found such that where f (m 2 t , ) is the vector of canonical master integrals, chosen to be free of -poles, dÃ(m 2 t ) =Ã(m 2 t ) dm 2 t , and the -dependence is factorised so thatÃ(m 2 t ) depends only on the kinematic invariants. Moreover, we define the canonical basis elements f i (m 2 t , ) as the product of the so-called pre-canonical integrals I i (m 2 t , ), the algebraic pre-factors , which are functions of the kinematic invariants, and the dimensional regularisation parameter to the n i -th power: The canonical basis was obtained by the use of the following semi-algorithmic approach: • The pre-canonical integrals I i (m 2 t , ) were chosen by modifying the powers of the denominators of the master integrals with the aim of maximising their symmetries.
The Feynman diagram representation of the pre-canonical integrals is given in Appendix B, while their analytic expressions are provided in Appendix D.
Starting from eq. (3.21), the change of variable m 2 t → m 2 H (1 − τ 2 )/4 was performed in order to linearise the square-root dependence on m 2 t . Hence, the canonical master integrals of both the real and the virtual corrections to the top-to-Higgs FF satisfy the following system of differential equations: (3.25) The differential equation matrix can be written in dlog form: where the matrices M i are 8 × 8 matrices with purely rational entries. The arguments α i and the matrices M i are given in Appendix C.
The integration of the master integrals was performed order by order in in terms of MPLs. The integration constants are m H and z dependent; therefore the boundary conditions were computed by calculating the master integrals in the m H → 0 limit.

Simplifying polylogarithmic expressions using symbols
As explained in Refs. [57][58][59], the symbol map is an operation which maps polylogarithmic functions, e.g. the classical polylogarithms and MPLs, onto a direct product of weight-1 objects, i.e. logarithms. Note that the constants π n and ζ(n) count as weight-n objects as well. For example, the symbol of a classical polylogarithm is given by 1 (3.27) A few properties of symbols are important for the present discussion. First, the symbol map is linear, i.e. symbol(a + b) = symbol(a) + symbol(b) . (3.28) Second, the operator ⊗ is distributive, i.e.
giving the important relation Third, the symbol of a product of two functions is given by the shuffle product of the symbols: where Σ(n, m) is the set of all permutations σ of 1, ..., n + m satisfying σ(1) < ... < σ(n) , σ(n + 1) < ... < σ(n + m) . (3.32) Another property relevant here is that ζ(n) and higher powers of π are mapped to zero. The final property needed is that if the symbols of two functions are identical, then the functions are identical up to contributions in the kernel of the symbol map, i.e. terms containing factors of ζ(n) or π n for n > 1. Using these properties, and with the help of tools implemented in PolyLogTools, one can now simplify the initial form of the master integrals as resulting from the solution of the differential equations via the following steps. First, one takes the symbol of the initial form. This allows one to work with logarithms and their simple relations, i.e. eq. (3.30), reducing the number of distinct logarithms to a smaller set with simpler arguments. One then searches for a linear combination of classical polylogarithms with simple arguments with the same symbol. This linear combination must be identical to the initial form of the master up to terms in the kernel of the symbol map, of which there are only a few at low weights and whose numerical coefficients are rational numbers. These terms can thus be obtained by fitting their coefficients over the rational numbers. Finally, a high-precision numerical comparison of the initial form and the simplified one over the entire relevant region of the parameter space 2 is used to confirm the fitted coefficients. Note that it is well known that through weight three, all MPLs can be written as linear combinations of classical polylogarithms, while this is not the case for higher weights. Since the present problem does not contain higher-weight MPLs, classical polylogarithms will always be sufficient.
Finding a function with the same symbol is the most laborious part of the procedure above. Fortunately, as explained in Ref. [58], different types of contributions satisfy different symmetry relations with respect to the various entries of the symbol tensor, allowing one to isolate a specific type of contribution, find its simplified form and subtract it from the full master, after which one can move on to the next type of contribution. For example, at weight two, there are four types of contributions: (3.33) Initially ignoring contributions containing factors of π leaves Li 2 (a) and ln(a)ln(b). The symbol of the latter is symmetric with respect to the first and the second entry of the symbol: whereas the symbol of Li 2 (a) is not (cf. eq. (3.27)). Anti-symmetrising with respect to these entries thus only leaves a part given by dilogarithmic contributions, which can now be determined more easily. After subtracting the dilogarithmic part, all contributions to the symbol are of the form ln(a)ln(b) or π ln(a) and can essentially be directly read off from the symbol. The reason why the contributions involving π are initially ignored is that π behaves in an exceptional way within the symbol, making (anti-)symmetrisation impossible. The only contributions left must be proportional to π 2 and can thus be trivially determined via a numerical comparison. At weight three, there are 8 contributions: Anti-symmetrising with respect to the last two entries of the symbol leaves only ln(a)Li 2 (b) and π Li 2 (a), with the first entry determined by the logarithms and the last two entries by the dilogarithms. After subtracting these contributions from the full result, antisymmetrisation with respect to the first two entries of the symbol (ignoring again terms containing π) leaves only the trilogarithms. After those have been determined and subtracted, only the purely logarithmic terms remain and can be read off directly. The terms of the type π 2 ln(a), π 3 and ζ(3) must again be fitted. Here, the knowledge of which logarithms showed up in the simplification so far is used as a guide for the Ansatz: indeed, no logarithms were needed for these fits that were not already needed for the contributions that can be determined using the symbol. In theory, terms proportional to π 3 or ζ(3) can only be separated by requiring their rational coefficients to be reasonable, after which the high-precision comparison of the initial and final forms is used to confirm their correctness.
In practice, at least for the present computation, the π 3 -terms are always accompanied by the imaginary unit i, while the ζ(3)-terms are not, allowing for a much simpler separation.

Results
The non-zero NLO parton-to-Higgs splitting functions are As required, they are independent of m t and m H , which is an important test of their validity. They have been extracted from the full calculation, as well as from simplified calculations where either m t or m H is set to zero. The results agree in all three cases, providing a further cross-check of the results.
Defining the shorthand notation the NLO corrections to the fragmentation functions read .  H > 0). However, the coefficients of the divergences are numerically tiny, so that they only become relevant for values of z extremely close to 0 or 1, as is apparent from fig. 1. The fragmentation function D g→H is also logarithmically divergent for z → 0, while for z → 1 it tends to 0. As the magnitude of the gluon-to-Higgs FF is consistently much smaller than the top-to-Higgs FF, gluon-induced contributions can safely be neglected as long as their suppression is not compensated for by an enhanced partonic cross section, in analogy to the non-singlet approximation commonly used for heavy-quark fragmentation. While for hadron collisions the gluon-induced cross section is indeed enhanced with respect to the top-quark production cross section, an approximation where gluon-induced contributions are neglected is appropriate for lepton collisions.
One might be tempted to compute higher-order corrections to the Higgs FFs in the m H = 0 limit as a first approximation. While m H m t is not satisfied in nature, this has nonetheless proven to be a surprisingly good approximation in other calculations, see e.g. Ref. [61] for a review. However, the LO result already demonstrates that this is not the case here. Indeed, the form of the LO FF suggests that m H /(zm t ) would be a more appropriate expansion variable. The same behaviour holds at NLO, suggesting that the m H = 0 limit is only a good approximation for large z 0.7. This is born out by the numerical results for the NLO fragmentation functions corresponding to m H = 0 and m H > 0 shown in fig.  1.

Conclusions
In this paper we have presented the fragmentation and splitting functions for the final-state transition from a top-quark to a Higgs boson or from a gluon to a Higgs boson through NLO in QCD, i.e. through O(y 2 t α s ). These new results can be used in two different ways to obtain higher-precision theoretical predictions for any process involving these transitions. First, any such process can be described at NLO in QCD, up to power corrections in the top-quark or Higgs-boson mass, by performing the convolution of massless coefficient functions with the fragmentation functions. As setting all masses to zero can significantly simplify the computation of coefficient functions, this may allow to obtain higher-order results which would be prohibitively complicated for the massive case. Second, the results can be used to perform a resummation of logarithms of the type ln(p T /m) at NLL accuracy using the DGLAP equations, which is necessary to maintain the perturbative convergence of the predictions for top-Higgs associated production at large transverse momentum.
To obtain our results, we started from the field-theoretic definitions of the fragmentation functions given in Ref. [32]. All employed computational techniques (e.g. reverse unitarity, IBP reduction and the solution of master integrals using differential equations) have become increasingly standard for multi-loop calculations in recent years, and most of them were already applied in Refs. [30,31] to the computation of the NNLO corrections to the heavy-quark fragmentation functions. It has been demonstrated that in order to obtain accurate results, one cannot neglect the Higgs-boson mass in the computation of the Higgs fragmentation functions. At the same time, the results suggest that neglecting fragmentation contributions from gluons may be justified at lepton-colliders. In subsequent work we plan to study the size of missing power corrections, the effect of NLL DGLAP resummation and the size of gluoninitiated contributions for top-Higgs associated production at the LHC, to explore the phenomenological importance of those effects. I 0,0,0,1,1,2,1 I 0,0,0,1,2,1,1 I 0,0,1,1,1,2,1 I 0,0,1,2,1,1,1 I 0,0,1,1,1,1,1 I 0,1,2,0,1,1,1 I 0,1,1,0,1,1,2 I 0,1,1,0,1,1,1 Figure A5: The pre-canonical master integrals for the real corrections to D t→H . The black dots represent squared propagators. Figure A6: The scalar integral topology for the virtual corrections to D t→H defined in eq. (D.12).
C Canonical form

C.1 Canonical form for real corrections
The differential equation matrix in dlog form for the real corrections is:

C.2 Canonical form for virtual corrections
The differential equation matrix in dlog form for the virtual corrections is: (C.14)

D.1 Masters for the real corrections
The master integrals for the real corrections are defined as follows: where the prefactor −16πiµ 4 has been chosen to match the factor resulting from the use of reverse unitary and the couplings in the full result. Note that the numerator is merely a constant introduced to eliminate overall factors of n · p H . A propagator with a subscript c denotes a cut propagator: (D.2) In this notation, the masters are: Their expansions in through the order required for the present calculation are given below. They are also included as ancillary files with the e-print version of this paper.