# On the interference of ggH and c\bar{c}H Higgs production mechanisms and the determination of charm Yukawa coupling at the LHC

FFebruary 9, 2021

Prepared for submission to JHEP

TTP21-003, P3H-21-010

On the interference of ggH and c ¯ cH Higgs productionmechanisms and the determination of charm Yukawacoupling at the LHC

Wojciech Bizo´n, a,b

Kirill Melnikov, a J´er´emie Quarroz a a Institute for Theoretical Particle Physics (TTP), Karslruhe Institute of Technology, D-76128 Karl-srue, Germany b Institute for Astroparticle Physics, Karslruhe Institute of Technology, D-76344 Eggenstein-Leopoldshafen,Germany

Abstract:

Higgs boson production in association with a charm-quark jet proceeds throughtwo diﬀerent mechanisms – one that involves the charm Yukawa coupling and the other thatinvolves direct Higgs coupling to gluons. The interference of the two contributions requires ahelicity ﬂip and, therefore, cannot be computed with massless charm quarks. In this paper, weconsider QCD corrections to the interference contribution starting from charm-gluon collisionswith massive charm quarks and taking the massless limit, m c →

0. The behavior of QCDcross sections in that limit diﬀers from expectations based on the canonical QCD factorization.This implies that QCD corrections to the interference term necessarily involve logarithms ofthe ratio M H /m c whose resummation is currently unknown. Although the explicit next-to-leading order QCD computation does conﬁrm the presence of up to two powers of ln( M H /m c )in the interference contribution, their overall impact on the magnitude of QCD correctionsto the interference turns out to be moderate due to a cancellation between double and singlelogarithmic terms. a r X i v : . [ h e p - ph ] F e b ontents O (ln m c ) -enhanced contributions in the real corrections 18 A.1 Integration of the soft-quark subtraction terms 20A.2 Integration of the quasi-collinear subtraction terms 20A.3 Numerical checks 23

B Soft-quark integrals 24 – 1 –

Introduction

Studies of Yukawa couplings play an important role in the veriﬁcation of the mechanism ofelectroweak symmetry breaking as described by the Standard Model. By now, Higgs couplingsto bottom and top quarks, as well as to tau leptons and muons, have been measured to aprecision of about twenty percent [1–6]. Within the error bars, the measured values for allfour Yukawa couplings are consistent with the Standard Model predictions.However, the Yukawa couplings to lighter fermions have not been studied experimentally.Although it is generally agreed that the Yukawa couplings of electrons and up, down andstrange quarks can be observed if and only if they enormously deviate from their StandardModel values, the situation with the charm Yukawa is not so hopeless. In fact, it appearsthat with the full LHC luminosity, the charm Yukawa coupling can be measured if its valuedeviates from the Standard Model expectation by an order one factor [7]. Diﬀerent observablesto measure the charm Yukawa coupling at the LHC have been proposed; they include inclusive( H → c ¯ c ) and exclusive ( H → J/ψ + γ and similar) decays of the Higgs boson [8–10], themodiﬁcations of the Higgs transverse momentum distribution [11] in the gg → H + X processand, ﬁnally, Higgs boson production cross section in association with a charm jet [12].In this paper we focus on the latter process, pp → H + jet c . At leading order in perturba-tive QCD, Higgs bosons are produced in association with charm jets in the partonic process cg → Hc . The amplitude of this process receives contributions proportional to the charmYukawa coupling and to an eﬀective ggH coupling M ∼ g Yuk M + g ggH M , (1.1)see Figure 1. As a result, the pp → H + jet c cross section contains the interference term σ Hc ∼ g ˜ σ + g ggH ˜ σ + g Yuk g ggH ˜ σ Int . (1.2)It can be expected that a reliable description of Higgs boson production in association witha charm jet can be obtained by systematically computing the diﬀerent terms in Eq. (1.2)to higher orders in perturbative QCD. In fact, it is emphasized in Ref. [12] that the largesttheoretical uncertainty in using H + jet c production cross section to constrain charm Yukawacoupling is related to perturbative QCD uncertainties so that it seems natural to computehigher order QCD corrections to σ Hc in Eq. (1.2).However, pursuing this program for the interference term in Eq. (1.2) is quite subtle aswe now discuss. Indeed, perturbative computations in QCD are performed with masslessincoming partons. In case of the massless charm quark that, however, has non-vanishingYukawa coupling to the Higgs boson, the interference term in Eq. (1.2) vanishes and weobtain lim m c → σ Hc ∼ g ˜ σ + g ggH ˜ σ . (1.3)This happens because the Yukawa interaction ﬂips charm’s helicity but the gluon-charminteraction conserves it; hence, the two contributions in Eq. (1.1) cannot interfere if m c = 0.– 2 – H (a) Yukawa coupling H (b) Eﬀective ggH vertex

Figure 1:

