Measurements of decay branching fractions of ${{H \to b\bar b/c\bar c/gg}}$ in associated ${{(e^ + e^ - /\mu ^ +\mu ^- )H}}$ production at the CEPC

The discovery of a scalar boson with a mass around 125 GeV at the Large Hadron Collider (LHC) [1, 2] completed the final piece of the standard model (SM). This particle, interpreted as the Higgs boson, plays a crucial role in the Electroweak Spontaneous Symmetry Breaking (EWSB), known as the Higgs mechanism [3-5]. The Higgs mechanism allows the W boson and Z boson to be massive while keeping the $ SU(2)_{\rm L} \times U(1)_{\rm Y} $ gauge invariance. As a consequence of this mechanism, the fermions like quarks and charged leptons get their masses from their couplings to the Higgs field. The masses of the fermions ($ m_{f_i} $) in the SM are proportional to their Yukawa couplings ($ h_i $) to the Higgs field: $ m_{f_i} = vh_i/\sqrt{2} $, where $ v\approx 246\;{\rm{GeV}} $ is the vacuum expectation value of the Higgs field. Thus measuring the Yukawa couplings between the Higgs boson and the SM fermions is essential to understand the origin of the fermions' masses. The deviation of these couplings from SM prediction would indicate new physics.

The dominant Higgs boson decays into fermionic final states are $H \to b\bar b$, $H \to {\tau ^ + }{\tau ^ - }$ and $H \to c\bar c$, the decay branching fractions of which are predicted to be 57%, 6% and 2.7% respectively [6-17]. In addition, the Higgs boson can decay to a pair of gluons via heavy quark loops. The large coupling between the Higgs boson and the top quark leads to considerably large branching fraction of $H \to gg$ which is estimated to be about 9% [18-21].

Until now, all the Higgs boson measurements are performed in hadron colliders. The leading fermionic Higgs boson decay, $H \to b\bar b$, was studied in both ATLAS and CMS, using the LHC Run-I data, which contains about 5 fb^-1 and 20 fb^-1 of $ pp $ collision data with $ \sqrt{s} $ of 7 TeV and 8 TeV respectively. These measurements include several Higgs boson production channels: $ VH $ [22, 23], $ t\bar{t}H $ [24-26] and VBF [27, 28] processes. The $H \to b\bar b$ were also studied at Tevatron [29] in $ VH $ production, using 9.7 fb^-1 $ p\bar{p} $ collision data with $ \sqrt{s} $ of 1.96 TeV. The $ H \to b\bar b $ signal strength, defined as the ratio of the measured cross section to the corresponding SM prediction, is estimated to be $ 0.70\pm0.29 $ according to the combination of ATLAS and CMS analysis of run-I data [30]. In 2018 observations of $ H \to b\bar b $ decay in $ VH $ production were declared by ATLAS [31] and CMS [32], using 79.8 and 41.3 fb^-1 $ pp $ collision data with $ \sqrt{s} $ of 13 TeV. The signal strength is $ 1.16\pm 0.16(\rm stat.)^{+0.21}_{-0.19}(\rm sys.) $ and $ 1.01\pm0.22 $ respectively. The $ H \to c\bar c $ was also studied using 36.1 fb^-1 data with $ \sqrt{s} $ of 13 TeV in ATLAS [33], giving a upper limit about 100 times higher than the SM prediction with 95% confidence level. The precision of those results is limited by large QCD production background, which is inevitable in hadron colliders.

A lepton collider has significant advantage in precise Higgs measurements as it's free of QCD production background and a has precise and tunable initial energy. Several future lepton colliders have been proposed with the capability of precise measurement of Higgs boson parameters, such as the International Linear Collider (ILC) [34], the $ e^{+}e^{-} $ Future Circular Collider (FCC-ee) [35], the Compact Linear Collider (CLIC) [36] and the Circular Electron Positron Collider (CEPC) [37]. The CEPC is a proposed electron-positron collider by the Chinese high energy physics community. It can be operated with $ \sqrt{s} $ of 240 GeV to 250 GeV and the designed instaneous luminosity is 2 × 10³⁴cm⁻²s⁻¹. The cross section of Higgs production is about 0.2 pb in CEPC with $ \sqrt{s} $ of 250 GeV. The primary production process is via $ ZH $ production (96.6%), which is often referred as Higgs-strahlung [38-40] process, while the fraction of production via $ WW $-fusion [41-45] and $ ZZ $-fusion is 3.06% and 0.29% respectively. After ten years of running, one million of the Higgs boson (5000 fb⁻¹ collision data) are expected to be collected at the CEPC.

