Next-to-leading-order QCD corrections to $e^+e^-\to H+\gamma$

The associated production of Higgs boson with a hard photon at lepton collider, i.e., $e^+e^-\to H\gamma$, is known to bear a rather small cross section in Standard Model, and can serve as a sensitive probe for the potential new physics signals. Similar to the loop-induced Higgs decay channels $H\to \gamma\gamma, Z\gamma$, the $e^+e^-\to H\gamma$ process also starts at one-loop order provided that the tiny electron mass is neglected. In this work, we calculate the next-to-leading-order (NLO) QCD corrections to this associated $H+\gamma$ production process, which mainly stem from the gluonic dressing to the top quark loop. The QCD corrections are found to be rather modest at lower center-of-mass energy range ($\sqrt{s}<300$ GeV), thus of negligible impact on Higgs factory such as CEPC. Nevertheless, when the energy is boosted to the ILC energy range ($\sqrt{s}\approx 400$ GeV), QCD corrections may enhance the leading-order cross section by $20\%$. In any event, the $e^+e^-\to H\gamma$ process has a maximal production rate $\sigma_{\rm max}\approx 0.08$ fb around $\sqrt{s}= 250$ GeV, thus CEPC turns out to be the best place to look for this rare Higgs production process. In the high energy limit, the effect of NLO QCD corrections become completely negligible, which can be simply attributed to the different asymptotic scaling behaviors of the LO and NLO cross sections, where the former exhibits a milder decrement $\propto 1/s$ , but the latter undergoes a much faster decrease $\propto 1/s^2$.


Introduction
After the historical discovery of the 125 GeV boson by the ATLAS and CMS collaborations at LHC in 2012 [1,2], a great amount of efforts have been devoted to unravelling its nature.Numerous evidences have been accumulating to indicate that this new particle is just the long-sought Higgs boson, which plays the pivotal role in mediating spontaneous electroweak symmetry breaking.To date, the measured couplings between the Higgs boson and heavy fermions/gauge bosons are compatible with the Standard Model (SM) predictions within 10% − 20% accuracy [3].
One of the central goals of contemporary high-energy physics is to precisely nail down the properties of the Higgs boson, and to search for the footprint of new physics in the Higgs sector.In contrast to the hadron colliders that are plagued with enormous background events, lepton colliders appear to be much more appealing options for conducting precision measurements on Higgs properties.Three promising next-generation e + e − colliders, International Linear Collider (ILC) in Japan [4,5], Future Circular Collider (FCC-ee) at CERN [6] (formerly called TLEP), and Circular Electron-Positron Collider (CEPC) in China [7,8], have been proposed in recent years.
All these three e + e − colliders plan to operate at center-of-mass (CM) energy around 250 GeV, where the dominant Higgs production channel is via the so-called Higgsstrahlung process, e + e − → HZ.For this reason, these e + e − machines are collectively referred to as Higgs factory.This golden process has been intensively studied theoretically in the past decades, e.g., the leading order (LO) prediction was first given in Refs.[9,10,11], while the next-to-leading order (NLO) electroweak corrections were analyzed in Refs.[12,13,14].Very recently, the mixed electroweak-QCD next-to-next-to-leading order (NNLO) corrections have also been investigated by two groups [15,16], which yield the state-of-the-art prediction for σ(HZ) about 230 fb around √ s = 240 GeV.The motif of this work is to study another type of Higgs production process at future lepton colliders, the associated production of Higgs boson with a hard photon, that is, e + e − → Hγ.Owing to the exceedingly small electron Yukawa coupling and the absence of Hγγ, HγZ couplings at tree level, this process first arises at one-loop order in SM.The LO prediction to σ(Hγ) in SM is available long ago [17,18,19], which turns out to be several orders-of-magnitude smaller than that of σ(HZ) around the Higgs factory energy.
Due to the highly suppressed production rate predicted in SM, the discovery prospect of the e + e − → Hγ process appears to be rather obscure in the future e + e − colliders.On the other hand, this may turn into a virtue, since this rare Higgs production process can serve as a sensitive probe for new physics search.Had the couplings between the Higgs boson and the gauge bosons/top quark been notably modified by some beyond-SM models, or had some hypothesized heavy charged particles been strongly coupled to the Higgs boson, the production cross section for σ(Hγ) might be substantially enhanced relative to its SM value, so that the e + e − → Hγ process could even possibly be observed at Higgs factory.The impact of possible new physics scenarios on this process has been extensively investigated in literatures [20,21,22,23,24,25,26].
Needless to say, an accurate SM account for the e + e − → Hγ process is crucial and mandatory for confidently interpreting the potential experimental signal in future.The goal of this work is to conduct a detailed study on the NLO QCD corrections to this process in various energy range.In a sense, the required two-loop calculation is similar to the previous works about NLO QCD corrections to H → γγ [27,28,29,30,31,32,33,34,35,36] and H → Zγ [37,38,39].Nevertheless, the situation in our case is more involved than these Higgs decay processes, due to the occurrence of a new energy scale, √ s.We are also interested in inferring how the effects of NLO QCD corrections vary with √ s.The rest of this paper is organized as follows.In Section 2, we establish the notations and briefly describe the procedure of our NLO calculation.Section 3 is dedicated to presenting our numerical results.Finally we present a summary and outlook in Section 4. For the convenience of the readers, the compact expressions for the LO amplitude are collected in Appendix.

