Mixed electroweak-QCD corrections to $e^+e^-\to \mu^+\mu^- H$ at CEPC with finite-width effect

The associated production of Higgs boson with a muon pair, $e^+e^-\to \mu^+\mu^- H$, is one of the golden channels to pin down the properties of the Higgs boson in the prospective Higgs factories exemplified by \textsf{CEPC}. The projected accuracy of the corresponding cross section measurement is about per cent level at CEPC. In this work, we investigate both ${\mathcal O}(\alpha)$ weak correction and the ${\mathcal O}(\alpha\alpha_s)$ mixed electroweak-QCD corrections for this channel, appropriately taking into account the effect of finite $Z^0$ width. The $\mu^+\mu^-$ invariant mass spectrum is also predicted. The mixed electroweak-QCD correction turns out to reach 1.5\% of the Born-order result, and thereby must be included in future confrontation with the data. We also observe that, after including higher-order corrections, the simplified prediction for the integrated cross section employing the narrow-width-approximation may deviate from our full result by a few per cents.


I. INTRODUCTION
Ever since the ground-breaking discovery of the Higgs boson at Large Hadron Collider (LHC) in 2012 [1,2], one of the highest priorities of particle physics is to nail down the properties of the Higgs boson as precise as possible. Unlike the hadron colliders, which suffer from severe contamination due to the copious background events, the electron-positron colliders provide an ideal platform to precisely measure various Higgs couplings [3]. In recent years, three next-generation e + e − colliders have been proposed for dedicated study of Higgs boson: International Linear Collider (ILC) [4], Future Circular Collider (FCC-ee) [5], and Circular Electron-Positron Collider (CEPC) [6], all of which plan to operate at the centerof-mass energy around 240 ∼ 250 GeV.
There emerge several Higgs production mechanisms at e + e − colliders: Higgsstrahlung, W W fusion and ZZ fusion, etc.. Around √ s ≈ 240 GeV, which is the projected energy range of CEPC, the Higgs production is dominated by the Higgsstrahlung channel e + e − → ZH, Higgs production associated with a Z 0 boson. It is anticipated that, with the aid of very high luminosity and the recoil mass technique, CEPC can measure the Higgs production cross section with an exquisite sub-per-cent accuracy. Needless to say, it is indispensable for theoretical predictions for the Higgsstrahlung channel to be commensurate with the projected experimental precision. The leading order (LO) prediction for e + e − → ZH was first considered in 70s [7][8][9][10]. In the early 90s, the next-to-leading order (NLO) electroweak correction for this process has also been addressed by three groups independently [11][12][13], which turns out to be significant. Very recently, the mixed electroweak-QCD next-to-next-to-leading (NNLO) corrections were also be independently calculated by two groups [14,15]. The O(αα s ) correction may reach 1% of the LO prediction, thereby must be included when confronting the future measurement. Recently, the ISR effect of this process has also been carefully analyzed [16].
From the experimental angle, it is the decay products of the Z 0 boson, rather than the Z 0 itself that are tagged by detectors in the Higgsstrahlung channel, since the Z 0 is an unstable particle. Therefore, in order to get closer contact with experiment, it is advantageous to make precise predictions directly for the process e + e − → (Z * →)ff +H, where f represents leptons or quarks. Among a flurry of Higgs production channels associated with various Z decay products, the e + e − → µ + µ − H process occupies a unique place for probing Higgs properties, because it is a very clean channel and possesses large cross section. The production cross section for this individual channel can be measured with 0.9% precision at CEPC [6,17]. Combining several other channels, CEPC is anticipated to measure the Higgs production rate with the accuracy of 0.51%.
The LO contribution to the e + e − → ff + H process was first considered in 70s [18]. The initial-state-radiation (ISR) correction to these types of processes was addressed in 80s [19]. There exist a flurry of higher-order studies for the process e + e − → ννH, where both Higgsstrahlung and W W fusion mechanisms contribute [20][21][22][23][24][25][26][27]. To our knowledge, there appears no dedicated work to investigate the NLO weak correction to e + e − → µ + µ − H. Nevertheless, the NLO weak correction to a similar process e + e − → e + e − H were calculated by the GRACE group more than a decade ago [28,29]. One can extract the corresponding NLO weak correction to e + e − → µ + µ − H by singling out a subset of diagrams in [28,29].
The purpose of this work is to conduct a systematic investigation on the higher-order radiative corrections to the process e + e − → µ + µ − H, to match the projected experimental precision at CEPC. We first compute the NLO weak correction to e + e − → µ + µ − H, then proceed to include the O(αα s ) mixed electroweak-QCD correction. Besides the integrated cross section, we also study the impact of radiative corrections to various kinematic distributions such as the µ + µ − invariant mass distribution. For this purpose, the finite Z 0 width effect must be consistently taken into account. It is also instructive to examine how our results deviate from those obtained by invoking the narrow width approximation (NWA).
The rest of the paper is structured as follows. In section II, adopting the Breit-Wigner ansatz for the resonant Z 0 propagator, we recapitulate the LO prediction for e + e − → µ + µ − H and also show the corresponding NWA result. In section III, we specify our strategy of implementing the finite Z 0 -width effect in higher-order calculation. In section IV, we present the calculation for the NLO weak correction to this channel. In section V, we describe the calculation for the mixed electroweak-QCD corrections. In section VI, we present the numerical results and phenomenological analysis. Finally we summarize in section VII. We are considering the process