The work presented here focuses on the Higgs production in association with a pair of leptons ($ e^{+}e^{-} $ or $ \mu^{+}\mu^{-} $), in which Higgs decays to a pair of b-quark, c-quark or gluons. The leptons are either from Z-boson decay in Higgs-strahlung process, or ZZ-fusion which takes place only in $ e^{+}e^{-}H $ channel. The measurements of signal cross sections, denoted as $ \sigma^{H\to b\bar{b}}_{l^{+}l^{-}H} $, $ \sigma^{H\to c\bar{c}}_{l^{+}l^{-}H} $ and $ \sigma^{H\to gg}_{l^{+}l^{-}H} $, are described. The branching fraction of $ H\to b\bar{b}/c\bar{c}/gg $ can be derived once the cross sections of the $ l^{+}l^{-}H $ production, $ \sigma_{e^{+}e^{-}H} $ and $ \sigma_{\mu^{+}\mu^{-}H} $, are determined from other measurements. The work presented in this paper is partly inspired by the $ H\to b\bar{b}/c\bar{c}/gg $ analysis in ILC [46], and it is the subsequent study in $ H\to b\bar{b}/c\bar{c}/gg $ analysis presented in CEPC Higgs white paper [47] with improvement in background estimation.

This paper is organized as follows. After the introduction, a brief description of the detector will be presented in Section 2. The MC samples and event selections are described in Section 3 and Section 4 respectively. In Section 5, the flavor tagging and flavor-template-recoil-mass-fit, which is the procedure to extract signal event yields, is presented. In Section 6, the uncertainties of the signal cross sections and signal branching fractions are discussed, and finally in Section 7 a short summary is provided.

The detailed description of the proposed CEPC detector can be found elsewhere [37]. A vertex detector with high pixel resolution is located in the inner most part of the detector. The 6 layers of sensors are laid coaxially, in radius from 16 mm to 60 mm, covering 97% - 90% in the range of the polar angle. The spatial resolution in a single layer is 2.8 µm in the 2 inner layers and 4 µm in the 4 outer layers. The overall IP resolution can be represented as:

The first and second term in the right side of equation depict the resolution from finite single point position and the resolution due to multi-scattering respectively. These two types of resolution are parameterized as $ a $ and $ b $, which are estimated to be 5 and 10 respectively in CEPC. The parameters $ p $ and $ \theta $ are the momentum and polar angle of the reconstructed charged particle and $ \oplus $ denotes summation in quadrature. The high spatial resolution is essential to track impact parameter (IP) measurements and vertex reconstruction, on which the identification of the flavor of jets (flavor tagging) primarily relies.

A Time Project Chamber (TPC) is located outside of the vertex dector to take the major task of track measurement. It covers the solid angle up to $ \cos \theta = 0.98 $. When being operated in designed magnetic field of 3.5 T, the momentum resolution is $ \sigma(1/p_T) = 10^{-4}\;{\rm{GeV}} $.

The calorimeter system includes two sub-systems: the electromagnetic calorimeter (ECAL) and the hadronic calorimeter (HCAL). They are designed to have high energy resolution as well as high spatial resolution. The ECAL is a silicon-tungsten-based detector, which uses tungsten as absorber and silicon as sensor. It contains 30 layers of sampling structures. Each layer is divided into cells of $ 5\;{\rm{mm}}\times 5\;{\rm{mm}} $ in size. The HCAL includes 40 layers. Each layer contains a 20 mm thick stainless steel as absorber layer, a 3 mm thick glass resistive plate chamber (GRPC) and a 3 mm thick readout electronics in $ 1\;{\rm{cm}}\times 1\;{\rm{cm}} $ readout pad. The overall jet energy resolution(JER) is 3%-4%, and the two-jet system invariant mass resolution is required to be around 3%-4%, to distinguish the final states from $ Z $ boson and $ W^{\pm} $ boson hadronic decay. The high granularity of the ECAL and HCAL is crucial for application of the particle flow algorithm [48], which intends to reconstruct and identify each particle individually by combining information from all sub-detectors. The details of the ECAL and HCAL design and performance can be found in Ref. [49].

The muon system is mounted at the outermost part of the detector. The baseline design of muon detector has 8 sensitive layers in barrel and endcap region. On average, muons with momentum 2 GeV are expected to hit the first layer while those with momentum over 4 GeV can penetrate all the 8 layers. Muons with momentum above 5 GeV can be detected with standalone muon identification efficiency above 95% while the fake rate of pions is less than 1%. In order to improve the precision of the muon momentum measurement, the longitude and transverse position resolution are required to be $ \sigma _{z} = 1.5\; {\rm{cm}} $ and $ \sigma_{r\phi} = 2.0\; {\rm{cm}} $ respectively.

Both background and signal events are generated using Whizard [50] configured as no-polarization electron-positron collision with the $ \sqrt{s} $ of 250 GeV. The mass of the Higgs boson is assumed to be 125 GeV and the couplings are set to the SM predictions. The fragmentation and hadronization are performed by implementing PYTHIA6 [51]. All the MC datasets are normalized to the expected yields in data with integrated luminosity of 5000 fb^-1, by assigning a weight to the events of each process. The details of the event generation in CEPC can be found at Ref. [52].