Leading-order Feynman diagrams contributing to the pp → Hc process. Wedistinguish two separate production mechanisms: one that is driven by the Yukawa coupling(left) and the other one that requires direct coupling of Higgs to gluons (right).For the massive charm quark the interference does not vanish and is proportional to the charmmass in the ﬁrst power. Whether or not the interference contribution is negligible depends onthe relative magnitude of the two amplitudes in Eq. (1.1). Leading-order computations withmassless quarks show that the charm-Yukawa independent amplitude g ggH M in Eq. (1.1) islarger than the charm-Yukawa dependent one g Yuk M suggesting that the interference maybe non-negligible.It is straightforward to calculate the interference at leading order in perturbative QCD.Indeed, the interference requires one helicity ﬂip on a charm line that connects initial andﬁnal states; this ﬂip is accomplished by a single mass insertion. This implies that one cancompute the interference of the two amplitudes using massive charm quarks, take the m c → m c . Since we require acharm jet in the ﬁnal state, none of the kinematic invariants of the cg → Hc process can besmall. Hence, once one power of m c is extracted, the rest of the leading-order calculation ofthe interference contribution can be performed using the standard approximation of massless(charm) quarks. Such calculation, that we describe in Section 5, shows that the leading-orderinterference amounts to about ten percent of the contribution to the H + jet c cross sectionthat is proportional to the Yukawa coupling squared.Although the interference contribution is not large, it is worth thinking about it at next-to-leading order (NLO) in perturbative QCD since there are reasons to believe that theinterference contribution is perturbatively unstable , at variance with the two other contribu-tions to pp → H + jet c cross section. Indeed, even if we require an energetic charm jet in theﬁnal state, soft and collinear kinematic conﬁgurations lead to logarithmic sensitivity of theinterference to the charm mass m c . Hence, before the m c → m c have to be extracted fromboth real and virtual corrections to the interference part of the production cross section.One may argue that, since the ﬁnite charm mass provides yet another way to regulatecollinear divergences, it is to be expected that the procedure described above will lead toa familiar picture of (quasi)-collinear factorization of QCD amplitudes. If so, all ln( m c )–dependent terms should disappear once infrared safe cross sections and distributions arecomputed using short-distance quantities, including conventional parton distribution func-tions (PDFs). However, we will show that for the interference contribution this expectation– 3 –s invalid and that well-known formulas that describe collinear factorization of mass singu-larities are not applicable in that case. We will also show that the helicity ﬂip leads to anappearance of soft-quark singularities that, interestingly, make jet algorithms logarithmically-sensitive to m c .There are two consequences of the above discussion. First, the problem of estimating themagnitude of the interference contribution to the production of a Higgs boson in associationwith a charm jet turns into an interesting problem in perturbative QCD that borders onsuch important issues as soft and collinear QCD factorization for mass power corrections[13–17]. Second, a more complex pattern of this factorization, as compared to the canonicalcollinear and soft cases [18], implies that NLO QCD corrections to leading-order interferenceare enhanced by up to two powers of a large logarithm ln Q/m c where Q is a typical hardscale in the process pp → H + jet c . For this reason NLO QCD corrections to the interferencemay be expected to be signiﬁcant and it becomes essential to explicitly compute them. Thisis what we set out to do in this paper.The rest of the paper is organized as follows. In the next section we derive a relationbetween MS-regulated and mass-regulated parton distribution functions at O ( α s ) using theprocess of Higgs boson production in c ¯ c annihilation. We use the established relation toremove “conventional” collinear logarithms from NLO QCD corrections to the interferencecontribution to the production of Higgs boson in association with charm jet. In Section 3we discuss factorization of mass singularities in the interference contribution to cg → Hc process and show that it works diﬀerently as compared to the standard case [18]. In Section 4we brieﬂy describe the technical details of the calculation of NLO QCD corrections to theinterference contribution. In Section 5 we present phenomenological results and discuss therelative importance of logarithmically-enhanced terms. We conclude in Section 6. Additionaldiscussion of soft and collinear limits of the interference contributions as well as some relevantsoft integrals can be found in several appendices. It is well-known that quark masses screen collinear singularities. For this reason we canthink about small quark masses as a particular choice of a collinear regulator. Since, whendescribing “leading-twist” inclusive partonic processes, collinear sensitivity either cancels outor is absorbed into parton distribution functions, it is possible to derive relations betweenparton distribution functions that are used for computations with nearly massive and strictlymassless quarks by requiring that predictions for physical processes are independent of acollinear regulator. To derive a relation between “massive” and “massless” PDFs, we start with the produc-tion of a Higgs boson in an annihilation of two massive charm quarks and write the diﬀerential We note that a derivation of the initial condition for the electron structure function in QED was recentlypresented in Ref. [19]. There is a strong conceptual overlap of the discussion in that reference and thecomputation reported in this section. – 4 –ross section as d σ pp → H = (cid:88) ij (cid:90) d x d x f ( m ) i ( x ) f ( m ) j ( x )dˆ σ ( m ) ij → H + X . (2.1)Here f ( m ) i are parton distribution functions and the superscript m implies that all relevantquantities should be computed using quark masses as collinear regulators. Also, ˆ σ ( m ) ij → H + X is the partonic diﬀerential cross-section. At leading order i ( j ) = c, ¯ c ; at higher orders otherchannels also contribute.Calculation at leading order in α s is straightforward since the leading-order cross section σ ( m ) c ¯ c → H has a regular m c → α s there is no diﬀerencebetween f ( m ) i and conventional MS parton distribution functions, i.e. f ( m ) i = f MS i .The situation becomes more complicated at next-to-leading order where the charm quarkmass screens collinear singularities; hence, our goal is to re-write the NLO QCD contributionsto the cross section c ¯ c → H + X in such a way that logarithms of m c are extracted explicitly.We begin by considering the process c ( p ) + ¯ c ( p ) → H + g ( p ) and treating charmquarks as massive. Kinematic regions that lead to soft and (quasi-)collinear singularities arewell understood. The behavior of matrix elements in these limits is described by conventionalfactorization formulas [18]. We can deﬁne a hard m c -independent cross section by subtractingthe singular limits. When the subtracted terms are added back and integrated over unresolvedparts of the Hg phase space, logarithms of the charm mass appear. This procedure is iden-tical to methods developed to extract infrared and collinear singularities from real emissioncontributions to partonic cross sections. Its application in the present context allows us toexplicitly extract logarithms of the charm mass.To organize the calculation, we follow the nested soft-collinear subtraction scheme [20–23]which, at next-to-leading order, is equivalent to the FKS scheme [24, 25]. We use dimensionalregularization to regularize soft singularities and the charm mass to regularize the collinearones. Using notations from Ref. [20], we write the partonic cross section for c ( p ) + ¯ c ( p ) → H + g ( p ) as2 s · dˆ σ c ¯ c → H + g = (cid:90) [d g ] F LM (1 c , ¯ c ; 3 g ) ≡ (cid:104) F LM (1 c , ¯ c ; 3 g ) (cid:105) = (cid:104) S F LM (1 c , ¯ c ; 3 g ) (cid:105) + (cid:104) ( C + C )( I − S ) F LM (1 c , ¯ c ; 3 g ) (cid:105) + (cid:104) ( I − C − C )( I − S ) F LM (1 c , ¯ c ; 3 g ) (cid:105) . (2.2)The key observation is that since the fully-regulated (last) term in Eq. (2.2) is free ofboth soft and quasi-collinear singularities, the limit m c → p · p ∼ p · p →