Descriptions of the NLO calculation
Owing to the Lorentz covariance and electromagnetic current conservation, the amplitude for e + e − → Hγ can be decomposed into the linear combination of four distinct Lorentz structures [19]: where C ± 1,2 are scalar coefficients, and where p + and p − signify the momenta of the incoming positron and electron, p γ and p H signify the momenta of the outgoing photon and Higgs boson, and ε γ denotes the photon polarization vector.For simplicity, we shall neglect the tiny electron mass throughout this work.Substituting (2) into (1), squaring, averaging upon the e ± spins and summing over photon helicity, we can express the unpolarized differential cross section for e + e − → Hγ as where θ denotes the polar angle between the the outgoing Higgs and the incoming This process first occurs at one-loop order in SM, with some typical diagrams shown in Figure 1.The calculation of the LO amplitude has been comprehensively described in Refs.[17,18,19].
To assess the impact of the NLO QCD corrections, it is convenient to expand the scalar form factors C ± 1,2 in powers of the strong coupling constant: where C ±(0) 1,2 designate the LO contributions, and C ±(1) 1,2 encode the relative order-α s corrections.Substituting (4) back into (3), employing , one then readily identifies the NLO QCD corrections to the differential cross section.The LO expressions for C ±(0) 1,2 are well recorded in literatures [17,18,19].The central challenge of this work is then to compute the four form factors C ± 1,2 through order α s .
In our calculation, we employ the Feynman gauge in both electroweak and QCD sectors.Moreover, we adopt the dimensional regularization to regularize both UV and IR divergences.The LO and NLO Feynman diagrams and corresponding amplitudes are generated by FeynArts [40].We use the Mathematica packages FeynCalc [41,42]/FeynCalcFormLink [43] to carry out the trace over Dirac/color matrices, and utilize the packages Apart [44] (and another private code) and the C++ package FIRE [45] to carry out the partial fraction together with the integrationby-parts (IBP) reduction for tensor integrals, finally end up with a number of master integrals (MIs).
Similar to the loop-induced processes H → γγ, Zγ, the LO amplitude for e + e − → Hγ is also UV finite.The corresponding C ±(0) 1,2 can then be expressed in terms of the linear combinations of a number of one-loop Passarino-Veltman functions.These form factors have been presented in this format in previous work [19], while retaining a non-vanishing electron mass to regulate the potential collinear divergence (which can arise from the last diagram in Figure 1).Of course, the ultimate amplitude is free from collinear singularity and one is allowed to take the smooth limit m e → 0 in the end.As stated before, we have set m e = 0 at the outset and used dimensional regularization to regularize both UV and collinear divergences.We feel that our procedure looks simpler both conceptually and technically.For the sake of completeness, and for readers' convenience, we present our condensed expressions for 1,2 in the Appendix.The involved Passarino-Veltman functions can be worked out analytically, or can be accurately evaluated by the packages LoopTools [46] and Collier [47] in a numerical manner.
Typical Feynman diagrams for the NLO QCD corrections to e + e − → Hγ.The cap signifies the insertion of the top quark mass counterterm δm t , as given in (7).
We wish to describe more details about the calculation of the NLO QCD corrections.At NLO, all the two-loop diagrams have simple s-channel topology, some of which are illustrated in Figure 2.They can simply be obtained by dressing the gluon to the top quark loop in all possible ways.For simplicity, we have suppressed the contributions from five lighter quarks, due to their much smaller Yukawa couplings.Thanks to the simple s-channel topology, the amplitude assumes a factorized form: where P ± = 1±γ 5 2 are the chirality projection operators.
couplings for the left-handed and right-handed fermion f , respectively, where Q f and I 3 f designate the corresponding electric charge and the third component of the weak isospin of the fermion f .c W (s W ) stands for the cosine (sine) of the Weinberg angle.The two terms in (5) stem from two distinct channels, e + e − → Z * → Hγ and e + e − → γ * → Hγ, respectively.T µν VγH represent the corresponding tensor form factors related to the VγH (V = γ, Z) vertex.As mentioned before, at NLO, we only need consider the contribution to the VγH vertex from the top quark loop.