The generated events undergo the detector simulation by Mokka [53], a GEANT4 [54] based detector simulator. The simulated hits are digitized and reconstructed with ArborPFA [55-57].

The charged particles are identified as electrons or muons by Lepton Identification in Calorimeter with High Granularity(LICH) algorithm, a dedicate lepton identification algorithm designed for Higgs factories, as described in Ref. [58]. The algorithm use the $ {\rm d}E/{\rm d}X $ information measured by TPC, together with shower and hit information in high granularity calorimeter, as discrimination variables. A Boosted Decision Trees algorithm [59] with Gradient boosting (GBDT) method is implemented to further extract the discriminative characteristics of the variables. The overall efficiencies for electron and muons are 99.7% and 99.9% respectively, with the rate of electron and muons misidentified as each other smaller than 0.07%. The rate of particles like $ \pi^{\pm} $ identified as electrons or muons is 0.21% and 0.05% respectively.

Jets reconstruction and flavor tagging are essential to this analysis. They are done with the LCFIPLUS [60] software package, integrating the functionality of vertex finding, jet reconstruction and jet flavor tagging. Before the jet clustering, the secondary vertices are identified based on the reconstructed tracks. Jets are reconstructed by Durham algorithm [61]. This algorithm begins with jet cluster candidates, which are either single reconstructed particles, or compound objects like reconstructed secondary vertices. The procedure iteratively pairs the clusters and calculates the distance between them, defined as $ y_{ij} = \min\{E^2_i,E^2_j\}(1-\cos\theta_{ij})/E^2_{vis} $, where $ E_i $ and $ E_j $ are the energy of $ i $-th and $ j $-th cluster, and $ \theta_{ij} $ refers to the angle between them. $ E_{vis} $ are the sum of energy of all the clusters in the event. Clusters with minimum $ y_{ij} $ are merged, reducing the cluster number by 1, until the ceasing criteria are met. The ceasing criteria can be either a minimum $ y_{ij} $-value threshold, or the remaining clusters number equals to the required jet multiplicity. In this analysis, the $ y_{ij} $ threshold is set to 0, and each event is forced to have two jets reconstructed. The minimum value $ y_{ij} $ can be denoted as $ Y_k $, in which $ k $ is the multiplicity of cluster candidates. When $ k $ is larger than 2, the signal events have relatively smaller $ Y_k $ than the backgrounds with at least $ k $ primary parton^①. This is because in a signal event, the closet two clusters, among more of them, are likely to be from the same parton. They tend to be collinear, which lead to small $ Y_k $.

The final state of the signal contains two jets and two leptons with opposite charge and same flavor. There are two types of backgrounds according to the final states: the irreducible backgrounds and the reducible backgrounds. The irreducible backgrounds contain the same final states as that in signal. The semi-leptonic $ ZZ $ process, in which one $ Z- $boson decays to $ e^{+}e^{-} $ or $ \mu^{+}\mu^{-} $ and the other decays to quark pair, is a typical example and the major components of the irreducible backgrounds. The reducible backgrounds include all the other types of background which have different final states, such as hadronic or leptonic $ W^+W^- $ and $ ZZ $ production, the lepton pair or quark pair production, semi-leptonic $ W^+W^- $ production, and the Higgs production processes other than $ l^{+}l^{-}H $. There are backgrounds from $ e^{+}e^{-}H $ or $ \mu^{+}\mu^{-}H $ production, in which the Higgs boson decays to $ WW^* $ or $ ZZ^* $ and subsequently both vector bosons decay to quarks. They are reducible backgrounds since there are more quarks in final states than that in signal. But they behave more like irreducible background experimentally, so they will be discussed separately from the typical irreducible or reducible backgrounds, and they will be referred as $ e^{+}e^{-}H $ or $ \mu^{+}\mu^{-}H $ background in this analysis. There are also other irreducible backgrounds from $ l^{+}l^{-}H $, in which the Higgs boson decays to light quark pair, to $ \tau^+\tau^- $ or to photons when all the tauons or photons from the Higgs boson decay are misidentified as jets. Their contributions are expected to be very small and they are classified as reducible backgrounds instead of $ l^{+}l^{-}H $ backgrounds.