II. LEADING ORDER RESULTS AND NARROW WIDTH APPROXIMATION
where the momenta of the incoming and outgoing particles are specified in the parentheses. For future usage, we define s ≡ (k 1 + k 2 ) 2 , and s 12 ≡ (p 1 + p 2 ) 2 . For convenience, we also define the invariant mass of the muon pair by M µµ ≡ √ s 12 , which lies in the range At Higgs factory, lepton masses can be safely neglected owing to their exceedingly small Yukawa couplings. Consequently at the lowest order, there is only a single s-channel diagram as depicted in Fig. 1. The LO amplitude reads where c W ≡ cos θ W , s W ≡ sin θ W , with θ W the Weinberg angle, M Z represents the mass of the Z 0 boson.
is the coupling of the gauge boson and the charged lepton. Specifically speaking, The chirality structure of the neutral current demands that, e + and e − (also µ + and µ − ) must carry opposite helicity in order to render a non-vanishing amplitude.
As can be readily seen from Fig. 1, it is possible for the µ + µ − pair to be resonantly produced from the on-shell Z 0 boson, consequently the amplitude in (2) blows up at s 12 = M 2 Z , which reflects that fixed-order calculation breaks down near the Z pole. To tame the singularity in the limit s 12 → M 2 Z , it is customary to replace the second Z boson propagator in (2) with the Breit-Wigner form, which amounts to include the Dyson summation for the Z boson self-energy diagrams. Retaining finite Z width would effectively cutoff the IR singularity. One may define a new amplitude: where F is a rescaling factor, and Γ Z signifies the width of the Z 0 boson. The LO cross section is then given by where the three-body phase space in the center-of-mass (CM) frame can be conveniently parameterized as where (|p * 1 |, Ω * 1 ) signifies the 3-momentum of the µ − in the rest frame of the dimuon system, |p H |, θ H represent the magnitude of the momentum and the polar angle of the Higgs boson in the laboratory frame, respectively. Upon neglecting masses of the electron and muon, one obtains |p * In deriving (5), we have utilized the axial symmetry to eliminate the trivial dependence on the azimuthal angle of the outgoing Higgs boson.
Squaring (2), summing over µ + µ − helicities, and averaging upon the e + e − polarizations, one observes that the squared amplitude bears a factorized structure, thanks to the simple s-channel topology. Substituting it into (4), integrating over the solid angle Ω * 1 , one then arrives at the following double differential cross section: with α ≡ e 2 4π the electromagnetic fine structure constant. Integrating (6) over the polar angle, one obtains the Born-order spectrum of the invariant mass of µ + µ − : Since Γ Z ≪ M Z , one naturally expects that the NWA should be fairly reliable for the process under consideration. Inserting the limiting formula into (6), and integrating over s 12 , we obtain the angular distribution: where is the angular distribution of the Z(H) in the process e + e − → ZH at Born order, with In (9), the Born-order partial width and branching fraction of Z → µ + µ − are given by From (9), one readily obtains the LO integrated cross section in the NWA ansatz: where Note that the unpolarized LO cross section σ 0 (ZH) in (13) decreases rather mildly (∝ 1/s) in the high energy limit, reflecting the dominance of producing the longitudinally polarized Z in large √ s. However, at moderate energy such as √ s = 250 GeV at CEPC, the longitudinallypolarized cross section only comprises of 42% of the total unpolarized cross section.