Lorentz covariance enables us to parameterize the vertex form factors T µν
VγH (V = γ, Z) as where k = p γ + p H is the 4-momentum of the virtual gauge boson V, and T V,i (i = 1, . . ., 6) are various scalar form factors as the functions of s and m H . Furry's theorem, equivalently the charge conjugation symmetry, forbids the axial vector coupling of Ztt to yield a nonvanishing contribution (at least at this order), consequently T V,6 = 0. Electromagnetic current conservation implies that T V,2 = 0, and T V,3 = −T V,5 .Furthermore, after substituting (6) into (5), one finds that T V,1 and T V,4 do not contribute to the amplitude.Consequently, we just need to calculate a single scalar form factor T V,5 .We emphasize that, the mere QCD renormalization required in this work is to implement the top quark mass renormalization, with the mass counterterm: where the spacetime dimensions d = 4 − 2ǫ, and m t signifies the top quark pole mass.The counterterm diagrams are represented by the last two diagrams in Figure 2.
Having eliminated the UV divergences in T V,5 (V = Z, γ) that enter the NLO amplitude in (5), after some straightforward manipulations, one can identify the NLO coefficients C ±(1) 1,2 with where κ ≡ T Z,5 /T γ,5 = − It is easy to understand why T Z,5 and T γ,5 are proportional to each other, since only the vector part of the Ztt coupling survives in the amplitude.Plugging ( 8) into ( 4) and (3), we then deduce the desired NLO QCD corrections to the differential cross section.
We pause here to briefly describe our numerical strategy of computing the NLO corrections.After IBP reduction, we end up with about 34 two-loop MIs for T γ,5 .For each MI, we combine FIESTA [48]/CubPack [49] to perform the sector decomposition and the subsequent numerical integrations with quadruple precision.We also utilize the Fortran package Cuba [50] to crosscheck the numerical integrations for most MIs.
It is worth mentioning that our two-loop calculation is quite similar to the corresponding calculation for the rare decay process H → Zγ.The NLO QCD corrections to that process were recently accomplished by two groups in an analytic fashion [38,39].As a crosscheck, we have also revisited that process and found good agreement with the numerical results tabulated in [38].Having passed such a nontrivial test, we feel confident about the correctness of our NLO QCD calculation for e + e − → Hγ.
In principle, our T Z,5 can be obtained from its counterpart in H → Zγ by performing some sort of analytic continuation.Unfortunately, the analytic expressions recorded in Refs.[38,39] involve rather lengthy and complicated special functions, rendering analytic continuation a rather nontrivial task.Therefore, we are contented with only presenting the highly accurate numerical predictions in this Letter.