Each event must contain two isolated tracks with opposite charge, reconstructed as $ e^{+}e^{-} $ or $ \mu^{+}\mu^{-} $. The energy of each isolated lepton candidate must be above 20 GeV. The isolation criterion requires $ E^2_{{\rm{cone}}}<4 E_{{\rm{lep}}}+12.2 $, where $ E_{{\rm{lep}}} $ is the energy of the lepton, and $ E_{{\rm{cone}}} $ is the energy within a cone $ \cos\theta_{{\rm{cone}}}>0.98 $ around the lepton. Here $ E_{{\rm{lep}}} $ and $ E_{{\rm{cone}}} $ are measured in GeV. Events with additional isolated leptons are rejected (extra isolated lepton veto). The polar angle of lepton pair system is required to be in the range of $ |\cos\theta_{\mu^{+}\mu^{-}}| < 0.81 $ and $ |\cos\theta_{e^{+}e^{-}}| < 0.71 $. The angle between the two isolation tracks $ \psi $ is required to satisfy $ \cos\psi > -0.93 $ and $ \cos\psi > -0.74 $ for $ e^{+}e^{-} H $ and $ \mu^{+}\mu^{-} H $ channel respectively, to reject events from lepton pair production where leptons tend to be back to back. The invariant mass of lepton pair is required to be inside the Z-mass window, which is defined as 77.5 GeV to 104.5 GeV.

The remaining particles in the event are used to reconstruct exactly two jets with polar angle $ \theta_{\rm jet} $ in the range of $ |\cos\theta_{\rm jet}|<0.96 $. The two jets are required to contain at least 20 particles, each with energy no less than 0.4 GeV, according to the optimization of distinguishing the signal events against the events including fake jets from photons or leptons. The invariant mass of the pair of jets is required to be between 75 GeV and 150 GeV to reject the irreducible backgrounds. To suppress the $ e^{+}e^{-}H $ and the $ \mu^{+}\mu^{-}H $ backgrounds, $ Y_4 $, an indicator of primary parton multiplicity described in Section 3, was required to be less than 0.011.

The invariant mass of the lepton pair's recoil system, denoted as $ M^{l\bar{l}}_{{\rm{recoil}}} $, provides clear signature of the $ l\bar{l}H $ events. The definition of $ M^{l\bar{l}}_{{\rm{recoil}}} $ is:

in which $ \sqrt{s} = 250\;{\rm{GeV}}$, while $ E $ and $ \vec{P} $ stand for the energy and momentum of the leptons respectively. A Higgs mass window is defined by requiring $ M^{l\bar{l}}_{{\rm{recoil}}} $ between 124 GeV and 140 GeV. The signal and background yields during the event selections are summarized in Table 1 for $ \mu^{+}\mu^{-}H $ and $ e^{+}e^{-}H $ analysis, respectively.

$\mu^{+}\mu^{-}H\to \mu^{+}\mu^{-}+b\bar{b}/c\bar{c}/gg$ channel
	signals	$\mu^{+}\mu^{-}H$ background	other Higgs background	irreducible background	other SM background
original	2.45×104	1.10×104	1.01×106	1.05×106	4.96×108
lepton pair selection without recoil mass cut	1.51×104	6.56×103	227	1.09×104	2.79× 104
jets pair selection and lepton pair recoil mass cut for fit region	1.32×104	1.80×103	108	7.75×103	43.6
signal region	1.31× 104	1.80× 103	96.1	5.78× 103	38.4
$e^{+}e^{-}H\to e^{+}e^{-}+b\bar{b}/c\bar{c}/gg$ channel
	signals	$e^{+}e^{-}H$ background	other Higgs background	irreducible background	other SM background
original	2.63× 104	1.17× 104	1.01× 106	1.62× 106	4.95× 108
lepton pair selection without recoil mass cut	9.17× 103	3.53× 103	128	9.00× 103	7.11× 104
jets pair selection and recoil lepton pair mass cut of fit region	7.14× 103	917	56.1	8.63× 103	69.4
signal region	7.13× 103	913	36.4	4.14× 103	67.4

Table 1.  Event yields of cut flow. Signal events are $l\bar{l}+H\to l\bar{l}+b\bar{b}/c\bar{c}/gg$ combined. $\mu^{+}\mu^{-}H$ and $e^{+}e^{-}H$ background refers to the background which Higgs are produced associated with $\mu^{+}\mu^{-}$ and $e^{+}e^{-}$, but decay to final states other than $b\bar{b}/c\bar{c}/gg$. 'Other Higgs background' stands for the Higgs production process different from the signal. 'Irreducible SM background' is the $e^{+}e^{-}/\mu^{+}\mu^{-}+$jet pair process without Higgs productions. 'Other SM background' includes all the other background processes. The 'fit region' will be described in Section 5.