III. THE TREATMENT OF FINITE Z 0 WIDTH IN HIGHER-ORDER CORREC-TIONS
As mentioned before, in this work we are interested in addressing the NLO weak and mixed electroweak-QCD corrections for e + e − → µ + µ − H: For simplicity, in this work we have neglected the pure QED corrections (such as ISR and FSR effect), which can instead be simulated by the package Whizard [30]. As a consequence, a simplifying feature arises that the dominant higher-order diagrams resemble the s-channel topology as depicted in Fig. 1, which contains only one resonant Z propagator.
Once going beyond LO, it becomes a quite delicate issue to incorporate the finite Z width effect yet without spoiling gauge invariance and bringing double counting. Over the past decades, numerous practical schemes have been proposed to tackle the unstable particle, such as the pole scheme [31][32][33], factorization scheme [34,35], fermion-loop scheme [36,37], boson-loop scheme [38], complex mass scheme [39,40], etc.. It is worth mentioning that a systematic and model-independent approach, the unstable particle effective theory, has also emerged finally [41,42]. However, this approach is valid only near the resonance peak, and cannot be applied in the entire kinematic range.
Owing to the particularly simple s-channel topology of our process, it is most convenient to employ the factorization scheme [34,35], which is particularly suitable for such resonancedominated process. In this scheme, one rescales a gauge-invariant higher-order amplitude by a Breit-Wigner factor F , and subtracting the iM Z Γ Z terms which potentially generates double counting. The merit of this scheme is that gauge invariance is preserved, and can be readily implemented in automated calculation. Recently this scheme has also been used by Denner et al. to analyze the NLO electroweak correction to e + e − → ννH [25].
For our purpose, we specify the recipe of the factorization scheme closely following [25]: where n = 0, 1, M represents the fixed-order amplitude where the Z 0 is treated as a rigorously stable particle, F and M 0 have been defined in (3),Σ T ZZ (s) represents the transverse part of the renormalized one-particle irreducible self-energy diagrams for Z boson. M Z is the pole mass of the Z 0 boson, and throughout the work we take Γ Z as the experimentally determined Z 0 boson width 1 .
The second term in the right-hand side of (15) is included to subtract the double-counting term. Fortunately, due to its orthogonal phase, the interference of this term with M generates a purely imaginary contribution to the cross section, thus can be safely neglected.
Once the rescaled O(α) and O(αα s ) amplitudes are obtained, we then deduce the corresponding higher-order corrections to the differential cross section through with n = 0, 1. Note even for the mixed electroweak-QCD correction, we only need consider its interference with the Born-order amplitude. We conclude this section by stressing that, since the non-resonant diagrams are regular at s 12 = M 2 Z , the rescaling procedure in (15) enforces their contributions to the amplitude to vanish on the Z 0 pole. In the vicinity of the resonance, it is intuitively appealing that the non-resonant diagrams are much more suppressed relative to the resonant diagrams. As will be seen in section VI, our numerical predictions indeed confirm this anticipation.

IV. CALCULATION OF THE NLO WEAK CORRECTION
We now outline the calculation of the NLO weak correction to e + e − → µ + µ − H, with some representative diagrams depicted in Fig. 2 and 3. As stressed before, we will not 1 By default, the pole mass of the Z 0 is determined by the condition Re Σ ZZ (M 2 Z ) ≡ 0, whereas its width is inferred from the optical theorem, M Z Γ Z = Im Σ ZZ (M 2 Z ) . It was argued [32,33,43] that the pole mass and the corresponding width defined this way are gauge dependent. Nevertheless, the gaugedependent terms arise at order-α 3 in the 't Hooft-Feynman gauge, which is beyond the accuracy targeted in this work. Thus we will pretend M Z and Γ Z to be gauge-invariant quantities.  Fig. 3. Diagrams in the first two rows correspond to the "resonant" channel e + e − → (Z * /γ * →)µ + µ − + H, while those in the last row exhibit a completely different "non-resonant" topology. consider the ISR and FSR types of diagrams. It is obvious that the NLO diagrams can be separated into two gauge-invariant subgroups, with either "resonant" or "non-resonant" structures. For the former subset, the diagrams are very similar to those encountered in the previous NLO weak correction for e + e − → ZH, so are the corresponding calculations; for the latter, there emerges no singularity as s 12 → M 2 Z , so there is no need to include width effect for any particle routing around the loop.
The NLO amplitude is computed in Feynman gauge. Masses of all light fermions are neglected except the top quark. Dimensional regularization (DR) is employed to regularize UV divergence. The Feynman diagrams and the corresponding amplitude are generated by the package FeynArts [44]. Tensor contraction and Dirac/color matrices trace are conducted by using FeynCalc and FeynCalcFormLink [45][46][47]. Tensor integrals are further reduced to the Passarino-Veltman scalar functions, which are numerically evaluated by Collier [48] and LoopTools [49].
We also choose to use the standard on-shell renormalization scheme to sweep UV divergences, where various electroweak counterterms are tabulated in [50]. Depending on the specific recipe for the charge renormalization constant Z e , there are three popular subschemes of the on-shell renormalization: α(0), α(M Z ) and G µ schemes [13]. In the first scheme, the fine structure constant α is assuming its Thomson-limit value, whereas α(0) is replaced with in the α(M Z ) and G µ schemes, respectively. Differing from the α(0) scheme, these two schemes effectively resum either some universal large logarithms from the light fermion loop or some m 2 t -enhanced terms from the top quark loop. Once the M (α) is rendered finite after the renormalization procedure, we then employ (15) to obtain the rescaled amplitude M (α) , which encapsulates the finite Z-width effect. It is then straightforward to utilize (16) to infer the NLO weak correction to the differential cross section.