Numerical results
In order to make concrete predictions, we specify the following input parameters, in accordance with the latest compilation of the Particle Data Group [51]: In addition, the default renormalization scale µ = √ s is assumed for the QCD running coupling constant α s (µ).We take the initial value α (5)  s (M Z ) = 0.118, and use the two-loop running formula to evolve to another different energy scale.Practically, it is very convenient to utilize the package RunDec [52], which has automated the evolution of the QCD running coupling and taken care of the top quark threshold effect.
Before proceeding into phenomenological discussions, we should remind the readers that our LO predictions are subject to an intrinsic ambiguity due to the choice of the electroweak coupling constants.In this work, we have chosen to use the Thomson-limit value of QED coupling α(0).Certainly, it is equally legitimate to adopt different values for α, such as the widely-used α G µ and α(M Z ) [53].This uncertainty might be even more pronounced than the size of the NLO QCD corrections.This ambiguity can only be lessened once one incorporates the full NLO weak corrections, unfortunately which is currently unavailable.
In Figure 3 and 4, we plot the angular distributions of the Higgs boson at two benchmark CM energy points, √ s = 240 GeV at CEPC, and 500 GeV at ILC.As the CM energy increases, the angular distributions turn to be more and more concentrated in the forward/backward directions.As one can readily tell, the effect of NLO QCD corrections at CEPC energy, which only constitutes ∼ 0.2% of the LO contribution, is hardly discernible.In the meanwhile, the impact of the NLO QCD corrections becomes clearly visible at ILC energy range.
To ascertain the integrated cross section at each CM energy point, we utilize Cuba [50] to carry out the angular integration numerically.In Figure 5, the variation of the LO cross section with the CM energy is depicted over a wide range.There the line shape exhibits some interesting features.Starting from the threshold √ s = m H , the cross section rises from zero due to the increment of the phase space.The cross section reaches its first peak around √ s = 245 GeV, followed by a dip around √ s = 359 GeV, then slowly bounces up and reaches the second lower peak centered at √ s = 480 GeV.The cross section then decreases monotonically with increasing √ s.To understand the origin of the peak-dip-peak structure, we separate the contributions into two categories, those from the diagrams involving the top quark loop and from the diagrams involving electroweak gauge bosons.There arises a cusp from the former channel, exactly located at the tt threshold √ s = 2m t = 348.4GeV.From Figure 5, one observes that there arises a substantial destructive interference between these two channels, whose effects reach the maximum at √ s ≈ 352 GeV.This destructive interference seems to be largely responsible for the global minimum of the net LO line shape.It is interesting to note that the location of this dip deviates from the tt threshold by 10 GeV.
We wish to understand the asymptotic behavior of σ(Hγ) in the high energy limit.For this purpose, it is instructive to temporarily digress into the asymptotic scaling of the total cross section of a similar Higgs production channel, the Higgsstrahlung process.It is well-known that σ(HZ) scales as 1/s in the high energy limit [12,13,14].This scaling is saturated by the channel with a longitudinally-polarized Z, while the polarized cross section for a transverselypolarized Z declines at a much steeper rate, ∝ 1/s 2 [14,15].
Since the photon must be transversely polarized, one might naively anticipate that σ(Hγ) should decrease as 1/s 2 asymptotically.However, a careful numerical investigation reveals that this expectation is not true!At sufficiently high energy, the cross section turns out to decline with a much slower pace, ∝ 1/s.In addition to this power-law scaling, we have not observed any hint of logarithmic enhancement, to a decent numerical accuracy.Therefore, the To trace the origin of the nontrivial line shape, we deliberately isolate the contributions from two classes of diagrams.The dotted, dashed and dot-dashed lines represent the contribution from diagrams involving the top quark loop, that from all other diagrams involving weak gauge bosons in the loop, and their interference, respectively.The integrated cross sections of e + e − → Hγ in a variety range of CM energy, at both LO and NLO accuracy.We also enumerate the corresponding values of the scalar form factors T γ,5 in (8).For simplicity, we have eliminated the α-dependence by redefining T γ,5 ≡ α −3/2 T γ,5 .