After applying the object and event selection described in Section 4, it is necessary to get the component fractions of $ H\to b\bar{b} $, $ H\to c\bar{c} $ and $ H\to gg $ processes. It is achieved by the high performance multi-variable-based flavor tagging toolkit in LCFIPLUS. This toolkit is responsible for the flavor tagging algorithm training and implementing. The training is applied to the simulated $ Z\to q\bar{q} $ sample, produced with $ \sqrt{s} $ of 91.2 GeV. The reconstructed jets in the sample are classified to 4 categories according to multiplicity the secondary vertex and lepton in the jet: jets with both secondary vertex and lepton, jets with secondary vertex but without lepton, jets without secondary vertex but with lepton and jets without neither secondary vertex nor lepton. The details of the secondary vertex finding can be find in Ref. [60]. The leptons are selected according to lepton identification presented in Section 3 without any isolation criteria required. In each category, two types of training, one for the $ b $-tagging algorithm and the other for the $ c $-tagging algorithm, are implemented with GBDT method. Discrimination variables such as jets kinematic variables, impact parameters of tracks and secondary vertex parameters are included in the training. After the training, a $ b $-tagging model and a $ c $-tagging model are created. By invoking these models, a $ b $-jet likeness weight and a $ c $-jet likeness weight are calculated for each jet, representing the resemblance of the jet to a $ b $-jet or a $ c $-jet respectively. The likeness weights are in the range between 0 and 1, and a higher weight indicates higher likelihood of a jet to be $ b $-jet or $ c $-jet.

The $ b $-weight likeness of the two individual jets of any selected event, $ L_{b1} $ and $ L_{b2} $, can be used to construct the combined $ B $ likeness, defined as:

A combined $ C$-likeness $ X_C $ can be defined in similar way, by replacing the $ L_{b1} $ and $ L_{b2} $, by $ L_{c1} $ and $ L_{c2} $ respectively, in the right side of Eq. (3). The conservation of quark flavor in the Higgs boson decay guarantee that $ X_B $($ X_C $) is close to 1 if the Higgs boson decay to $ b\bar{b} $ ($ c\bar{c} $) while close to 0 otherwise. Thus the flavor of the jets in each event can be characterized by the two dimensional distribution of variable $ X_B $ and $ X_C $. Flavor templates are created for the $ (X_B,X_C) $ distributions of different processes. By fitting the data with these flavor templates, one can get the fraction of each process. This template fit approach was implemented in ILC $ H\to b\bar{b}/c\bar{c}/gg $ analysis [46].

Here the template fit is combined with the fit to $ M^{l\bar{l}}_{{\rm{recoil}}} $, defined in Eq. (2). Such a combined fit is motivated by separating the irreducible backgrounds from signal while extracting the flavor components in signal events. The fit is applied to three variables: $ M^{l\bar{l}}_{{\rm{recoil}}} $, $ X_B $ and $ X_C $ and simultaneously. The RooFit package is implemented to perform an unbinned likelihood fit to the weighted events. The overall likelihood function is constructed as:

The $ N^{{\rm{sig}}}_{H\to b\bar{b}} $, $ N^{{\rm{sig}}}_{H\to c\bar{c}} $ and $ N^{{\rm{sig}}}_{H\to gg} $ are the signal event yields of $ H\to b\bar{b}/c\bar{c}/gg $ respectively, which are the concern of this analysis. $ N^{{\rm{bkg}}}_{{\rm{irred}}\_bb} $, $ N^{{\rm{bkg}}}_{{\rm{irred}}\_cc} $ and $ N^{{\rm{bkg}}}_{{\rm{irred}}\_uds} $ are the event yields of the irreducible background with $ b\bar{b} $, $ c\bar{c} $ or light quark in final states, respectively. $ N_{{\rm{redu}}} $ is the event yield of reducible backgrounds. $ \vec{\theta_s} $ and $ \vec{\theta_b} $ are the shape parameters of $ M^{l\bar{l}}_{{\rm{recoil}}} $ spectrum of $ l\bar{l}H $ and irreducible background respectively. Functions like $ P^{{\rm{p}}}_{{\rm{flavor}}}(X_B,X_C) $ are the two-dimensional distributions of $ (X_B,X_C) $ of process $ p $, which are modeled as two-dimensional histograms with $ 20\times 20 $ bins from MC simulation. The recoil mass function of $ l^{+}l^{-}H $ events, denoted as $ P_{{\rm{sig}}}(M^{l\bar{l}}_{{\rm{recoil}}};\vec{\theta_s}) $, are described by a crystal ball function, combined with a double sided exponential function which has peak near the Higgs mass threshold, to describe the resolution effect of track energy and momentum measurements. The lepton pair recoil mass spectrum of irreducible background events, denoted as $ P_{{\rm{irred}}}(M^{l\bar{l}}_{{\rm{recoil}}};\vec{\theta_b}) $, is described by a first order Chebychev polynomial function ($\mu^{+}\mu^{-}H$ channel) or exponential functions ($e^{+}e^{-}H$ channel). $ P_{{\rm{redu}}} $ is the 3-dimensional distribution of $ M^{l\bar{l}}_{{\rm{recoil}}} $, $ X_B $ and $ X_C $ of the reducible backgrounds. It is depicted by MC simulation. The $ M^{l\bar{l}}_{{\rm{recoil}}} $ ranges in the fit are set to be between 120 GeV and 140 GeV for $\mu^{+}\mu^{-}H$ channel, and between 115 GeV and 140 GeV for $e^{+}e^{-}H$ channel. They are slightly wider than the signal regions to have better estimation of irreducible background. The signal and background yields in these regions are also summarized in Table 1. All the functions of $ P_{{\rm{sig}}}(M^{l\bar{l}}_{{\rm{recoil}}};\vec{\theta_s}) $, $ P_{{\rm{irred}}}(M^{l\bar{l}}_{{\rm{recoil}}};\vec{\theta_b}) $, $ P^{{\rm{p}}}_{{\rm{flavor}}}(X_B,X_C) $ and $ P_{{\rm{redu}}}(M^{l\bar{l}}_{{\rm{recoil}}},X_B,X_C) $ are normalized to 1 in the fit range.