0. It reads [18] S F LM (1 c , ¯ c ; 3 g ) ≈ g s C F (cid:18) p · p )( p · p )( p · p ) − m c ( p · p ) − m c ( p · p ) (cid:19) F LM (1 c , ¯ c ) , (2.3) The space-time dimension d is parametrized as d = 4 − (cid:15) . – 5 –here g s is the unrenormalized strong coupling constant.Since the soft gluon decouples from the function F LM (1 c , ¯ c ) we can integrate Eq. (2.3)over the gluon phase space. We work in the center-of-mass frame of the colliding charmpartons and parametrize their energies as E = E = E . The center of mass energy squaredin the massless approximation is then s = 4 E . We also cut integrals over gluon energy at E = E max , cf. Ref [20]. Integrating Eq. (2.3) over gluon phase space [d g ] and taking the m c → (cid:104) S F LM (1 c , c ; 3 g ) (cid:105) = − C F [ α s ] E − (cid:15) max (cid:15) (cid:2) I m ( E ) − I m ( E ) (cid:3) (cid:104) F LM (˜1 c , ˜2 ¯ c ) (cid:105) , (2.4)where [ α s ] = g s Ω ( d − / (2(2 π ) d − ) and Ω ( d − is the solid angle of the ( d − The notation ˜1 c (˜2 ¯ c ) implies that the corresponding four-momenta should be taken inthe massless approximation. The two integrals I m (2 m ) in Eq. (2.4) read I m ( E ) = (cid:90) − d cos θ (sin θ ) − (cid:15) − β cos θ ≈ − − (cid:15) Γ (1 − (cid:15) ) (cid:15) Γ(1 − (cid:15) ) (cid:34) − Γ(1 + (cid:15) )Γ(1 − (cid:15) )Γ(1 − (cid:15) ) (cid:18) m c E (cid:19) − (cid:15) (cid:35) ,I m ( E ) = m c E (cid:90) − d cos θ (sin θ ) − (cid:15) (1 − β cos θ ) ≈ (cid:18) m c E (cid:19) − (cid:15) Γ(1 − (cid:15) )Γ(1 + (cid:15) ) , (2.5)where β = (cid:112) − m c /E and we neglected all power-suppressed terms when writing theresults.Collinear subtraction terms contain quasi-collinear singularities. The two collinear limitscorrespond to two distinct cases, p · p ∼ m c → p · p ∼ m c →

0. They read C i F LM (1 c , ¯ c ; 3 g ) = g s ( p i · p ) (cid:20) P qq ( z ) − C F m c z ( p i · p ) (cid:21) F ( i )LM (˜1 c , ˜2 ¯ c ; z ) z , (2.6)where F ( i )LM (˜1 c , ˜2 ¯ c ; z ) = δ i F LM ( z · ˜1 c , ˜2 ¯ c ) + δ i F LM (˜1 c , z · ˜2 ¯ c ) , (2.7)and P qq ( z ) = C F (cid:18) z − z − (cid:15) (1 − z ) (cid:19) (2.8)is the collinear splitting function. The variable z is deﬁned as z = ( E i − E ) /E i with i = 1 , I m (2 m ) shown The solid angle of the d -dimensional space is Ω ( d ) = 2 π d/ / Γ (cid:0) d/ (cid:1) . – 6 –n Eq. (2.5). Performing the soft subtraction of the collinear-subtracted cross section, we ﬁnd (cid:104) C ( I − S ) F LM (1 c , ¯ c ; 3 g ) (cid:105) = [ α s ] E − (cid:15) (cid:90) d z I ( E, z ) (cid:28) F LM ( z · ˜1 c , ˜2 ¯ c ) z (cid:29) , (2.9)where I ( E, z ) = I m ( E ) ¯ P qq ( z ) − I m ( E ) ¯ P ( m ) qq ( z ) , (2.10)and ¯ P qq = C F (cid:18) z (1 − z ) (cid:15) − (cid:15) (1 − z ) − (cid:15) + 1 (cid:15) δ (1 − z ) e − (cid:15)L (cid:19) , ¯ P ( m ) qq = C F (cid:18) z (1 − z ) (cid:15) + 12 (cid:15) δ (1 − z ) e − (cid:15)L (cid:19) , (2.11)and L = ln( E max /E ). A similar expression can be written for (cid:104) C ( I − S ) F LM (1 ,

2; 3) (cid:105) .It is straightforward to combine soft and collinear contributions and to expand them in (cid:15) . At this point, it is convenient to switch to a strong coupling constant renormalized at thescale µ . We obtain (cid:104) S F LM (1 ,

2; 3) (cid:105) + (cid:104) C ( I − S ) F LM (1 ,

2; 3) (cid:105) + (cid:104) C ( I − S ) F LM (1 ,

2; 3) (cid:105) == C F α s ( µ )2 π (cid:26) − L c (cid:15) − L c + 1 − π − L c + 2 L c ln M H µ (cid:27) (cid:104) F LM (1 , (cid:105) + C F α s ( µ )2 π (cid:90) d z (cid:26)(cid:20) z − z (cid:21) + L c + (1 − z ) (cid:27) (cid:28) F LM ( z · ˜1 c , ˜2 ¯ c ) z + F LM (˜1 c , z · ˜2 ¯ c ) z (cid:29) , (2.12)where L c = ln( M H /m c ) − c ¯ c → H cross section, we need to includevirtual corrections. They are computed in a standard way (see e.g. Ref. [26] where such acomputation is reported); the result is then expanded around m c = 0. We renormalize theYukawa coupling in the MS scheme at the scale µ . The result reads2 s · dˆ σ V = C F α s ( µ )2 π (cid:20) (cid:15) L c + 4 π L c + (2 L c + 3) ln µ M H + 3 L c (cid:21) (cid:104) F LM (˜1 c , ˜2 ¯ c ) (cid:105) . (2.13)Upon combining virtual, soft, collinear and fully-regulated terms, we obtain the followingNLO QCD contribution to c ¯ c → H + X cross section2 s · dˆ σ NLO = (cid:104) ( I − C − C )( I − S ) F LM (˜1 c , ˜2 ¯ c ; 3 g ) (cid:105) + C F α s ( µ )2 π (cid:18) π − M H µ (cid:19) (cid:10) F LM (˜1 c , ˜2 ¯ c ) (cid:11) + α s ( µ )2 π (cid:88) i =1 1 (cid:90) d z (cid:26) P APqq ( z ) ln µ m c + P ﬁn qq ( z ) (cid:27) (cid:42) F ( i )LM (˜1 c , ˜2 ¯ c ; z ) z (cid:43) , (2.14)– 7 –here P AP qq = C F (cid:20) z − z (cid:21) + and P ﬁn qq = C F (cid:20) z − z (cid:21) + (cid:18) ln M H µ − (cid:19) + C F (1 − z ) . (2.15)The ﬁrst term on the right hand side of Eq. (2.14) is the hard inelastic contribution; it canbe computed directly in the massless limit, m c = 0. The second term on the right hand sideof Eq. (2.14) is the soft-virtual piece; it describes kinematic conﬁguration that is equivalentto the leading-order one. The third term in Eq. (2.14) describes kinematic conﬁgurationsthat are boosted relative to the leading-order ones; we note that the residual logarithmicdependence on m c is present in these contributions only .A similar computation for massless charm partons requires, in addition, a collinear PDFrenormalization to make the partonic cross section collinear-ﬁnite and independent of theregularization parameter (cid:15) . The result reads2 s · dˆ σ m c =0NLO = (cid:104) ( I − C − C )( I − S ) F LM (˜1 c , ˜2 ¯ c ; 3 g ) (cid:105) + C F α s ( µ )2 π (cid:20) π − M H µ (cid:21) (cid:104) F LM (˜1 c , ˜2 ¯ c ) (cid:105) + α s ( µ )2 π (cid:88) i =1 1 (cid:90) d z P ( (cid:15) ) qq ( z ) (cid:42) F ( i )LM (˜1 c , ˜2 ¯ c ; z ) z (cid:43) , (2.16)where P ( (cid:15) ) qq ( z ) = C F (cid:20) z − z (cid:21) + ln M H µ + 2 C F (cid:20) z − z ln(1 − z ) (cid:21) + + C F (1 − z ) . (2.17)The partonic cross sections in Eqs. (2.14) and (2.16) should be convoluted with diﬀerent parton distribution functions to obtain hadronic cross sections: in case of Eq. (2.16) we mustuse conventional MS PDFs whereas in case when the incoming charm quarks are massive aspecial set of PDFs is required. Nevertheless, since m c is just a collinear regulator, resultsfor short-distance hadronic cross sections should be the same, independent of whether onestarts with nearly massive or massless charm quarks. This requirement allows us to derive arelation between the “massive” and the MS PDFs. It reads f ( m ) a = ˆ O ab ⊗ f MS b , (2.18)where ˆ O ab ( z ) = δ ab δ (1 − z ) + (cid:16) α s π (cid:17) G ab ( z ) + . . . . (2.19)The computation that we just described allows us to determine the coeﬃcient G cc ( z ).We ﬁnd G cc ( z ) = − ln (cid:18) µ m c (cid:19) P AP qq ( z ) + C F (cid:20) z − z (1 + 2 ln(1 − z )) (cid:21) + . (2.20)– 8 –e can also compute “oﬀ-diagonal”coeﬃcients G ab that involve charm quark and gluonPDFs; they are important for removing mass singularities that arise in g → c and c → g transitions. Computations proceed along the same lines as described above except that weemploy other partonic processes for the analysis. Namely, we derive a relation for g → c transition by considering a cg → Hc process in a theory where only Yukawa coupling ispresent and no direct ggH coupling is allowed. To derive a relation for c → g transition,we again consider cg → Hc process but now only allow for the ggH coupling. In bothcases only (quasi)-collinear singularities are present; this simpliﬁes the required computationssigniﬁcantly. We ﬁnd G cg ( z ) = − ln (cid:18) µ m c (cid:19) P AP qg ( z ) ,G gc ( z ) = − (cid:20) ln (cid:18) µ m c (cid:19) − z ) − (cid:21) P AP gq ( z ) , (2.21)where P AP qg = T R [1 − z (1 − z )] , P AP gq = C F − z ) z . (2.22)The results for the functions G ab ( z ) reported in Eqs. (2.20,2.21) are important for thecalculation of NLO QCD corrections to the interference contribution to Higgs boson produc-tion in association with a charm jet. Indeed, as explained in the Introduction, to access theinterference, we need to start with the massive incoming charm quarks and carefully studythe massless limit. Since the charm mass serves as a collinear regulator, we are forced touse parton distribution functions f ( m ) i . We then use the relation Eq. (2.18) to express thesefunctions through the conventional MS PDFs and, in doing so, remove logarithms of m c thatare associated with the radiation by the incoming charm quarks. Because the interferencecontribution to pp → H + jet c involves a helicity ﬂip, standard collinear logarithms associatedwith initial state emissions are not the only logarithms of the charm mass that appear in thecross section. We elaborate on this statement in the next section. The discussion in the previous section shows that the dependence on m c disappears from hardcross sections provided that conventional factorization formulas for soft and quasi-collinearsingularities, Eqs. (2.3,2.6), hold true. However, since the interference contribution requiresa helicity ﬂip, its soft and quasi-collinear limits are diﬀerent from the conventional ones. Aswe explain below, such limits can still be described by simpler matrix elements but thesematrix elements do not always correspond to processes with reduced multiplicities of ﬁnalstate particles.To discuss and illustrate these subtleties in more detail, consider the process c ( p ) + g ( p ) −→ H ( p H ) + c ( p ) + g ( p ) . (3.1)– 9 –he ﬁrst point that needs to be emphasized is that, if we work with a ﬁnite charm mass, true soft and collinear limits of the process in Eq. (3.1) are, in fact, conventional. Thesecontributions can be extracted and combined with the virtual corrections to cg → Hc andrenormalized gluon parton distribution function giving a ﬁnite result for the partonic crosssection. Such a procedure is identical to what is usually done in NLO QCD computations[24, 25, 27] that are traditionally performed using dimensional regularization for soft andcollinear divergences. However, cancellation of “true” infrared and collinear divergences doesnot tell us anything about non-analytic dependence of partonic cross sections on the charmmass that we need to extract before taking the m c → m c → explicitly with massive charm quarks, for all relevant partonic processes, cg → Hcg , gg → Hc ¯ c , cq → Hcq , cc → Hcc , c ¯ c → Hc ¯ c . In that sense, we do not attempt todevelop an understanding of infrared and quasi-collinear factorization for generic processesthat involve a helicity ﬂip. However, to illustrate main diﬀerences with the conventionalcollinear factorization we discuss an emission of a collinear gluon oﬀ an incoming charmquark in case of the interference contribution in some detail.To this end, we consider the process in Eq. (3.1) in the quasi-collinear limit p · p ∼ m c .To describe this limit, we divide the amplitude for the full process cg → Hcg into two parts M = M sing + M ﬁn , (3.2)where M sing refers to diagrams where the gluon is emitted oﬀ the incoming quark with themomentum p and M ﬁn refers to the remaining diagrams. The ﬁrst contribution ( M sing ) issingular in the p · p ∼ m c → M ﬁn ) is not.Upon squaring the amplitude, Eq. (3.2), and summing over polarizations of initial andﬁnal state particles, we obtain (cid:88) pol |M| = (cid:88) pol |M sing | + (cid:88) pol (cid:16) M sing M † ﬁn + h . c . (cid:17) + . . . , (3.3)where the ellipsis stands for contributions that are ﬁnite in the p · p ∼ m c → not the case for the interference contributions considered in this paper.To extract the quasi-collinear singularities from the square of the amplitude in Eq. (3.3)we need to analyze the quasi-collinear kinematics; for this analysis there is no diﬀerencebetween helicity-conserving and helicity-ﬂipping contributions. Indeed, following the stan-dard approach, we re-write the four-momentum of the incoming charm quark and the four-momentum of the emitted gluon through massless momenta ˜ p and p and ﬁnd p = (cid:18) − m c s (cid:19) ˜ p + (cid:18) m c s (cid:19) p + O ( m c ) , p = (1 − z )˜ p + yp + p , ⊥ . (3.4)– 10 –n the above equation, s = 2˜ p · p and p , ⊥ · ˜ p = p , ⊥ · p = 0. We use the on-shell condition p = 0 and obtain y = − p , ⊥ (1 − z ) s . (3.5)It follows that 2 p · p ≈ s (cid:18) (1 − z ) m c s + y (cid:19) = 11 − z (cid:0) (1 − z ) m c − p , ⊥ (cid:1) . (3.6)We conclude that the kinematic region where p , ⊥ ∼ m c provides unsuppressed contributionsto the cross section in the m c → M sing = − g s t a i c i ˆ A i c sing ˆ p − ˆ p + m c p · p ) ˆ (cid:15) u ( p ) , (3.7)and use the decomposition of the four-momenta p , given in Eq. (3.4) to obtain M sing = g s t a i c i p · p ) ˆ A i c sing (cid:20) − p , ⊥ · (cid:15) − z + ˆ (cid:15) ( m c (1 − z ) + ˆ p , ⊥ + κ ˆ p ) (cid:21) u ( p ) . (3.8)In Eq. (3.8), we introduced a parameter κ deﬁned as κ = − m c (1 − z ) s − p , ⊥ s (1 − z ) . (3.9)We note that in deriving Eq. (3.8) we have used p · (cid:15) = 0 and the gauge ﬁxing condition p · (cid:15) = 0.The result for the singular contribution, Eq. (3.8), is generic; it does not distinguishbetween helicity-conserving and helicity-ﬂipping contributions. However, it is easy to see thatthere is a diﬀerence between the two. For example, since the helicity-ﬂipping contributionrequires one additional power of m c , one can convince oneself that a combination of termslabeled as κ in Eq. (3.9) may contribute to the collinear limit of the interference but it cannotcontribute to the collinear limit of the helicity-conserving amplitudes.Hence, performing standard manipulations and paying attention to subtleties indicatedabove, we obtain the contribution to the interference that is non-analytic in the m c → p · p → Int (cid:2) M (1 c , g ; 3 c , g ) (cid:3) = g s (cid:34) (cid:18) P qq ( z )( p · p ) − C F m c z ( p · p ) (cid:19) Int (cid:20) |M ( z c , g ; 3 c ) | z (cid:21) − C F m c (1 − z ) z ( p · p ) Int (cid:104) Tr (cid:104) ˆ A i c ( z · ˜1 c , g ; ˜3 c ) ˆ A i c , † ( z · ˜1 c , g ; ˜3 c ) (cid:105)(cid:105) (cid:35) + g s C F (1 − z ) m c p · p ) Int (cid:104) Tr (cid:104) ˆ p A i c ( z · ˜1 c , g ; ˜3 c ) ˆ A i c , † ﬁn (˜1 c , g ; ˜3 c , (1 − z ) · ˜1 g )ˆ (cid:15) + h . c . (cid:105)(cid:105) . (3.10)– 11 –t is understood that Int[ ... ] extracts the interference contributions from the relevant quanti-ties; sums over colors and polarizations are implicit. The quantities A are related to ampli-tudes M in the following way M (˜1 c , g ; ˜3 c ) = ˆ A i (˜1 c , g ; ˜3 c ) u (˜ p ) , M ﬁn (˜1 c , g ; ˜3 c , g ) = ˆ A i ﬁn (˜1 c , g ; ˜3 c , g ) u (˜ p ) . (3.11)They have to be computed in the m c = 0 limit. It is instructive to discuss the origin of the diﬀerent terms in Eq. (3.10). The ﬁrst termon the right-hand side of Eq. (3.10) contains the leading-order interference multiplied withthe standard massive collinear splitting function; this is the conventional quasi-collinear limitapplied to the interference. If only this term were present in Eq. (3.10), there would beno diﬀerences in the collinear factorization between helicity-changing and helicity-conservingcontributions.The second and the third terms on the right-hand side of Eq. (3.10) are new structures;they appear because the required helicity ﬂip can occur on the external charm quark lines .To illustrate this point, we square the ﬁrst equation in Eq. (3.11), sum over polarizations andﬁnd (cid:88) spins |M (1 c , g ; 3 c ) | = Tr (cid:104) (ˆ p + m c ) ˆ A i c ˆ A i c , † (cid:105) . (3.12)The interference requires a helicity ﬂip that is facilitated by a single mass insertion. The aboveequation shows that this mass insertion can occur either in the (ˆ p + m c ) density matrix ofthe external quark or “inside” the ˆ A ˆ A † term. The structure that appears in the second termon the right hand side of Eq. (3.10) originates from the mass term in the density matrix.Once the mass term is extracted, the rest can be computed in the massless approximation.We ﬁnd Int (cid:104) Tr (cid:104) (ˆ p + m c ) ˆ A i c (1 c , g ; 3 c ) ˆ A i c , † (1 c , g ; 3 c ) (cid:105)(cid:105) → m c Int (cid:104) Tr (cid:104) ˆ A i c (˜1 c , g ; ˜3 c ) ˆ A i c , † (˜1 c , g ; ˜3 c ) (cid:105)(cid:105) . (3.13)The last term on the right hand side in Eq. (3.10) describes a quasi-collinear singularitythat originates from the interference of singular and regular contributions in Eq. (3.3); itis particular to the helicity-ﬂipping case and does not appear in the conventional collinearlimits. As a consequence, this contribution still depends on the part of the reduced matrixelement of the original 2 → We remind the reader that the notation ˜ i implies that a light-cone four-momentum of a particle i must beused in the computation. – 12 –nterference contribution as they do not follow canonical pattern and cannot be removed by atransition to MS parton distribution functions. In addition, we also ﬁnd that the interferencecontributions exhibit quasi-soft quark singularities that also lead to logarithms of the charmmass. Although it would be interesting to understand factorization of mass singularities inprocesses with the helicity ﬂip from a more general perspective, our strategy for now is to explicitly compute all relevant contributions within ﬁxed-order perturbation theory extractingall non-analytic m c -dependent terms along the way. Additional details of our approach canbe found in several appendices. In this section we brieﬂy describe calculation of one-loop and real emission contributions tothe interference part of the cg → Hc process. We begin with the discussion of the virtualcorrections.We compute the relevant one-loop diagrams keeping charm-quark masses ﬁnite. We em-ploy the standard Passarino-Veltman reduction [28] to express the cg → Hc amplitude interms of one-loop scalar integrals. After computing the one-loop contribution to the inter-ference, we expand the expression around m c = 0 and keep the leading O ( m c ) term in thisexpansion. The one-loop amplitudes contain ultraviolet and infrared singularities. The ul-traviolet singularities are removed by the renormalization. We closely follow the discussion inAppendix A of Ref. [26] where many of the required renormalization constants are presented.Similar to Ref. [26], we renormalize the charm-quark mass in the on-shell scheme but employthe MS renormalization for the Yukawa coupling constant. In addition to the discussion inthat reference, we require the one-loop renormalization constant of the eﬀective ggH vertexthat we take from Ref. [29]. After the ultraviolet renormalization is performed, the cg → Hc amplitude still contains 1 /(cid:15) poles of infrared origin. These poles satisfy the Catani’s formula[30] and cancel with similar poles in real emission contributions to the partonic cross section.According to the discussion in Section 3, factorization of quasi-collinear and quasi-softsingularities in the interference contribution is not canonical. This implies that even if wetake the m c → m c → We use