V. CALCULATION OF MIXED ELECTROWEAK-QCD CORRECTIONS
Finally we turn to the O(αα s ) mixed electroweak-QCD correction to e + e − → µ + µ − H. Since it is the quarks instead of leptons that can experience the strong color force, we only need retain those diagrams involving quark loop. Moreover, since the top quark couples the Higgs boson with the strongest strength, for simplicity we have neglected the masses of all lighter quarks, so we only retain those two-loop diagrams where only the top quark loop dressed by gluon. Some typical two-loop diagrams are shown in Fig. 2 and 3, bearing only the s-channel "resonant" structure. As indicated in Fig. 3, at this order, QCD renormalization is realized by merely inserting the one-loop top quark mass counterterm, δm t , into the internal top-quark propagator, as well as into the Htt vertex [15]. The calculation very much resembles our preceding work on O(αα s ) correction to e + e − → ZH [15], and we referred the interested readers to that paper for more details.
For the actual two-loop computation, we utilize the packages Apart [51] and FIRE [52] to perform partial fraction and integration-by-parts (IBP) reduction. We then combine FIESTA [53]/CubPack [54] to perform sector decomposition and subsequent numerical integrations for master integrals with quadruple precision.
Besides the finite renormalization of Zee vertex [15], the O(αα s ) amplitude can be ex-pressed in terms of the Born-order amplitude supplemented with an effective HV V vertex: where the sum is extended over V 1 , V 2 = Z 0 , γ, and −ieT µν HV 1 V 2 is the HV 1 V 2 effective vertex, which depends on K = k 1 + k 2 and P = p 1 + p 2 . The gauge boson V 1 is coupled with the incoming e + e − pair, whereas the gauge boson V 2 is affiliated with the outgoing µ + µ − pair. Γ µ V represents the coupling between the gauge boson and charged leptons, whose form has already been specified in the paragraph after (2). The electromagnetic coupling of lepton is chiral symmetric, g ± γ = 1. By Lorentz covariance, the renormalized vertex tensor T µν HV 1 V 2 can be decomposed as where T i (i = 1, · · · , 6) are Lorentz scalar solely depending on s, s 12 and M 2 V 1,2 . Furry theorem enforces that T 6 = 0, technically because C-invariance forbids a single γ 5 to emerge in the trace over the fermionic loop. Owing to the current conservation associated with massless leptons, it turns out that only the scalar form factors T 4,5 survive in the differential cross sections.
Substituting (19) into (18), utilizing the factorization scheme (15) to implement the finite Z 0 width effect, we then obtain the rescaled amplitude M (ααs) . From (16), we find the O(αα s ) mixed electroweak-QCD correction to the differential cross section to be where
We adopt the following values for the input parameters [22]: and fix µ = √ s/2 for the QCD coupling, and take α s ( √ s/2) = 0.1135. The impact of the O(α) and O(αα s ) corrections to the Higgs angular distribution is quite analogous to what is found in our preceding work on e + e − → ZH [15]. In Fig. 5, we also show the µ + µ − invariant mass spectrum at √ s = 240 GeV, including both NLO weak and mixed electroweak-QCD corrections. As expected, the spectrum develops a sharp Breit-Wigner peak around the Z resonance. The O(α) and O(αα s ) corrections play a very minor role except in the proximity of the Z pole. It is interesting for the future measurement of the di-muon spectrum at CEPC to examine our predictions.
In Table I, we supplement more details for the dimuon invariant mass spectrum. We divide the NLO weak correction into the contribution from the resonant diagrams and the one from non-resonant diagrams. As can be seen from the Table I, the O(α) correction is saturated by the resonant diagrams almost in the entire energy range, especially near the Z peak.
Our goal is to present to date the most comprehensive predictions for the e + e − → µ + µ − H process, taking into various sorts of theoretical uncertainties account. In Table II, we present  The results are provided with three renormalization sub-schemes. We also include the un-   certainty inherent in the input parameters (first error) and the uncertainty due to the QCD renormalization scale (second error). To assess the parametric uncertainty, we vary the values of M W and m t , and ∆α (5) had around the central PDG values within the 1σ bands. For the QCD scale uncertainty, we slide the µ in α s from M Z to √ s.
From Table II, we observe a very similar pattern of scheme and parametric dependence of higher-order corrections as [15]. While the parametric and scale uncertainties of the NNLO predictions in the α(0) and α(M Z ) schemes are both about 0.5% of the NNLO results, the relative errors are somewhat reduced in the G µ scheme (≈ 0.2%). We also find that in the G µ scheme, the mixed electroweak-QCD corrections only amount to 0.4% of LO cross section, which might be attributed to the fact that in addition to the running of α, universal corrections to the ρ parameter are also absorbed into the LO cross section. As can also be seen in Table II, though the predicted LO cross sections from three renormalization schemes differ significantly, including the NLO weak correction significantly help them converge to each other. Including mixed electroweak-QCD correction appears not to further reduce the scheme dependence. To yield a scheme-insensitive prediction, it appears to be imperative to continue to compute the NNLO electroweak correction, which is certainly an extremely daunting task.
Since Γ Z ≪ M Z , and the production rate is predominantly saturated by the Z 0 resonance. It may seem natural to anticipate that the NWA remains valid even after including higher order corrections. Under the assumption of NWA, one may approximate the LO cross section and the higher-order radiative corrections by where σ(ZH) represents the Higgsstrahlung cross section, with σ 0 (ZH) given in (13). Br 0 is defined in (11), and the radiative corrections Br (αα n s ) (n = 0, 1) can be read off from Since the width of the Z 0 is held fixed, the perturbative expansion for the branching fraction of Z 0 → µ + µ − amounts to the expansion for the corresponding partial width.  In Table III, we compare the predicted e + e − → µ + µ − H cross section from the literal full calculation with that from NWA. For the sake of concreteness, we take √ s = 240 GeV, and employ the α(0) scheme. At LO, the NWA prediction is about 3% higher than the full prediction, while O(α) and O(αα s ) corrections are observed to be only slightly different. As a consequence, the NWA prediction to the total cross section at NNLO accuracy turns out to be about 4% higher than the full NNLO prediction.