The shape parameters of the crystal ball function in $ \vec{\theta_s} $ and all the parameters in $ \vec{\theta_b} $ are free in the fit. The tail distribution of signal events, in $ M^{l\bar{l}}_{{\rm{recoil}}} $ spectrum near the Higgs mass threshold, are fixed according to a signal-only-fit. The event yields parameters such as $ N^{{\rm{sig}}}_{H\to b\bar{b}} $, $ N^{{\rm{sig}}}_{H\to c\bar{c}} $, $ N^{{\rm{sig}}}_{H\to gg} $, $ N^{{\rm{bkg}}}_{{\rm{irred}}\_b\bar{b}} $, $ N^{{\rm{bkg}}}_{{\rm{irred}}\_c\bar{c}} $ and $ N^{{\rm{bkg}}}_{{\rm{irred}}\_uds} $ are also free in the fit, while $ N^{{\rm{bkg}}}_{l^{+}l^{-}H} $ and $ N_{{\rm{redu}}} $ are fixed to MC predictions. The 3D distribution of reducible background $ P_{{\rm{redu}}}(M^{l\bar{l}}_{{\rm{recoil}}},X_B,X_C) $ is also fixed as that predicted by MC simulation.

The fit program minimizes the negative logarithm of the likelihood function $ -\sum_i w_i \log L(M_i,X^i_B,X^i_C) $, in which $ L $ is the likelihood function presented in Eq.(4); $ w_i $, $ X^i_B $, $ X^i_C $ and $ M_i $ are event weight^②, $ X_B $, $ X_C $ and $ M^{l\bar{l}}_{{\rm{recoil}}} $ of event $ i $ respectively. The summation is applied by including all the events in the fit region. The fit results of the simulated data are shown in Fig. 1(a), in which one can find that the model describes the data very well.

Figure 1.  (color online) 3D-fit result projected on three dimension: (a) fit result projected on recoil mass distribution in $ \mu^{+}\mu^{-}H $ channel, (b) fit result projected on recoil mass distribution in $ e^{+}e^{-}H $ channel, (c) fit result projected on B-likeness distribution in $ \mu^{+}\mu^{-}H $ channel, (d) fit result projected on B-likeness distribution in $ e^{+}e^{-}H $ channel, (e) fit result projected on C-likeness distribution in $ \mu^{+}\mu^{-}H $ channel, (f) fit result projected on C-likeness distribution in $ e^{+}e^{-}H $ channel.

The statistical uncertainty was estimated by using toyMC method that includes 6000 ($\mu^{+}\mu^{-}H$ channel) and 10000 ($e^{+}e^{-}H$ channel) iterations. In each iteration, the 'data' is filled in a 3D histogram with dimension of $ (M^{l\bar{l}}_{{\rm{recoil}}},X_B,X_C) $. In each bin of the histogram, the event yields fluctuated according to a Poisson distribution. Binned fit with the model described in Section 5 is then applied to the fluctuated histogram. The statistical uncertainty of signal event yields can be estimated according to dispersion of fit results of the toy MC test. The results of toy MC test for $ H\to b\bar{b} $, $ H\to c\bar{c} $ and $ H\to gg $ are represented in Fig. 2(a), in terms of pull of fitted signal events number for toy MC samples in the signal region which conform well to the standard normal distribution.

Figure 2.  (color online) Toy MC test result in terms of pull of signal strength in each channel. The pull distributions are fitted with Gaussian function.