√ s(GeV
ln n s (n = 1, 2) terms, if exist, must be accompanied with the power-suppressed terms (∝ 1/s 2 ).It is interesting to trace the origin for the 1/s scaling of σ(Hγ) at high energy.Note there is one important difference between Higgsstrahlung and the process considered in this Letter, e.g., the former starts at tree level, but the latter starts at one loop order.The s-channel diagrams in Figure 1 indeed only yield a suppressed contribution ∝ 1/s 2 , modulo possible logarithms of s.Nevertheless, the box diagrams in Figure 1, which are responsible for the sharp concentration of the differential cross sections in the forward/backward zones, will bring forth an enhancement factor ∝ s M 2 (here M refers to a generic mass scale of order M Z,W or M H ) upon angular integration.A closer analytic examination also supports the absence of logarithms of s accompanying this leading-power scaling.
In Table 1 we enumerate the predicted cross sections over a wide range of CM energy, at both LO and NLO accuracy.In addition, we also tabulate the values of the scalar form factors T γ,5 at various CM energies.At CEPC energy range, √ s ≈ 250 GeV, the NLO QCD corrections have a negligible effect.In contrast, if the ILC and FCC-ee operate at √ s ≈ 400 GeV, the magnitude of the NLO QCD corrections could reach roughly 20% of the LO cross section.
As can be noticed in Table 1, as √ s > 2m t , the two-loop form factor T γ,5 develops an imaginary part due to the opening of the tt threshold.Through an exhaustive threshold scan, we observe that the real part of T γ,5 develops a logarithmic singularity ∝ ln |β| (β ≡ 1 − 4m 2 t s ) in the vicinity of the threshold, while the imaginary part approaches a constant above from the tt threshold.This is a clear sign that the nearly-on-shell tt pair in the loop carries the dominant quantum number of 3 S 1 .It is well-known in such a case, the Coulomb singularity ∝ 1 β will arise for even higher-order diagrams containing more Coulomb gluon exchange, so that fixed-order QCD perturbative expansion inevitably breaks down near the threshold.For inclusive tt production near the threshold, a vast amount of literatures have been devoted to providing a systematic description in the threshold region, by resumming Coulomb singularities to all orders and incorporating the finite t width effect, in the context of non-relativistic effective field theory NRQCD [54,55,56,57,58].For our exclusive H + γ process, one may closely follow the ansatz adopted for the pseudosalar Higgs decay into two photons when the pseudoscalar Higgs mass is about twice the top quark mass, where a tt pair is in the 1 S 0 state in the loop [59,60].The resummation of the Coulomb gluon ladder diagrams, as well as including top quark width, could also be readily fulfilled, so that one can obtain reliable predictions for σ(Hγ) in the vicinity of √ s = 2m t .Such a "nonperturbative" treatment is beyond the scope of this Letter, and we plan to present a comprehensive threshold analysis for this process in the near future.In any rate, our fixed-order predictions recorded in Table 1, should still be viewed as trustworthy, as long as √ s does not lie in proximity to the tt threshold.For the sake of clarity, in Figure 6 we also plot the integrated cross section as a function of the CM energy, explicitly including the NLO QCD corrections.In most energy ranges, the impact of the NLO QCD corrections seems rather modest, especially when √ s < 290 GeV, with a relative effect typically less than 1%.Nevertheless, with the increasing CM energy, the QCD corrections appear to become more pronounced, and reach the maximum about 20% at √ s ≈ 400 GeV.The effect of NLO corrections diminishes again as the CM energy further increases, completely negligible after √ s > 1 TeV.This can be readily understood in light of the preceding discussions: due to the dominance of the electroweak box diagrams, the LO cross section exhibits a mild 1/s asymptotic scaling, while the NLO QCD corrections fall with a much faster pace ∝ 1/s  Finally in Table 2, we further study the dependence of the NLO QCD corrections on the Higgs boson mass, with √ s held at 240 GeV.The LO cross section decreases monotonously with the rising m H , partly due to the shrinking phase space.The magnitude of the NLO QCD corrections varies with m H with a modest pace, nevertheless always less significant than 1%.

Summary and Outlook
The LO predictions for the production rate of e + e − → Hγ in SM are rather small at the next-generation e + e − colliders, which may render it a sensitive probe for the elusive beyond-SM signals.An accurate account of this process within SM appears mandatory, if one wishes to confidently confront the potential new physics signals in the future.In this Letter, we have, for the first time, calculated the NLO QCD corrections to this associated H + γ production process, and conducted a comprehensive numerical study covering a wide range of CM energy.The key finding is that, the NLO QCD corrections to this process at CEPC can be safely neglected, but may have sizable impact at ILC energy range.Nevertheless, the e + e − → Hγ process appears to possess the maximal production rate, σ max ≈ 0.08 fb, around √ s = 250 GeV, therefore CEPC appears to be the ideal place to look for this rare Higgs production process.
It may also be valuable to reanalyze this process for polarized electron beam, as is useful for the ILC experiment.
As was mentioned before, in order to obtain a precise prediction for the e + e − → Hγ process, it seems mandatory to lessen the large ambiguity inherent in the choice of the electroweak coupling constants.This symptom can only be relieved by further including the NLO weak corrections.As was witnessed in the Higgsstrahlung process, incorporating the NLO weak corrections does significantly stabilize the predicted production rate for HZ [14,15].This is particularly important for our process, since the ambiguity due to the choice of α might even bring forth greater uncertainty than the magnitude of NLO QCD corrections per se.Despite its daunting difficulty, we hope that future fulfillment of the NLO weak corrections to this process will significantly reduce the α-scheme dependence.Combined with the NLO QCD corrections calculated in this Letter, we will then be able to provide a fairly reliable guidance for the future experimental search of this rare Higgs production channel.
Finally, we think it is also of some interest to conduct a full-fledged study for σ(Hγ) in the vicinity of the tt threshold, which necessitates to include resummation of the Coulomb gluon ladder diagrams and top quark width effect.
Appendix: The compact expressions for C ±(0) 1,2 For reader's convenience, here we list the compact expressions for the LO coefficients functions C ±(0) 1,2 appearing in (3), with m e set to zero at the outset.They can be split into three parts: where i = 1, 2, and where We use the same conventions for the Passarino-Veltman scalar functions as in most modern literatures (see, for instance, Ref. [47]).Note that the functions (12e), which originate from the last diagram in Figure 1, are by themselves free from any collinear singularity.However, if the IBP reduction is further implemented to (12e), one will end up with a linear combination of a number of one-loop scalar MIs, some of which are IR divergent.Fortunately, one can use LoopTools [46] and Collier [47] to directly calculate these Passarino-Veltman functions in (12).

Figure 3 :Figure 4 :
Figure 3: Angular distributions of the Higgs boson in the e + e − → Hγ process at √ s = 240 GeV.The right panel embodies the relative magnitude of the NLO QCD corrections.

Figure 5 :
Figure 5: The LO cross section as a function of √ s (the solid line).To trace the origin of the nontrivial line shape, we deliberately isolate the contributions from two classes of diagrams.The dotted, dashed and dot-dashed lines represent the contribution from diagrams involving the top quark loop, that from all other diagrams involving weak gauge bosons in the loop, and their interference, respectively.

Figure 6 :
Figure 6: The total cross section as a function of √ s, both at LO and NLO in α s .The vertical band with √ s = 2m t ± 5 GeV signifies the threshold region, inside which the perturbative expansion is expected to break down and our fixed-order predictions become invalid.

Table 1 :
2, since the corresponding diagrams possess only s-channel topology.

Table 2 :
Cross sections of e + e − → Hγ as a function of Higgs boson mass, with √ s fixed at 240 GeV.