FeynCalc [31, 32] for a cross-check of our computation. We employed the Package-X program [33] for numerical checks of scalar integrals and their m c → – 13 –DFs, and the virtual corrections. The only diﬀerence with respect to the canonical procedurefor NLO QCD computations is that in our case the subtraction terms are directly obtainedfrom the squared amplitude and are not written in terms of easily recognizable universalfunctions; see Section 3 and Appendix A for further details. In particular, even contributionsthat are enhanced by two powers of a logarithm of the charm mass, O ( α s ln ( m c )), do notappear to be proportional to the leading-order interference contribution to the cross section. To present numerical results we consider proton-proton collisions at 13 TeV. We take M H =125 GeV for the Higgs-boson mass and m c = 1 . m c ( M H ) = 0 .

81 GeV. We use

NNPDF31 lo as 0118 and

NNPDF31 nlo as 0118 parton distribution functions [34, 35]for leading and next-to-leading order computations, respectively. The value of the strongcoupling constant α s is calculated using dedicated routines provided with NNPDF sets.To deﬁne jets we use standard anti- k ⊥ algorithm with ∆ R = 0 .

4; charm jets are requiredto contain at least one c or ¯ c quark. For numerical computations we require at least one charmjet with p t,j >

20 GeV and | η j | < .

5. Moreover, we demand that the charm parton insidethe charm jet carries at least 75% of the jet’s transverse momentum. The latter requirementremoves kinematic cases where a soft charm is clustered together with a hard gluon into acharm jet in spite of a large angular separation between the two. Since, as we explainedearlier, our calculation is logarithmically sensitive to soft emissions of charm quarks , deﬁningcharm jets with an additional cut on the charm quark transverse momentum allows us toavoid jet-algorithm dependent logarithms of m c that may appear otherwise. We note thatwe apply all these requirements even in the subtraction terms where c and ¯ c momenta arecomputed in the collinear and/or soft approximations.We start by presenting ﬁducial cross sections for the three terms in Eq. (1.2) separately.Central values for all the cross sections presented below correspond to the renormalizationand factorization scales set to µ F = µ R = µ = M H ; subscripts and superscripts indicate shiftsin central values if µ = M H / µ = 2 M H are used in the calculation. At leading order,we ﬁnd σ LO ggH = 176 . +47 . − . fb , σ LOYuk = 21 . +1 . − . fb , σ LOInt = − . +0 . − . fb , (5.1)for the ggH -dependent cross section, the Yukawa-dependent cross section and the interference,respectively. In Figure 2 we show the comparison between µ = M H cross sections and theinterference in dependence on the cut of the charm jet transverse momentum. We observethat the ratio of the interference to the Yukawa-dependent contribution is about ten percentfor all values of the p t,j -cut. We use program