The systematic uncertainties from luminosity, lepton identification and selection efficiencies, $ Z\to\mu^{+}\mu^{-}/e^{+}e^{-} $ modeling and ISR correction factor [62] are included in the measurements of $ \sigma^{H->b\bar{b}/c\bar{c}/gg}_{l^{+}l^{-}H} $. However, these uncertainties are also included in the measurement of inclusive Higgs boson production cross section associated with lepton pair $ \sigma_{l\bar{l}H} $, like that presented in Ref. [63]. To get the branching fractions of $ H\to b\bar{b}/c\bar{c}/gg $, one need to divide the measured $ \sigma^{H->b\bar{b}/c\bar{c}/gg}_{l^{+}l^{-}H} $ by $ \sigma_{l\bar{l}H} $. As a consequence the systematic uncertainties discussed above will be canceled. So here we can ignore these uncertainties.

The fit method described in Section 5 has two types of systematic uncertainty. The first kind is due to imperfect modeling of the PDF in the likelihood in Eq. (4). The inaccurate modeling of $ M^{l\bar{l}}_{{\rm{recoil}}} $ distribution and the bias in the prediction of $ F_{{\rm{flavor}}}(X_B,X_C) $ can lead to such kind of uncertainty. The latter will be discussed later and here we only focus on the recoil mass modeling. The recoil mass functions for signal and background are verified by fitting the signal and background datasets alone respectively. These results demonstrate that the function describes the shape very well. The shape parameters are left to float in the fit, except for those of the tail function in signal recoil mass. These tail function can be studied by comparing the track momentum and energy resolution in data and MC. So far we assume the resolution is well simulated in MC sample.

The other kind of systematic uncertainty in the fit comes from the uncertainty of fixed parameters, including the event yields of $ e^{+}e^{-}H $ or $ \mu^{+}\mu^{-}H $ background, as well as the yields of all irreducible background. We conservatively set $ H\to WW^* $ and $ H\to ZZ^* $ event yields 5% higher and lower than the MC prediction, and vary the yields of non-ZZ backgrounds by ±100%, to estimate the corresponding systematic uncertainty. The systematic uncertainties discussed above are included in the row of 'Fixed Background' in Table 2.

Higgs boson production	$\mu^{+}\mu^{-}H$				$e^{+}e^{-}H$
Higgs boson decay	$H\to b\bar{b}$	$H\to c\bar{c}$	$H\to gg$		$H\to b\bar{b}$	$H\to c\bar{c}$	$H\to gg$
statistic uncertainty	1.1%	10.5%	5.4%		1.6%	14.7%	10.5%
fixed background	−0.2%	+4.1%	7.6%		−0.2%	+4.1%	7.6%
	+0.1%	−4.2%			+0.1%	−4.2%
event selection	+0.7%	+0.4%	+0.7%		+0.7%	+0.4%	+0.7%
	−0.2%	−1.1%	−1.7%		−0.2%	−1.1%	−1.7%
flavor tagging	−0.4%	+3.7%	+0.2%		−0.4%	+3.7%	+0.2%
	+0.2%	−5.0%	−0.7%		+0.2%	−5.0%	−0.7%
combined systematic uncertainty	+0.7%	+5.5%	+7.6%		+0.7%	+5.5%	+7.6%
	−0.5%	−6.6%	−7.8%		−0.5%	−6.6%	−7.8%

Table 2.  Uncertainties on $\sigma^{b\bar{b}}_{l^{+}l^{-}H}$, $\sigma^{c\bar{c}}_{l^{+}l^{-}H}$ and $\sigma^{gg}_{l^{+}l^{-}H}$.

In signal events extra leptons (leptons other than the 2 primary leptons) comes from leptonic decay of heavy flavor quarks. Systematic uncertainty of extra isolated lepton veto efficiency for $ H\to b\bar{b} $ and $ H\to c\bar{c} $ can be estimated using the $ b\bar{b} $ and $ c\bar{c} $ events produced at $ Z $-pole respectively. By selecting the two-jets events, requiring both of the two jets tagged as $ b $-jets and looking for an isolated lepton, the efficiency of lepton veto can be studied with 2 billion $b\bar{b}$ events in the Z-pole sample. The veto efficiency can be studied in a precision so high that it has no visible impact to the analysis presented, due to the high statistics of Z-pole sample and very low rate of occurrence of extra leptons. The impact to $ H\to c\bar{c} $ can be studied in similar way. By requiring both jets to be $ c $-tagged in $ Z- $pole sample and looking for isolated leptons, one can also study the impact with a precision high enough to have any visible impact to the analysis. For $ H\to gg $ we assume the MC have ±50% of uncertainty in predicting the lepton-veto efficiency, and estimate the impact to be negligible in this analysis.

The jets' particle multiplicity, jet angular distribution and $ Y_k $ value can be studied in very high precision with high statistics $ Z- $pole data. Correspondingly, the systematic uncertainties of the efficiencies of jets' particle multiplicity cut, jet $ \cos \theta $ cut and $ Y_k $ are negligible.

