Charged lepton flavor violating Higgs decays at future $e^+e^-$ colliders

After the discovery of the Higgs boson, several future experiments have been proposed to study the Higgs boson properties, including two circular lepton colliders, the CEPC and the FCC-ee, and one linear lepton collider, the ILC. We evaluate the precision reach of these colliders in measuring the branching ratios of the charged lepton flavor violating Higgs decays $H\to e^\pm\mu^\mp$, $e^\pm\tau^\mp$ and $\mu^\pm\tau^\mp$. The expected upper bounds on the branching ratios given by the circular (linear) colliders are found to be $\mathcal{B}(H\to e^\pm\mu^\mp)<1.2\ (2.1) \times 10^{-5}$, $\mathcal{B}(H\to e^\pm\tau^\mp)<1.6\ (2.4) \times 10^{-4}$ and $\mathcal{B}(H\to \mu^\pm\tau^\mp)<1.4\ (2.3) \times 10^{-4}$ at 95\% CL, which are improved by one to two orders compared to the current experimental bounds. We also discuss the constraints that these upper bounds set on certain theory parameters, including the charged lepton flavor violating Higgs couplings, the corresponding parameters in the type-III 2HDM, and the new physics cut-off scales in the SMEFT, in RS models and in models with heavy neutrinos.


I. INTRODUCTION
The discovery of the Higgs boson [1,2] not only completes the Standard Model (SM) of particle physics, but also opens new windows for searches for new physics through the Higgs portal.
The CLFV Higgs decays are interesting, because their observation may provide insight into some fundamental questions in nature, e.g., whether there is a secondary mechanism for the electroweak symmetry breaking [41], why the neutrino masses are tiny [42], and whether there is an extra dimension responsible for the gauge hierarchy generation [43]. They have thus attracted a lot of attention from both theorists and experimentalists . The CMS collaboration reported the first hint of charged lepton flavor violation in the H → µ ± τ ∓ channel with a significance of 2.4 standard deviations [63,64]. Although this signal disappeared later [63,64], the CLFV Higgs decays are still worthy to be studied with higher precision. On one hand, the so-called flavor anomalies, which indicate lepton flavor non-universality, reported by the B factories and the LHCb collaboration [65][66][67][68], in some sense also imply lepton flavor violation (when the lepton mass matrices are diagonalised to obtain the physical states, unequal diagonal couplings with different leptons will lead to off-diagonal couplings). On the other hand, the B meson decay channels in which the flavor anomalies are observed are always polluted by complicated strong dynamics, while the much cleaner CLFV Higgs decay channels will provide a better chance to study the mechanism generating the lepton flavor violation or non-universality once they are discovered. The potential to search for such decay channels at the High Luminosity LHC (HL-LHC) has been estimated in [69]. By this paper, we study the sensitivity of the three lepton colliders in measuring the CLFV Higgs decays based on the detector simulation of the signal events of the decay channels and the corresponding background. There are already some studies on the ILC measurement of H → µ ± τ ∓ [69][70][71][72], and the difference between our paper and theirs will be discussed.
The paper is organized as follows: in Section II, we perform the detector simulation of the signal and background events for the three CLFV Higgs decay channels at the CEPC and obtain the upper bounds on the three decay rates, based on which we estimate the corresponding upper bounds expected to be given by the FCC-ee and the ILC; we derive in Section III the constraints on theory parameters including the CLFV Higgs couplings, the relevant parameters in the type-III two-Higgs-doublet-model (2HDM) and the new physics cut-off scales in the SM effective field theory (SMEFT), in Randall-Sundrum (RS) models and in models with heavy neutrinos; we summarize by Section IV.