RunDec [36, 37] to compute the value of the running charm quark mass. If more than one c or ¯ c parton is clustered into a jet, we apply this requirement to the hardest of them. – 14 – p min t,j c [ GeV ]10 − | σ L O | [ f b ] LO (ggH)LO (Yuk)LO (Int) (a)

Fiducial cross section (LO)

20 40 60 80 100 120 p min t,j c [ GeV ]8910111213 − σ L O ( I n t) / σ L O ( Y u k ) · % (b) Interference/Yukawa ratio

Figure 2:

Leading-order cross sections computed for diﬀerent values of the charm-jet p t -cut.We show Yukawa-like (green) and ggH -like (orange) contributions as well as the absolute value of the interference (blue). In the right pane the ratio of the interference and the Yukawa-likeﬁducial cross section are shown.∆ σ NLO [ fb ] cg cq gg cc c ¯ c PDF sum const − .

63 0 .

13 2 .

33 0 . − .

01 0 .

11 0 . L . − − . − .

04 0 .

01 1 . − . L − . − .

66 0 . − . − . .

54 0 . − . − . − .

08 1 .

76 1 . Table 1:

The NLO QCD corrections to the interference split according to partonic channels.The results are given in femtobarns. The column marked “PDF” refers to the PDF-schemechange discussed in Section 2. For each partonic channel we show O ( L ), O ( L ) and O ( L )contributions where L = ln( M H /m c ).At next-to-leading order the ﬁducial cross section for the interference term becomes σ NLOInt = − . +0 . − . fb . (5.2)It follows from Eqs. (5.1,5.2) that the NLO QCD corrections decrease the absolute value ofthe leading-order interference by about ﬁfty percent. The scale uncertainty appears to bereduced by about a factor 2. The NLO result for the interference is outside the leading-orderscale-uncertainty interval, c.f. Eqs. (5.1,5.2), emphasizing the fact that the appearance ofthe logarithms of the charm mass in NLO QCD corrections to the interference makes thescale variation uncertainty of the leading-order result a very poor indicator of the theoreticaluncertainty in this case. – 15 –t is instructive to separate the NLO contributions to the interference into parts that areindependent of m c and parts that are logarithmically enhanced for all the partonic channels.The relevant results are shown in Table 1. We ﬁnd that the largest contribution at NLO comesfrom the gluon-gluon channel which is enhanced by the large gluon luminosity. Also, thecharm-gluon ( cg ) and the charm-quark channels ( cq ) provide relatively large contributions. Note that the ( cq ) channel is free of logarithmic contributions since there are no singular limitsthat involve charm quarks. Contributions related to the PDF transformation do not featurethe double-logarithmic part since the O (ln m c ) terms originate exclusively from soft-collinearlimits that involve c -quarks.It follows from Table 1 that double-logarithmic terms and single-logarithmic terms pro-vide nearly equal, but opposite in sign, contributions to the NLO QCD interference. Thiscancellation between terms with diﬀerent parametric dependence on m c should be consideredas an artifact but it does emphasize that studying only the leading logarithmic O (ln m c )contribution in this case is insuﬃcient for phenomenology. We also note that the O (ln m c )term in the cg channel is quite small reﬂecting the fact that there is a very strong – but incom-plete – cancellation between double-logarithmic contributions to real and virtual correctionsin this case. Finally, we emphasize that it is unclear to what extent these various cancella-tions persist in higher orders; for this reason, a resummation of charm-mass logarithms forthe interference contribution is desirable.We continue with the discussion of kinematic distributions. We focus on the transversemomentum and the rapidity distributions of Higgs bosons in the interference contribution to pp → Hc cross section. They are shown in Figure 3. We ﬁrst discuss the transverse momentumdistribution, Figure 3a; when interpreting this ﬁgure it is important to recall that the absolute value of both LO and NLO distributions is plotted there and that the LO distribution isalways negative . We observe in Figure 3a that the leading-order distribution (blue) is largeand negative at small p t,H ; as p t,H increases, the distribution goes to zero. The NLO QCDcorrections aﬀect the shape of the p t,H distribution. Indeed, a sharp edge at p t,H = 20 GeV,present at leading order, gets smeared at NLO. At moderate values of transverse momenta p t,H ∼

60 GeV the K -factor is equal to one, while there is a large O (+50%) reduction at p t,H ∼

100 GeV. Second, at p t,H ∼

150 GeV the NLO distribution goes through zero andbecomes positive for larger values of p t,H . Asymptotically, at even higher p t,H the LO andNLO distributions appear to be equal in absolute values but opposite in sign. Of course, allthis happens at such high values of p t,H that are irrelevant for phenomenology, but it is quitea peculiar feature nevertheless. Here, by “quark” we mean any quark of a ﬂavor other than c . There is a subtlety related to the b -quarkcontribution because b -quarks have stronger interactions with Higgs bosons as compared to charm quarks. Suchcontributions can, presumably, be dealt with using b anti-tagging. When presenting results for the interferencewe decided to include contributions of bottom quarks, setting bottom Yukawa coupling to zero, but we didcheck that the ﬂavor-excitation topologies with b in the initial state change the results for ( cq ) channel byabout three percent only. “Reduction” in this case means that the distribution becomes less negative. – 16 – − − − − − | d σ / d p t , H | [ f b / G e V ] LO (Int)NLO (Int)0 100 200 300 400 500 p t,H [ GeV ] − R a t i o (a) Higgs-boson transverse momentum − . − . − . − . − . . d σ / d | y H | [ f b ] LO (Int)NLO (Int)0 1 2 3 4 | y H | . . . R a t i o (b) Higgs-boson rapidity