The systematic uncertainty of the efficiency of jet pair invariant mass cut can be estimated from the jet energy resolution. We apply a smearing on jet pair mass distribution according to a gaussian distribution corresponding to the jet energy resolution. We take the value of 4% as the jet energy resolution from the CEPC pre-CDR [37] and calculate the uncertainty on the event yields of $ H\to b\bar{b} $, $ H\to c\bar{c} $ and $ H\to gg $ are $ ^{+0.68\%}_{-0.20\%} $, $ ^{+0.43\%}_{-1.08\%} $ and $ ^{+0.71\%}_{-1.68\%} $ respectively. The uncertainties of extra lepton veto, the jet angular and reconstructed particles multiplicity and jet pair mass resolution are included in the row of 'Event Selection' in Table 2.

Since flavor tagging method is implemented via flavor template fit, the flavor tagging systematic uncertainty is directly caused by the difference between the templates from the MC prediction and the templates in data. Evaluating such differences demands delicate flavor tagging commissioning and calibration. Although no such commissioning or calibration has been done yet, we can estimate the systematic uncertainty by assuming a difference between data and MC after the calibration was applied, and this difference is subsequently studied in terms of its impact on the $ H\to b\bar{b}/c\bar{c}/gg $ branching fractions measurement. We select $ ZZ\to q\bar{q}+\mu^{+}\mu^{-} $ events as a control sample, which has a purity of 99.6%, and assuming a data-MC comparison has been done on the template distribution on this control sample. The estimation of the data-MC agreement is limited by the statistic uncertainty of the control sample, and the knowledge of the flavor components of hadronic $ Z- $decay. For example, more than 80% of the $ b\bar{b} $ events is concentrated in the region with highest $ b $-likeness($ b $-likeness>0.95) and lowest $ c $-likeness($ c $-likeness<0.05). If, due to some kinds of $ b- $tagging systematic uncertainty, the $ Z\to b\bar{b} $ events fraction in this most concentrate bin changed, the data-MC disagreement would be increase. There are $ 1.92\times 10^4 $ $ Z\to b\bar{b} $ in this bin, which has a statistic uncertainty 0.72%. The current combined measurements of $ R_b $, defined as $ Br(Z\to b\bar{b})/Br(Z\to q\bar{q}) $, has the uncertainty of 0.31% [64]. So the data and MC can be compared in the precision of $ \sqrt{0.72{\text{%}}^2+0.31{\text{%}}^2} = 0.78{\text{%}} $. Scaling the contents in this bin up and down by 0.78% in the $ b\bar{b} $ template, one can estimate the uncertainty to $ H\to b\bar{b} $, $ H\to c\bar{c} $ and $ H\to gg $ are $ ^{-0.4\%}_{+0.2\%} $, $ ^{+3.7\%}_{-5.0\%} $ and $ ^{+0.2{\text{%}}}_{-0.7{\text{%}}} $ respectively. Had we use a much larger control sample (hadronic events at the $ Z $-pole), and had a better understanding on the relationship between the flavor tagging variables and kinematic feature, the uncertainty will be further reduced.

The statistical uncertainty and systematic uncertainty discussed in this section are summarized in Table 2.

The measurements of branching fractions of $ H\to b\bar{b}/c\bar{c}/gg $ are studied in $ \mu^{+}\mu^{-}H $ and $ e^{+}e^{-}H $ process, in the scenario of analysing 5000 fb^-1 $ e^{+}e^{-} $ collision data with $ \sqrt{s} $ of 250 GeV in CEPC. The statistical uncertainty on $ \sigma^{b\bar{b}}_{l^{+}l^{-}H} $, $ \sigma^{c\bar{c}}_{l^{+}l^{-}H} $ and $ \sigma^{gg}_{l^{+}l^{-}H} $ measurements are estimated around 1.1%, 10.5% and 5.4% respectively in $ \mu^{+}\mu^{-}H $ channel, and 1.6%, 14.7% and 10.5% respectively in $ e^{+}e^{-}H $ channel. The systematic uncertainties on the branching fraction measurements are also studied, which are around 0.6%, 6% and 8% for $ b\bar{b} $, $ c\bar{c} $ and $ gg $ final states respectively. The high precision of this measurement benefits from the distinct signature of events with the Higgs boson and clean background in electron-positron collider, as well as the model independent analysis. As a comparison, by combining the extrapolated results in ATLAS and CMS in the scenario of the High Luminosity LHC(HL-LHC) [65], the $ H\to b\bar{b} $ branching fraction precision is expected to be 4.4%, in which the statistic uncertainty, the systematic uncertainty and theoretical uncertainty are 1.5%, 1.3% and 4.0% respectively. This study demonstrates the feasibility of precise measurement of Higgs Yukawa coupling to quarks at the CEPC.

We would like to thank the CEPC higgs physics working group for the valuable discussions as well as the software and physics work without which this work couldn't have been accomplished.