VII. SUMMARY
Higgsstrahlung is the leading Higgs production mechanism at CEPC. The mixed electroweak-QCD correction to e + e − → ZH has recently become available [14,15]. This piece of NNLO correction appears to be surprisingly large, about 1% of the Born-order result, therefore must be considered when matching the exquisite experimental accuracy.
To make closer contact with the actual experimental measurement, in this work we have investigated both NLO weak and mixed electroweak-QCD corrections to one of the golden mode in CEPC, i.e. e + e − → µ + µ − H, with the finite Z 0 width properly accounted. At √ s ≈ 240 GeV, the NLO weak correction may reach 6% of the Born order cross section, while the NNLO mixed electroweak-QCD correction can reach 1.5% of the LO cross section, greater than the projected experimental accuracy of 0.9%. We also present numerical predictions to various differential cross sections at NNLO accuracy, in particular we predict the µ + µ − invariant-mass spectrum of the Breit-Wigner shape. We have also compared our full predictions with those based on the NWA, and found the agreement within a few percents. It is interesting to await the future experiment to examine our predictions.
We also carefully address the issue about scheme-dependence of our predictions, at various levels of perturbative accuracy. Employing three popular renormalization sub-schemes, we find that the predicted LO cross sections substantially differ from each other. Including the NLO weak correction is crucial to stabilize the predictions from different schemes, however including mixed electroweak-QCD correction seems not to help. To yield a schemeinsensitive prediction, it appears to be compulsory to continue to include the NNLO electroweak correction,