II. SIMULATION AND ANALYSIS
In this section, we evaluate the possible reach of the three colliders with √ s = 240 -250 GeV, the CEPC, the FCC-ee and the ILC, in measuring the branching ratios of the three CLFV Higgs decays. Investigating the dominant Higgs production process e + e − → ZH, with considered only the Z boson decaying hadronically into a quark pair, we expect that the signal events each contain two charged leptons of different flavors and two jets. Therefore, the four-fermion processes will form the major SM background. The signal processes e + e − → Z(→ qq)H(→ e ± µ ∓ , e ± τ ∓ , µ ± τ ∓ ) are simulated via MadGraph v2.5.2 [73], with the corresponding three CLFV Higgs vertices implemented. The background events are generated via WHIZARD 2.5.0 [74,75]. PYTHIA 6.4 [76] is then used to manage hadronization and parton showers for both the signal and background events.
Finally, Delphes 3.4.1 [77,78] is adopted for detector simulation. Note that we only carry out the above simulation procedure for the CEPC with √ s = 240 GeV and an integrated luminosity of 5 ab −1 . We suppose that the FCC-ee with also √ s = 240 GeV has a same integrated luminosity and similar detector performance, and hence the CEPC results will apply to the FCC-ee. Different from the CEPC and the FCC-ee, the ILC is planned to run at We find that the beam polarization will considerably change neither the angular distribution of the signal final states nor the statistical uncertainties owing to the background, which makes it possible to estimate the ILC sensitivity to the CLFV Higgs decays based on the CEPC simulation, as discussed in detail for each channel below in this section.
dddd, ddbb, bbbb, uūss, uūbb and ssbb; SW ql includes e + νecs; WW ql includes τ + ντūd, µ −ν µcs and µ −ν µud. For each CLFV Higgs decay channel, the production of 10000 signal events at the CEPC are simulated. We simulate the four fermion background events by WHIZARD with an integrated luminosity of 5 ab −1 , giving the cross sections and the numbers of events for different categories in Table I. The four fermion processes are classified into different categories according to their final states [79] as follows. We divide the processes into two groups. The first group contains the processes without (anti-) electron or (anti-) electron neutrino in their final states, and the processes in the second group have at least one (anti-) electron or (anti-) electron neutrino in their final states.
We start with a process whose final state is constituted by two pairs of mutually charge conjugate fermions like uūµ + µ − or uūe + e − that can not arise from decays of two W bosons. If this process belongs to the first group, we will classify it into the "ZZ" category; if it belongs to the second group, it will be classified into "Single Z" (SZ). A process belonging to the first or the second group with two pairs of mutually charge conjugate fermions in the final state will be classified into the "ZZ or WW" (ZZWW) category or the "Single Z or Single W" (SZSW) category, respectively, if its final state can arise from decays of two W bosons. The remaining processes in the first and second group are classified into "WW" and "Single W" (SW), respectively. More details can be found in [12]. We further divide each category of processes into different sets, according to whether the final states contain only quarks ("qq"), only leptons ("ll") or both of them ("ql"). For example, uūµ + µ − is classified into "ZZ ql", while ν µνµ µ + µ − is classified into "ZZWW ll". Note we only consider the four fermion processes listed in the footnote of TABLE I. The neglected processes do not contribute to the background after the chosen cuts are made.
In the following, for each CLFV Higgs decay channel, we set appropriate cuts on both the generated signal and background events. Owing to the fact that the signal events of different channels need to be reconstructed differently and thus suffer from different backgrounds, different event selection cuts are determined through the significance optimization for different channels.
Based on the signal detection efficiencies and the background event numbers, the CEPC bounds on the decay branching ratios are evaluated first, which also apply to the FCC-ee. After that, we estimate the corresponding ILC bounds by scaling the luminosity times the signal and background cross sections, i.e., scaling the event numbers, for each polarization option.
To reconstruct the signal events of e + e − → HZ → e ± µ ∓ Z, we select the events with final states containing one electron, one muon, and two jets to reconstruct the Z boson, with the di-lepton In with B(Z → qq) ≈ 70% and N H = 1.05 ×10 6 the number of the Higgs bosons to be produced by the CPEC, we evaluate the upper bound on the H → e ± µ ∓ branching ratio, This bound is also accepted for the FCC-ee.
As for the ILC, the ratios of the ZH production cross sections with the four polarization options to that of the CEPC are Recalling the integrated luminosities for the four options, we obtain that about 5.21 × 10 5 Higgs bosons will be produced at the ILC. Through a tree-amplitude analysis of the e + e − → Z(→ jj)H(→ ± ∓ ) process, we find that the beam polarization will not change the angular distribution of the final state, and thus the signal detection efficiency = 41.15% is also valid under the same event selection condition. The only one background event listed in Table II arises from the process e + e − τ + τ − . We calculate the its production cross section with each polarization option as what we did for the signal, and find that the background event number after 2 ab −1 data collected is One might be unsatisfied with that we estimate the background by scaling the total event numbers, since the beam polarization might change the e + e − τ + τ − angular distribution. We then test the polarization impact by assuming an extreme situation where the role of the polarization is so important that the background event number at the ILC is 0. Even in such a case, the upper bound is then only reduced to 1.9×10 −5 , by less than 10% of magnitude. The key point is that the CEPC background is already very small, so further reduction of the background will not optimize the upper bound significantly. Therefore, we conclude that, regarded as estimates, our results for the ILC upper bounds (also for the other two channels) are acceptable.
The upper bounds on B(H → e ± µ ∓ ) expected to be given by the CEPC, the FCC-ee and the ILC are displayed in FIG. 1, together with the present upper bound B(H → e ± µ ∓ ) < 0.035% at 95% CL reported by the CMS collaboration [81]. We find that the three future lepton colliders are To reconstruct the signal events of e + e − → HZ → e ± τ ∓ Z, we again first select the events with two jets, one electron and one muon in their final states. In each event, the Z boson is reconstructed using the two jets as in the H → e ± µ ∓ case, the tau lepton is reconstructed from the muon and the missing energy, and the Higgs boson is finally reconstructed from the electron and the tau lepton. Our tau reconstruction method has an efficiency that will allow us to estimate the CEPC sensitivity to B(H → e ± τ ∓ ), since the branching ratio of a tau lepton decaying into a muon and two neutrinos B(τ → µνν) is nearly 20%. Our method also greatly suppresses the background, e.g., if we choose to reconstruct the tau lepton from electron and missing energy, then the processes e + e − → e + e − qq would give a large background. The above event selection method requires the following cuts: 66 GeV < m jj < 94 GeV, m µE M < 4 GeV and 121 GeV < m eτ < 130 GeV, where m µE M is the invariant mass of the muon and missing energy, and m eτ is the invariant mass of the electron, muon and missing energy. We further set a cut on the electron pseudorapidity of |η e | < 2 to suppress the background arising from the SZ processes. The cuts are summarized in To reconstruct the signal events of e + e − → HZ → µ ± τ ∓ Z, we still select the events containing one electron, one muon and two jets in their final states. In each event, the Z boson is reconstructed from the two jets, the tau lepton is reconstructed from the electron and the missing energy, and the Higgs boson is reconstructed from the muon and the reconstructed tau lepton. Analogous to the H → e ± τ ∓ case, we do not consider other ways to reconstruct the tau lepton. We employ the According to another study of H → µ ± τ ∓ at the ILC [71], the upper bound on the branching ratio is given as B(H → µ ± τ ∓ ) < 2.9 × 10 −5 at 95% CL if 90% signals survive the event selection.
This bound is more stringent than that obtained in this work. A possible reason is that we only considered the cleanest method to reconstruct the tau lepton, as described previously, which makes our evaluated upper bound very conservative. It is also necessary to point out that [71] has assumed a tau lepton reconstruction efficiency as high as 70%, muon and jet detection efficiencies as high as 100% [82]. We also suspect that only the e + e − → qqµ ± τ ∓ν ν processes are considered as possible  background sources in [71] is too optimistic. The potential of the ILC to search for the H → µ ± τ ∓ channel has also been studied in [72], where the Z bosons are reconstructed using lepton pairs.
There it is discussed that a signal with 3σ statistical significance at the ILC with an integrated luminosity of 1 ab −1 requires H → µ ± τ ∓ to have a branching ratio larger than 4.09 × 10 −3 .