Figure 3:

The transverse momentum and rapidity distributions of the Higgs boson calculatedat LO (blue) and NLO (red) for central scale choice. We only consider the interferencecontribution. We note that the absolute value of d σ Int / d p t,H is displayed in the left panel.This implies that this distribution actually changes sign at around p t,H ∼

150 GeV. Thelower panels show ratios to the LO interference contribution.Compared to Higgs transverse momentum distribution, the rapidity distribution of theHiggs boson in the interference contribution is much less volatile. Indeed, it follows fromFigure 3b that the diﬀerence between leading and next-to-leading-order distributions is well-described by a constant K -factor all the way up to | y H | ∼

2. Beyond this value of the rapidity,the NLO distribution goes to zero faster than the LO one.

Production of Higgs bosons in association with charm jets at the LHC is mediated by twodistinct mechanisms, one that involves the charm Yukawa coupling and the other one thatinvolves an eﬀective ggH vertex. Their interference involves a helicity ﬂip and, for this reason,vanishes in the limit of massless charm quarks.Since partonic cross sections are routinely computed for massless incoming partons andsince the charm quark appears in the initial state in the main process cg → Hc , it is interestingto understand how to circumvent the problem of having to deal with a massive parton in theinitial state and to provide reliable estimate of the interference contribution.We have addressed this problem by studying the m c → m c → O (ln m c ) contributions by– 17 –xpressing results through conventional MS parton distribution functions valid for masslesspartons. Nevertheless, given an unconventional behavior of the interference in quasi-soft andquasi-collinear limits, logarithms of the charm quark mass survive in the ﬁnal result for theNLO QCD corrections.We have found that the absolute value of the leading-order interference is reduced byabout ﬁfty percent once NLO QCD corrections are accounted for. This signiﬁcant but still“perturbatively acceptable” reduction is the result of a very strong cancellation between termsthat involve double and single logarithms of the charm quark mass. We have observed that theNLO QCD corrections to the interference are kinematics-dependent and may change shapesof certain kinematic distributions in a signiﬁcant way.Higgs boson production in association with a charm jet is a promising way to studycharm Yukawa coupling at the LHC [12]. The interference contribution, that is estimated to beabout 10 percent of the Yukawa contribution at leading order, could have been perturbativelyunstable given the required helicity ﬂip and an unconventional pattern of quasi-soft andquasi-collinear limits. We addressed this question by performing a dedicated NLO QCDcomputation for the interference term and did not ﬁnd a strong indication that this mightbe the case. Nevertheless, the moderate size of the NLO QCD corrections is the consequenceof a very strong cancellations between double and single logarithms of the charm mass. Itis unclear if this cancellation persists in higher orders. Hence, resummation of O (ln m c )-enhanced terms for this process is quite desirable. Acknowledgments

This research is partially supported by the Deutsche Forschungsgemein-schaft (DFG, German Research Foundation) under grant 396021762 - TRR 257.

A Extraction of the O (ln m c ) -enhanced contributions in the real corrections In this appendix we describe a procedure to extract O (ln m c ) contributions to real emissioncorrections. They arise because of the quasi-singular behavior of real emission amplitudes inthe soft or collinear limits involving charm quarks. The potential singularities in these limitsare regulated by the charm mass leading to an appearance of O (ln m c )-enhanced terms whenintegrated over relevant phase spaces. To extract logarithms of m c , we subtract approximateexpressions from exact matrix elements that make the diﬀerence integrable in the m c → m c .As an example, we consider the gluon-gluon partonic channel, i.e. g ( p ) + g ( p ) −→ H ( p H ) + c ( p ) + ¯ c ( p ) , (A.1)and discuss the extraction of O (ln m c ) terms in detail. This channel is suitable for sucha discussion since, if the charm-quark mass is kept ﬁnite, it is free of soft and collineardivergences. Hence, all relevant contributions can be computed numerically for small but– 18 –nite m c , and used to validate formulas where logarithms of m c have been extracted and m c → massive charm quarks is new and requires an explanation.We focus on the interference contribution between the Yukawa-like and the ggH -likeproduction mechanisms in the process Eq. (A.1). The interference term is non-zero only ifhelicity ﬂip on the charm line occurs. Furthermore, the presence of such a helicity ﬂip causesthe usual factorization formulas to break down and the singular limits need to be explicitlyanalyzed. We note that, thanks to the symmetry of the squared amplitude for the processin Eq. (A.1) under the exchange of c and ¯ c , we can consider only the case where ¯ c quarkbecomes soft or quasi-collinear to one of the other partons. The case when both c and ¯ c become unresolved is impossible since we require a charm jet in the ﬁnal state.The quasi-singular limits which appear in this channel are related to the soft-quark limit S with E ∼ m c and the three collinear limits C i with i = 1 , , p · p i ) ∼ m c .Performing an iterative subtraction of these singular limits, we ﬁnd (cid:104) F LM (1 g , g ; 3 c , ¯ c ) (cid:105) = (cid:88) i =1 (cid:104) (1 − C i )(1 − S ) ω ( i )123 F LM (1 g , g ; 3 c , ¯ c ) (cid:105) + (cid:104) C i (1 − S ) ω ( i )123 F LM (1 g , g ; 3 c , ¯ c ) (cid:105) + (cid:104) S F LM (1 g , g ; 3 c , ¯ c ) (cid:105) , (A.2)where the ﬁrst term on the right-hand side denotes the fully-regulated contribution and thesecond and third terms are the collinear and the soft integrated subtraction terms. Thefactors ω ( i )123 are the weights that describe various collinear sectors. They read ω ( i )123 = 1 ρ i · (cid:20) ρ + 1 ρ + 1 ρ (cid:21) − , (A.3)with ρ i = 1 − cos θ i . We note that, since all m c → m c → O (ln m c )-terms and constants which survive the m c → .1 Integration of the soft-quark subtraction terms Consider the soft-quark subtraction term (cid:104) S F LM (1 g , g ; 3 c , ¯ c ) (cid:105) . (A.4)To compute it, we need to know the behavior of the amplitude in the limit p ∼ m c → p .Although, normally, soft (gluon) emissions factorize into a product of an eikonal factorand a tree-level matrix element squared, a similar formula for soft-quark emission, relevantfor helicity-ﬂipping processes, does not exist. Hence, we determine the soft-quark limit of theinterference by studying an explicit expression of the amplitude for the process in Eq. (A.1)in the limit p ∼ m c →

0. We ﬁnd S Int (cid:2) |M (1 g , g ; 3 c , ¯ c ) | (cid:3) ∼ (2 C F − C A )( p · p )( p · p )( p · p ) F ( p , p , p )+ C A ( p · p )( p · p )( m c + p · p ) F ( p , p , p )+ C A ( p · p )( p · p )( m c + p · p ) F ( p , p , p ) , (A.5)where functions F ij depend on the momenta p , p and p only . We emphasize that thesefunctions are diﬀerent from the leading-order interference contribution. The massless limit, m c →

0, can be now taken everywhere except for the eikonal factors and the phase-spacemeasure of the unresolved parton [d p ].We write the integrated soft subtraction term as follows (cid:104) S F LM (1 g , g ; 3 c , ¯ c ) (cid:105) = (2 C F − C A ) (cid:104) F ( p , p , p ) · (cid:90) [d p ]( p · p )( p · p )( p · p ) (cid:105) + C A (cid:104) F ( p , p , p ) · (cid:90) [d p ]( p · p )( p · p )( m c + p · p ) (cid:105) + C A (cid:104) F ( p , p , p ) · (cid:90) [d p ]( p · p )( p · p )( m c + p · p ) (cid:105) , (A.6)where (cid:104) . . . (cid:105) denotes the phase space integration and the relevant soft integrals can be found inAppendix B. We stress that soft integrals are ﬁnite in four dimensions since they are naturallyregulated by the charm-quark mass m c . A.2 Integration of the quasi-collinear subtraction terms

In this subsection, we discuss how to deﬁne and compute the soft-subtracted quasi-collinearlimits of the interference contribution using the process in Eq. (A.1) as an example. We focuson the sector 43 where c and ¯ c become collinear to each other. The relevant quantity reads (cid:104) C (1 − S ) F LM (1 g , g ; 3 c , ¯ c ) (cid:105) . (A.7) We note that weight factors introduced in Eq.(A.3) do not appear in the collinear limits. – 20 –o proceed further, we split the above formula into collinear and soft-collinear terms (cid:104) C (1 − S ) F LM (1 g , g ; 3 c , ¯ c ) (cid:105) = (cid:104) C F LM (1 g , g ; 3 c , ¯ c ) (cid:105) − (cid:104) C S F LM (1 g , g ; 3 c , ¯ c ) (cid:105) . (A.8)We ﬁrst discuss the collinear subtraction term (cid:104) C F LM (1 g , g ; 3 c , ¯ c ) (cid:105) deﬁned as follows (cid:104) C F LM (1 g , g ; 3 c , ¯ c ) (cid:105) = (cid:88) i =1 (cid:90) [d p H ][d p ][d p ](2 π ) δ ( p − p H − p − ¯ p ) × p + p ) p i · p p i · ¯ p C i (1 g , g , c , ¯4 ¯ c ) , (A.9)In the above equation, the functions C , depend on the momenta p , p , p and ¯ p . The barover momentum p indicates that the relevant collinear limit has been taken, i.e. C p = ( E , β (cid:126)n ) ≡ ¯ p . (A.10)Note that in Eq. (A.10) β = (cid:112) − m c /E is the velocity of ¯ c and (cid:126)n is a unit vector pointingin the direction of momentum (cid:126)p . We note that in Eq. (A.9) the massless limit m c → C , are regular in the soft-quark limit, E → m c ∼ O (ln m c ) terms arising from Eq. (A.9) and take the masslesslimit after that. To do so, we add and subtract the soft limits of the functions C i C i (1 g , g , c , ¯4 ¯ c ) = (cid:2) C i (1 g , g , c , ¯4 ¯ c ) − C i, soft (1 c , g , c ) (cid:3) + C i, soft (1 c , g , c ) , (A.11)where C i, soft (1 c , g , c ) = C i (1 g , g , c , regulated integral that contains the expression in the square bracket inEq. (A.11) and the soft part. In the regulated part, the soft divergence at E = 0 has beenregularized. This implies that, after integrating 1 / ( p + p ) over the relative angle between p and p and extracting logarithms of m c from this angular integral, we can set m c to zeroeverywhere else right away. We obtain (cid:104) C F LM (1 g , g ; 3 c , ¯ c ) (cid:105) reg = 1(2 π ) (cid:88) i =1 (cid:90) [d p H ][d p ](2 π ) δ ( p − p H − p ) (cid:90) z d z − z (cid:104) C i (1 g , g , z , (1 − z )34) − C i, soft (1 c , g , z (cid:105) × (cid:104) ln(2 E /m c ) + ln(1 − z ) + ln( z ) (cid:105) , (A.12)where we have used the fact that in m c → p = zp and ¯ p = (1 − z ) p for p = 0. – 21 –e will now discuss the soft part of the collinear subtraction term. It reads (cid:104) C F LM (1 g , g ; 3 c , ¯ c ) (cid:105) soft = (cid:88) i =1 (cid:90) [d p H ][d p ][d p ](2 π ) δ ( p − p H − p − ¯ p ) × p + p ) p i · p p i · ¯ p C i, soft (1 c , g , c ) . (A.13)We emphasize that this term still contains soft singularity and, for this reason, the m c → (cid:104) C S F LM (1 g , g ; 3 c , ¯ c ) (cid:105) , c.f. Eq. (A.8); if this is done, the required com-putations simplify signiﬁcantly.The soft-collinear integrated subtraction term in sector 43 reads (cid:104) C S F LM (1 g , g ; 3 c , ¯ c ) (cid:105) = (cid:88) i =1 (cid:90) [d p H ][d p ][d p ](2 π ) δ ( p − p H − p ) × C A F i ( p , p , p )( p + p ) p i · p p i · ¯ p . (A.14)To derive this result we used the soft-limit of the interference amplitude reported in Eq. (A.5).We emphasize that, since the soft operator is present on the left hand side in the aboveequation, the soft-quark momentum p is removed from the energy-momentum conservingdelta-function. Moreover, since C i, soft (1 c , g , c ) = 2 C A F i ( p , p , p ) , (A.15)the two integrals in Eqs. (A.13,A.14) appear to be the same up to the argument of thedelta-functions. We combine the two integrals and ﬁnd (cid:104) C F LM (1 g , g ; 3 c , ¯ c ) (cid:105) soft − (cid:104) C S F LM (1 g , g ; 3 c , ¯ c ) (cid:105) == (cid:88) i =1 (cid:90) [d p H ][d p ][d p ](2 π ) (cid:104) δ ( p − p H − p − ¯ p ) − δ ( p − p H − p ) (cid:105) × C A F i ( p , p , p )( p + p ) p i · p p i · ¯ p . (A.16)To proceed further, we note that it is straightforward to integrate over directions of thequark with momentum p but integration over its energy is more involved. It is convenientto split the E integration into two regions by introducing an auxiliary parameter σ E − σ ) + Θ( σ − E ) . (A.17)We choose σ to satisfy the following inequality m c (cid:28) σ (cid:28) E . For small energies, E < σ (cid:28) E , we can drop the momentum ¯ p from the energy momentum conserving delta-functionwhich leads to (cid:104) δ ( p − p H − p − ¯ p ) (cid:2) Θ( E − σ ) + Θ( σ − E ) (cid:3) − δ ( p − p H − p ) (cid:105) == (cid:104) δ ( p − p H − p − ¯ p ) − δ ( p − p H − p ) (cid:105) Θ( E − σ ) + O ( σ/E ) . (A.18)– 22 –his relation implies that the integrand in Eq. (A.16) is non-vanishing only in the high-energydomain where E > σ (cid:29) m c and, therefore, the limit m c → (cid:104) C F LM (1 g , g ; 3 c , ¯ c ) (cid:105) soft − (cid:104) C S F LM (1 g , g ; 3 c , ¯ c ) (cid:105) == C A (2 π ) (cid:88) i =1 (cid:90) [d p H ][d p ](2 π ) δ ( p − p H − p ) (cid:40) z max (cid:90) z min d z − z ln (cid:18) E (1 − z ) zm c (cid:19) (cid:104) zθ ( z ) F i ( p , p , zp ) − F i ( p , p , p ) (cid:105) + z max (cid:90) z min d z − z ln (cid:0) (2 − z ) z (cid:1) F i ( p , p , p ) (cid:41) . (A.19)To arrive at Eq. (A.19) we introduced the four-momentum p = p + ¯ p and a variable z such that p = zp in terms that contain the delta-function δ ( p − p H − p − ¯ p ). In termsthat contain the delta-function δ ( p − p H − p ), we set (1 − z ) = E /E , rename p into p and set σ →

0. The lower integration boundary z min is given by z min = 1 − E max /E < (cid:104) C (1 − S ) F LM (1 g , g ; 3 c , ¯ c ) (cid:105) is given by a sumof expressions in Eqs. (A.12,A.19). We describe a numerical check of validity of this result inthe following section. A.3 Numerical checks

Since the cross section of the gluon-gluon channel, Eq. (A.1), is ﬁnite as long as we keepthe non-zero charm mass, analytical results derived in the previous section can be checkednumerically by computing σ gg → Hc ¯ c explicitly for small values of the charm mass without anyapproximation.The comparison is shown in Figure 4. We use ﬁducial cuts described in the main text andcompare hadronic contributions to the interference for gg partonic channel computed in twodiﬀerent ways. Green points (rectangles) show the results of the computation without anyapproximation, i.e. by directly integrating the matrix element squared. Blue points (circles)show the results of the computation that relies on the expansion around m c → m c . Theupper panel of Figure 4 shows the absolute values of the interference cross section in the gg partonic channel obtained with the two methods, while their diﬀerence is shown in the lowerpanel. We see a better and better agreement between the two results as we mover to smallerand smaller values of the charm-quark mass. This indicates that the m c -dependence of theinterference contribution is properly reconstructed.– 23 – − − − charm-quark mass [GeV]0510 R e l a t i v e d i ﬀ e r e n ce [ % ] ( σ real − σ real ) /σ real ( σ rec − σ real ) /σ real − − C r o ss - s ec t i o n [ f b ] σ rec σ real Figure 4:

The cross section of the gg → Hc ¯ c process calculated by a direct integration of thematrix element with non-zero charm-quark mass, σ real (green rectangles), and reconstructedusing procedure described in previous subsections, σ rec (blue circles). We employ the sameparameters and kinematic constraints as in the main text. B Soft-quark integrals

In this section we list integrals that are required for the integrated soft-quark subtractionterms. We need a number of integrals depending on the type and conﬁguration of the emitters p a and p b as well as the propagator appearing in the eikonal factor.We note that we are interested only in the terms that contain logarithms of the charm-quark mass and constant terms, but we drop all power-suppressed terms which vanish in the m c → d = 4 dimensions since all singularities arenaturally regulated by the charm-quark mass.The phase-space measure for a massive-quark emission, p = m c , reads[d p ] = k d k E dΩ (3) (2 π ) Θ( E max − E )Θ( E − m c ) , (B.1)where k is the length of (cid:126)p momentum, dΩ (3) denotes angular integration and E max is theusual energy cutoﬀ of the nested soft-collinear subtraction scheme. In the remaining partof this section, we list soft-quark integrals that are needed to obtain integrated soft-quarksubtraction terms, see Section A.1 for details. Similar expressions to those in Section A.1 can be derived for other partonic channels featuring soft-quark – 24 – wo massless emitters:

Two emitters a, b have four-momenta p a = E a (1 , (cid:126)n a ) and p b = E b (1 , (cid:126)n b ), respectively. Both four-momenta are light-like p a = p b = 0. Vectors (cid:126)n a and (cid:126)n b describe direction of ﬂight of the emitters; we refer to the opening angle between (cid:126)n a and (cid:126)n b as θ ab .The soft integral reads (cid:90) [d p ]( p a · p b )( p a · p )( p b · p ) = 1(2 π ) (cid:20) ln (2 s ab E max /m c ) − π

12 + 12 Li ( c ab ) (cid:21) , (B.2)where we used s ab = sin( θ ab /

2) and c ab = cos( θ ab / One massive and one massless emitters:

Two emitters a, b have four-momenta p a = E a (1 , (cid:126)n a ) and p b = E b (1 , β b (cid:126)n b ), respectively. They satisfy p a = 0 and p b = m c . We refer tothe opening angle between (cid:126)n a and (cid:126)n b as θ ab . We require three soft integrals of this type (cid:90) [d p ]( p a · p b )( p a · p )( m c + p b · p ) = 1(2 π ) (cid:20) ln (2 s ab E max /m c ) − π

12+ Li ( − E max /E b ) + 12 Li ( c ab ) (cid:21) , (cid:90) [d p ]( p a · p b )( p a · p )( p b · p ) = 1(2 π ) (cid:20) ln (2 s ab E max /m c ) + 14 Li ( − E max2 /E b ) − π

12 + 12 Li ( c ab ) (cid:21) , (cid:90) [d p ]( p a · p b )( p a · p )( m c − p b · p ) = 1(2 π ) (cid:20) − ln (2 s ab E max /m c ) + Li (1 − E max /E b ) −

12 Li ( c ab ) + ln( E max /E b ) ln( E max /E b − − π (cid:21) , (B.3)where we used s ab = sin( θ ab /

2) and c ab = cos( θ ab / Two massive emitters:

Two emitters a, b have four-momenta p a = E a (1 , β a (cid:126)n a ) and p b = E b (1 , β b (cid:126)n b ), respectively. They satisfy p a = p b = m c . We refer to the opening angle between (cid:126)n a and (cid:126)n b as θ ab . In this case, we use E max = E a . We ﬁnd (cid:90) [d p ]( p a · p b )( m c − p a · p )( m c + p b · p ) = 1(2 π ) (cid:20) − ln (2 s ab E max /m c ) − π −

12 Li ( c ab ) − Li ( − E max /E b ) (cid:21) , (cid:90) [d p ]( p a · p b )( m c − p a · p )( m c − p b · p ) = 1(2 π ) (cid:20) ln (2 E max /m c ) + π (cid:21) , (B.4)where we used s ab = sin( θ ab /

2) and c ab = cos( θ ab / singularities, i.e. cc and c ¯ c . – 25 – eferences [1] M. Aaboud et al. [ATLAS], Phys. Lett. B (2018), 59-86[2] A. M. Sirunyan et al. [CMS], Phys. Rev. Lett. (2018) no.12, 121801[3] M. Aaboud et al. [ATLAS], Phys. Rev. D (2019), 072001[4] A. M. Sirunyan et al. [CMS], JHEP (2019), 093[5] G. Aad et al. [ATLAS], Phys. Lett. B (2021), 135980[6] A. M. Sirunyan et al. [CMS], Phys. Rev. Lett. (2019) no.2, 021801[7] G. Perez, Y. Soreq, E. Stamou and K. Tobioka, Phys. Rev. D (2016) no.1, 013001[8] A. L. Kagan, G. Perez, F. Petriello, Y. Soreq, S. Stoynev and J. Zupan, Phys. Rev. Lett. (2015) no.10, 101802[9] T. Modak and R. Srivastava, Mod. Phys. Lett. A (2017) no.03, 1750004[10] M. K¨onig and M. Neubert, JHEP (2015), 012[11] F. Bishara, U. Haisch, P. F. Monni and E. Re, Phys. Rev. Lett. (2017) no.12, 121801[12] I. Brivio, F. Goertz and G. Isidori, Phys. Rev. Lett. (2015) no.21, 211801[13] A. A. Penin, Phys. Lett. B (2015), 69-72 [erratum: Phys. Lett. B (2015), 596-596;erratum: Phys. Lett. B (2017), 633-634][14] K. Melnikov and A. Penin, JHEP (2016), 172[15] Z. L. Liu and M. Neubert, JHEP (2020), 033[16] Z. L. Liu, B. Mecaj, M. Neubert and X. Wang, JHEP (2021), 077[17] E. Laenen, J. Sinninghe Damst´e, L. Vernazza, W. Waalewijn and L. Zoppi, [arXiv:2008.01736].[18] S. Catani, S. Dittmaier and Z. Trocsanyi, Phys. Lett. B (2001), 149-160[19] S. Frixione, JHEP (2019), 158.[20] F. Caola, K. Melnikov and R. R¨ontsch, Eur. Phys. J. C (2017) no.4, 248[21] F. Caola, K. Melnikov and R. R¨ontsch, Eur. Phys. J. C (2019) no.5, 386[22] F. Caola, K. Melnikov and R. R¨ontsch, Eur. Phys. J. C (2019) no.12, 1013[23] K. Asteriadis, F. Caola, K. Melnikov and R. R¨ontsch, Eur. Phys. J. C (2020) no.1, 8[24] S. Frixione, Z. Kunszt and A. Signer, Nucl. Phys. B (1996), 399-442[25] S. Frixione, Nucl. Phys. B (1997), 295-314[26] A. Behring and W. Bizo´n, JHEP (2020), 189[27] S. Catani and M. H. Seymour, Nucl. Phys. B (1997), 291-419 [erratum: Nucl. Phys. B (1998), 503-504][28] G. Passarino and M. J. G. Veltman, Nucl. Phys. B (1979), 151-207[29] R. Mondini, U. Schubert and C. Williams, JHEP (2020), 058[30] S. Catani, Phys. Lett. B (1998), 161-171 – 26 –

31] R. Mertig, M. Bohm and A. Denner, Comput. Phys. Commun. (1991), 345-359[32] V. Shtabovenko, R. Mertig and F. Orellana, Comput. Phys. Commun. (2020), 107478[33] H. H. Patel, Comput. Phys. Commun. (2015), 276-290[34] A. Buckley, J. Ferrando, S. Lloyd, K. Nordstr¨om, B. Page, M. R¨ufenacht, M. Sch¨onherr andG. Watt, Eur. Phys. J. C (2015), 132[35] R. D. Ball et al. [NNPDF], Eur. Phys. J. C (2017) no.10, 663[36] K. G. Chetyrkin, J. H. Kuhn and M. Steinhauser, Comput. Phys. Commun. (2000), 43-65[37] F. Herren and M. Steinhauser, Comput. Phys. Commun. (2018), 333-345(2018), 333-345