Constraints on the SMEFT
We also consider the constraints on the new physics cut-off scale Λ implicated by the improved bounds on the CLFV Higgs decay rates in the SMEFT [84,85], with f i,j the mass eigenstates and i = j. Assuming C ij ∼ 1, the expected CEPC (ILC) constraint on the H → e ± µ ∓ branching ratio will give the most stringent lower bound on Λ, Λ 25 (22) TeV. However, the order of C ij depends crucially on flavor structures beyond the SM. If we adopt the Cheng-Sher ansatz C ij ∼ √ m i m j /v [90], the H → µ ± τ ∓ channel will set the most stringent lower bound on Λ, which reads Λ 0.6 (0.5) TeV.

Constraints on the type III 2HDM
Although one Higgs field is enough for electroweak symmetry breaking and giving masses to gauge bosons and fermions as in the SM, there are several motivations for introduction of two Higgs doublets [91], including requirement of supersymmetry [92], axion models [93] and baryogenesis [94] (see e.g. [95]). In a two-Higgs-doublet-model (2HDM), couplings of the Higgs boson with other particles are modified compared to the SM. Especially, the type III 2HDM naturally introduces tree-level CLFV Higgs couplings. In the type III 2HDM, two doublets Φ 1 = 1/ √ 2(..., v 1 + ρ 1 + ...) T and Φ 2 = 1/ √ 2(..., v 2 + ρ 2 + ...) T with hypercharge +1 couple to fermions freely. We can rotate the scalar doublets such that the vacuum expectation is entirely in the first doublet, with v = v 2 1 + v 2 2 and tan β = v 2 /v 1 . While the mass eigenstates are given by another rotation, and equivalently Diagonalizing the mass matrix automatically diagonalize the H 0 1 coupling matrix, so the CLHV vertices only come from H 0 2 . Therefore, the CLFV couplings with H is given by Under the Cheng-Sher ansatz [90], we define ξ ij = λ ij 2m i m j /v, where λ ij are of order one. On the other hand, the flavor conserving H couplings receive contributions from both H 0 1 and H 0 2 . For example, the expression for the Hbb coupling is given by y b sin(α − β) + ξ b cos(α − β) with According to the current measurements, no obvious deviation of Higgs couplings from the SM has been found [97]. This is further confirmed by the recent observation of the H → bb channel by the CMS [98] and the ATLAS [96] collaborations. It indicates that the parameters λ i(j) and α − β are strictly constrained, which also makes (8) approximately valid. Here we study how the bounds on the CLFV Higgs decay rates together with the Hbb coupling set constraints on the relevant parameters λ and cos(α − β) in the 2HDM. The H → µ ± τ ∓ channel is taken as an example.
From the ATLAS measurement of the H → bb signal strength [96], it can be extracted that the ratio of the Hbb coupling to the SM expectation is 1.005 ± 0.10. We choose the order-one λ b to be 1 and 0.5 as two benchmarks, and display the corresponding 1-and 2-sigma allowed ranges of cos(α − β) in the two panels of FIG. 2 by the orange and light orange contours. In FIG. 2, we also show on the cos(α − β)-λ µτ plain the constraints (the regions above the curves are excluded) set by the upper bounds on B(H → µ ± τ ∓ ). It is observed that cos(α − β) still has a large living space especially in the λ b = 0.5 case, and that in the region when | cos(α − β)| is large, the upper bounds on B(H → µ ± τ ∓ ) given by the the future lepton colliders restrict λ µτ to be smaller than O(0.1).

Constraints on RS models
In RS models [99,100] in which the fermions are allowed to propagate in the extra dimension, the large fermion mass hierarchies and the tiny neutrino masses can be explained [101,102]. In order to generate the observed structure in the lepton sector, which means the hierarchies between charged lepton masses, the neutrino masses with a similar size and the large neutrino mixing angles, one can, in the assumption of Dirac neutrinos, set same or similar profiles in the fifth dimension for the lepton doublets and neutrino singlets and set different profiles for the charged lepton singlets (see e.g. [103]). In such a case, we find that the coefficient in (11) C ij ∼ m τ /v (see Section 4.1 of [104]).
Therefore, in such models, the most stringent lower bound on Λ, or the famous Kaluza-Klein scale M KK [105,106], is set by the expected CEPC (ILC) constraint on the H → µ ± τ ∓ branching ratio as Λ 2.5 (2.2) TeV. Since the masses of the lightest Kaluza-Klein particles are approximately 2.45M KK , nondiscovery of H → µ ± τ ∓ at the CEPC (ILC) will excluded Kaluza-Klein particles with masses smaller than 6.1 (5.4) TeV.

Constraints on models with heavy neutrinos
Lepton flavor violation may originate from heavy neutrinos at one-loop level (see e.g. [107]).
Taking the Inverse Seesaw Model as an example with the right-handed neutrino masses M R close to the TeV scale, we have approximately the off-diagonal charged lepton Yukawa couplings according to (25) of [108], with Y ν the neutrino Yukawa coupling matrix, g the SU (2) gauge coupling constant, m i the mass of the ith generation charged lepton, m W the W boson mass and m H the Higgs boson mass. The function r(λ) is given by (26) of [108]. If we assume a benchmark neutrino Yukawa coupling matrix following [108], Y ν = {(0.1, 0, 0), (0, 1, 0), (0, 1, 0.014)}, a rough calculation indicates that the lower bound on M R set by the expected measurements of the CLFV Higgs decay channels at the three future lepton colliders will be M R 0.3 GeV. Since such a small right-handed neutrino mass would not satisfy the perturbation condition, we conclude that the expected improved bounds on the CLFV Higgs decay rates put no constraint on the right-handed neutrino masses. The complete expression for the off-diagonal charged lepton Yukawa couplings, which does not rely on the expansion in inverse powers of M R , can be found in Appendix C of [109]. As discussed in [108], using (16) always overestimates values of Y ij , so using the complete formula will give an even looser constraint on M R , which will not change our main conclusion.

IV. SUMMARY
The future e + e − colliders, the CEPC, the FCC-ee and the ILC, as Higgs factories, are ideal machines for precise studies of Higgs properties. In this paper, we evaluate the potential of the three lepton colliders for searching for the CLFV Higgs decays. We find that the expected upper bounds given by the CEPC or the FCC-ee (the ILC) on the branching ratios of H → e ± µ ∓ , e ± τ ∓ and µ ± τ ∓ are 1.2 (2.1) × 10 −5 , 1.6 (2.4) × 10 −4 and 1.4 (2.3) × 10 −4 at 95% CL, respectively.
The resulting constraints on certain theory parameters are also given, including the CLFV Higgs couplings, the relevant parameters in the type-III 2HDM, and the cut-off scales in the SMEFT and in RS models.