Novel probes for electron-muon flavor violation from exotic Higgs decays

TL;DR

This paper introduces new collider signatures involving Higgs decays to a light pseudoscalar that then decays into an electron-muon pair, offering a promising way to detect lepton flavor violation.

Contribution

It proposes novel collider signatures for electron-muon flavor violation via exotic Higgs decays involving a light pseudoscalar in the type-III Two-Higgs-doublet-model.

Findings

Proposed signatures lead to stronger constraints on lepton flavor violating couplings than current low-energy measurements.

Analyzed parameter space where the light pseudoscalar is consistent with existing constraints.

Collider searches can complement low-energy experiments in probing new physics.

Abstract

In this paper, we propose two novel signatures of Higgs decays to search for electron-muon flavor violation. These signatures arise from the presence of a light pseudoscalar into which the 125-GeV Higgs boson decays. The pseudoscalar subsequently decays into an electron-muon pair, leading to multilepton final states, which are relatively clean signatures to search for at the LHC. As a benchmark, we consider the type-III Two-Higgs-doublet-model. We analyze both low-energy and collider constraints on the model and identify regions of parameter space where the light pseudoscalar is viable. Our proposed signatures yield stronger constraints on the lepton flavor violating couplings than current low-energy precision measurements. Taken together, our findings suggest that collider-based probes of exotic Higgs decays provide a powerful complement to precision experiments in the quest to uncover new physics.

Tables

Table 1. Table 1 : Numerical values for f ( q , N ) f^{(q,N)} . Note that Q Q stands for heavy quarks, c , b , t c,b,t .

Nucleon	$f^{(u,N)}$ Bishara et al. (2016)	$f^{(d,N)}$ Bishara et al. (2016)	$f^{(s,N)}$ Junnarkar and Walker-Loud (2013)	$f^{(Q,N)}$ Harnik et al. (2013)
$p$	0.018$\pm$0.005	0.034$\pm$0.011	0.043$\pm$0.011	0.067$\pm$0.001
$n$	0.016$\pm$0.005	0.038$\pm$0.011	0.043$\pm$0.011	0.067$\pm$0.001

Equations

$$V=$$

$$+\frac{1}{2}{\lambda }_1({\Phi }_1^{†}{\Phi }_1^{}{)}^2+\frac{1}{2}{\lambda }_2({\Phi }_2^{†}{\Phi }_2^{}{)}^2+{\lambda }_3({\Phi }_1^{†}{\Phi }_1^{})({\Phi }_2^{†}{\Phi }_2^{})+{\lambda }_4({\Phi }_1^{†}{\Phi }_2^{})({\Phi }_2^{†}{\Phi }_1^{})$$

$$+\left\{\frac{1}{2}{\lambda }_5({\Phi }_1^{†}{\Phi }_2^{}{)}^2+\left[{\lambda }_6({\Phi }_1^{†}{\Phi }_1^{})+{\lambda }_7({\Phi }_2^{†}{\Phi }_2^{})\right]({\Phi }_1^{†}{\Phi }_2^{})+\text{h.c.}\right\}.$$

$${\Phi }_1=\left(\begin{array}{c}{G}^{+}\\ \frac{v+{\varphi }_1+i{G}^0}{\sqrt{2}}\end{array}\right),{\Phi }_2=\left(\begin{array}{c}{H}^{+}\\ \frac{{\varphi }_2+iA}{\sqrt{2}}\end{array}\right),$$

$$M_{11}^2=-\frac{1}{2}{\lambda }_1{v}^2,M_{12}^2=\frac{1}{2}{\lambda }_6{v}^2.$$

$$m_{H^{\pm }}^2$$

$$m_A^2$$

$${\mathcal{M}}^2=\left(\begin{array}{cc}{\lambda }_1{v}^2& {\lambda }_6{v}^2\\ {\lambda }_6{v}^2& m_{H^{\pm }}^2+\frac{1}{2}({\lambda }_4+{\lambda }_5)v^2\end{array}\right).$$

$$\left(\begin{array}{c}{\varphi }_1\\ {\varphi }_2\end{array}\right)=\left(\begin{array}{cc}{c}_{\alpha }& s_{\alpha }\\ -s_{\alpha }& c_{\alpha }\end{array}\right)\left(\begin{array}{c}h\\ H\end{array}\right),$$

$$s_{2\alpha }=\frac{2{\lambda }_6{v}^2}{m_H^2-m_h^2},$$

$$m_{h,H}^2=\frac{1}{4}\left(2m_{H^{\pm }}^2+(2{\lambda }_1+{\lambda }_4+{\lambda }_5)v^2\mp \sqrt{{\left(2m_{H^{\pm }}^2+(-2{\lambda }_1+{\lambda }_4+{\lambda }_5)v^2\right)}^2+16{\lambda }_6^2{v}^4}\right).$$

$$m_h^2$$

$$m_H^2$$

$${\lambda }_1,{\lambda }_2$$

$${\lambda }_3$$

$${\lambda }_3+{\lambda }_4-|{\lambda }_5|$$

$${\lambda }_1+{\lambda }_2+2({\lambda }_3+{\lambda }_4+{\lambda }_5)$$

$$3({\lambda }_1+{\lambda }_2)+\sqrt{9({\lambda }_1-{\lambda }_2{)}^2+(2{\lambda }_3+{\lambda }_4{)}^2}$$

$${\lambda }_1+{\lambda }_2+\sqrt{({\lambda }_1-{\lambda }_2{)}^2+4{\lambda }_4^2}$$

$${\lambda }_1+{\lambda }_2+\sqrt{({\lambda }_1-{\lambda }_2{)}^2+4{\lambda }_5^2}$$

$$|{\lambda }_3+2{\lambda }_4\pm 3{\lambda }_5|$$

$$|{\lambda }_3\pm {\lambda }_4|$$

$$|{\lambda }_3\pm {\lambda }_5|$$

$$\Delta T$$

$$F(x,y)=\frac{x+y}{2}-\frac{xy}{x-y}\ln \frac{x}{y}.$$

$${\mathcal{L}}_{Yuk}\supset -\frac{\sqrt{2}{m}_{\ell }^i}{v}{\delta }^{ij}{\bar{L}}_i{\ell }_{Rj}{\Phi }_1-\frac{\sqrt{2}{m}_U^i}{v}{\delta }^{ij}{\bar{Q}}_i{u}_{Rj}{\tilde{\Phi }}_1-\frac{\sqrt{2}{m}_D^k}{v}{V}^{ik}{\delta }^{kj}{\bar{Q}}_i{d}_{Rj}{\Phi }_1+\text{h.c.},$$

$$L=\left(\begin{array}{c}{\nu }_L\\ {\ell }_L\end{array}\right),Q=\left(\begin{array}{c}{u}_L\\ V{d}_L\end{array}\right),$$

$${\mathcal{L}}_{Yuk}\supset -\sqrt{2}{Y}_{ij}{\bar{L}}_i{\ell }_{Rj}{\Phi }_2+\text{h.c.}.$$

$${\mathcal{L}}_{eff}=$$

$$+(L\leftrightarrow R)+\text{h.c.},$$

$${\Gamma }_{\mu \to e\gamma }=\frac{{\alpha }_{em}{m}_{\mu }^5}{256{\pi }^4{v}^4}\left(|c_L{|}^2+|c_R{|}^2\right).$$

$$c_L^{(1)}=$$

$$c_R^{(1)}=$$

$$c_L^{(2t)}=\frac{s_{2\alpha }}{3}\frac{{\alpha }_{em}{Y}_{\mu e}}{\pi }\frac{m_t}{m_{\mu }}\left[\frac{m_{t}v}{m_h^2}f\left(\frac{m_t^2}{m_h^2}\right)-\frac{m_{t}v}{m_H^2}f\left(\frac{m_t^2}{m_H^2}\right)\right],$$

$$f(z)$$

$$⟨N|\bar{q}{\gamma }_{\rho }q|N⟩=n^{(q,N)}\bar{N}{\gamma }_{\rho }N,$$

$$⟨N|\bar{q}q|N⟩=f^{(q,N)}\frac{m_N}{m_q}\bar{N}N,$$

$$\Gamma (\mu \to e\text{ conv.})=m_{\mu }^5{\left|\frac{e{c}_{R}D}{32v^2{\pi }^2}+(2g_{LV}^u+g_{LV}^d)V^p+\sum _q\sum _{N=p,n}{g}_{LS}^q\frac{m_N}{m_q}{f}^{(q,N)}{S}^N\right|}^2+(L\leftrightarrow R),$$

Full text

1 Introduction

After a decade since the discovery of the Higgs boson, the era of precision Higgs measurements has opened up. The primary goal of these precise measurements is to determine if the Higgs sector strictly adheres to the Standard Model (SM) predictions, or harbors the sign of new physics. Although current measurements are in good agreement with the SM Aad and others (2016); Tumasyan and others (2022b); Aad and others (2022a), there are still room for small deviations that could point to extensions of the SM. Two major avenues are being pursued to probe such deviations: first, the high precision measurement of the SM-like Higgs couplings to the SM particles; and second, the search of exotic Higgs decays that are not predicted in the SM.

Among the exotic decays, the lepton-flavor violation (LFV) varieties offer particularly sensitive windows to the new physics. The SM strictly forbids LFV, so any observations of such a decay would be a clear indication of physics beyond the SM. Motivated by this, ATLAS and CMS collaborations have searched for the two-body final states LFV decays, i.e., $h\to\tau^{\pm}\mu^{\mp}$ and $h\to\tau^{\pm}e^{\mp}$ Sirunyan and others (2021); Aad and others (2023), as well as $h\to\mu^{\pm}e^{\mp}$ Aad and others (2020); Hayrapetyan and others (2023). These searches so far have yielded null results, and hence they place the upper limits on the LFV branching ratios of the Higgs boson of the order 0.1% for the first two channels and $\mathcal{O}\left(10^{-5}\right)$ for the latter mode.

These collider bounds must be interpreted in light of stringent constraints from the low-energy LFV observables. In the $e\mu$ case, for instance, the low-energy LFV processes $\mu\to e\gamma$ and $\mu\to e$ conversions in nuclei provide extremely tight limits, which are few orders of magnitude stronger than the bounds from $h\to\mu^{\pm}e^{\mp}$. Typically, the bounds are derived in the effective field theory (EFT) framework under the assumption that no new degrees of freedom exist around the electroweak scale Blankenburg et al. (2012); Harnik et al. (2013). However, such an assumption is invalid in general. In fact, the presence of new degrees of freedom around the electroweak scale can significantly affect both the low-energy and the high-energy observables Buschmann et al. (2016); Primulando et al. (2023); Koivunen and Raidal (2023); Afik et al. (2024); Crivellin et al. (2013). Thus, a consistent interpretation of LFV signals requires embedding the EFT into an ultraviolet (UV)-complete model.

In this work, we explore the LFV Higgs decays within the framework of the type-III Two-Higgs-Doublet-Model (2HDM), which extends the scalar sector by the introduction of the second Higgs doublet without imposing any discrete symmetries for flavor conservations. This allows for a tree-level LFV Higgs couplings. Unlike effective theories, the 2HDM provides a concrete UV structure linking LFV observables across energy scales. Importantly, it opens up additional decay channels for the SM-like Higgs bosons through its interactions with other scalar states.

We consider a specific scenario in which the pseudoscalar $A$ is lighter than the SM-like Higgs boson $h$. In this setup, the $h$ can decay via $h\to AZ$ and $h\to AA$, with the pseudoscalar subsequently decaying through the LFV channel $A\to\mu^{\pm}e^{\mp}$. These processes lead to distinctive multilepton signatures, such as $h\to\mu^{\pm}e^{\mp}\ell^{+}\ell^{-}$ or $h\to\mu^{\pm}e^{\mp}\mu^{\pm}e^{\mp}$, with little background from SM processes. Although no LHC searches have yet targeted these exact final states, existing analyses in the related multilepton or exotic Higgs channels can be repurposed to constrain them.

This paper presents a systematic study of these novel LFV signatures at the LHC. In section 2, we provide a brief summary of the type-III 2HDM. We then discuss its low-energy LFV observables in section 3. In section 4, we discuss the relevant collider constraints for the scenario where the pseudoscalar is lighter than the SM-like Higgs boson. We then propose novel collider signatures for the SM-like Higgs decays in section 5. We then conclude and discuss our results in section 6.

2 Type-III 2HDM

In this section, we give a brief overview of the Type-III 2HDM. We will closely follow the notation and convention of ref. Primulando and Uttayarat (2017). The two electroweak scalar doublets are denoted by $\Phi_{1}$ and $\Phi_{2}$. Their scalar potential is given by

[TABLE]

In principle, the parameters $M_{12}^{2}$, $\lambda_{5}$, $\lambda_{6}$ and $\lambda_{7}$ can be complex, which would lead to $CP$-violation in the scalar sector. In this work, for simplicity, we will assume that such parameters are real, so that the scalar sector is $CP$ symmetric. For the present discussion, it will be more convenient to work in the Higgs basis Georgi and Nanopoulos (1979), that is, the doublets $\Phi_{1}$ and $\Phi_{2}$ are expanded as

[TABLE]

where $v=246$ GeV is the electroweak vacuum expectation value, $G^{+}$ and $G^{0}$ are the would-be Goldstone bosons, $H^{+}$ is the physical charged Higgs boson, $A$ is the physical $CP$-odd Higgs boson and $\phi_{1,2}$ are two $CP$-even neutral scalars.

From the minimization of the scalar potential, one can express some parameters in terms of the others. For instance,

[TABLE]

From here, one can derive the masses of $H^{\pm}$ and $A$

[TABLE]

and the mass matrix of $\phi_{1}$ and $\phi_{2}$ in the basis of $(\phi_{1},\phi_{2})$

[TABLE]

The above mass matrix is diagonalized by rotating the basis into $h$ and $H$ mass eigenbasis

[TABLE]

where $c_{x},s_{x}$ denote $\cos x,\sin x$. The mixing angle is given by

[TABLE]

resulting in eigenvalues

[TABLE]

In this paper, we identify $h$ with the 125-GeV Higgs boson and $H$ with the heavy $CP$-even Higgs boson. The bound on the mixing angle $\alpha$ is determined from the combined Run 1 and Run 2 measurements of the ATLAS and CMS, which yields $s_{\alpha}\leq 0.15$ at 95% confidence level (CL) Aad and others (2016); Tumasyan and others (2022b); Aad and others (2022a). In the small $s_{\alpha}$ limit, the masses of $h$ and $H$ are approximated as

[TABLE]

The presence of $M_{22}^{2}$ makes it possible for $H$, $A$ and $H^{\pm}$ masses to be much greater than the weak scale. However, their mass-squared differences are controlled by quartic couplings, $\lambda_{i}$, theoretically constrained by perturbativity, vacuum stability and unitarity conditions. In our analysis, we take $|\lambda_{i}|<4\pi$ for perturbativity. The vacuum stability and unitarity constraints, in the presence of $\lambda_{6}$ and $\lambda_{7}$, are complicated and not illuminating. The coupling $\lambda_{6}$ is directly proportional to the mixing angle $s_{\alpha}$, so it is expected to be small. The coupling $\lambda_{7}$, on the other hand, does not play any role in our analysis. Hence, for simplicity, we will take $\lambda_{7}=0$. With these simplifications, the vacuum stability constraints read Ferreira et al. (2004); Ivanov (2007); Ferreira and Jones (2009)111The constraint involving $\lambda_{6}$ is only the necessary condition for the scalar potential to be bounded below.

[TABLE]

and the unitarity constraints are given by Ginzburg and Ivanov (2003)

[TABLE]

The scalar mass splittings are also constrained by the $T$ parameter, inferred from the electroweak precision measurements. The new physics correction to the $T$ parameter reads

[TABLE]

where $\alpha_{em}^{-1}=137$ is the electromagnetic fine-structure constant in the Thompson limit, whereas $F(x,y)$ is defined as

[TABLE]

One can notice that if one of the neutral scalar masses is degenerate with the charged one, $\Delta T$ will vanish. This is consistent with the value of $\Delta T$ inferred from the fit to the electroweak precision data, which is $\Delta T=0.00\pm 0.06$ Navas and others (2024).

In this work, we shall consider the case with $m_{A}\lesssim 100~\text{GeV}$ with heavier $m_{H}$ and $m_{H^{\pm}}$. To comply with the $T$ parameter constraint, it is instructive to assume $m_{H}=m_{H^{\pm}}$, which in turn implies $\lambda_{4}+\lambda_{5}=0$, see eq. (11). The value of $m_{H^{\pm}}$ is constrained by LHC searches. For leptophilic charged scalar decaying 100% into electron or muon, it is constrained to be heavier than 550 GeV. As a consequence, by virtue of equations. (5) and (14), we expect large quartic couplings. We will discuss it in detail in section 5.

The Higgs doublet $\Phi_{1}$ is responsible to generate fermion masses thorugh Yukawa interactions

[TABLE]

where $m_{\ell}$, $m_{U}$ and $m_{D}$ are the diagonal lepton, up-type quark and down-type quark mass matrices respectively, $V$ is the Cabibbo-Kobayashi-Maskawa matrix and $\tilde{\Phi}_{a}\equiv i\sigma_{2}\Phi^{\ast}_{a}$. In the above equation, the left-handed fermion doublets are taken to be

[TABLE]

where fields $\ell_{L(R)}$, $u_{L(R)}$ and $d_{L(R)}$ are defined in their respected mass eigenbases. The Yukawa couplings of $\Phi_{2}$ in general lead to flavor violation. In this work, we will focus on the leptonic couplings

[TABLE]

The coupling $Y_{ij}$ in principle can be complex, but for simplicity, we assume that they are real. In our analysis, we mainly focus on the $Y_{e\mu}$ and $Y_{\mu e}$ couplings. However, other components of the $Y$ matrix, e.g., $Y_{\tau\tau}$, can also be present. This observation will play an important role when considering the collider searches in multilepton channels to be discussed in section 5.

3 Low-energy lepton flavor-violating constraints

The Yukawa couplings in eq. (26), together with the neutral scalar mixing angle $s_{\alpha}$, lead to LFV decays of charged lepton. In particular, nonzero $Y_{e\mu}$ and $Y_{\mu e}$ give rise to $\mu\to e\gamma$ and $\mu\to e$ conversion in atomic nuclei. These low-energy LFV processes can be described by effective operators

[TABLE]

where $P_{R,L}=(1\pm\gamma_{5})/2$ are the right- and left-chirality projection operators, respectively. Terms with coefficients $c_{L,R}$ are dipole operators, which are relevant for both $\mu\to e\gamma$ and $\mu\to e$ conversion processes. The last two terms are effective vector and scalar current interactions between leptons and quarks which are relevant for $\mu\to e$ conversion.

The partial decay width for $\mu\to e\gamma$ is given by

[TABLE]

The coefficients $c_{L,R}$ are induced by quantum effects and are expressed as $c_{L,R}=c_{L,R}^{(1)}+c_{L,R}^{(2)}$, with the superscript denoting the loop level. The one-loop contributions are induced through Feynman diagram in figure 1 (left) and are given by Hisano et al. (1996)

[TABLE]

It has long been realized in the literature that numerically large contributions to the dipole operators can arise from two-loop Feynman diagrams due to possible large hierarchies among the couplings Weinberg (1989); Dicus (1990); Barr and Zee (1990). For example, the so-called Barr-Zee diagram with an internal top-quark loop, shown in the middle of figure 1, results in Chang et al. (1993); Davidson and Grenier (2010); Crivellin et al. (2014)

[TABLE]

where

[TABLE]

For $m_{H}$ below a TeV scale, the loop function $f(z)\lesssim 4\pi$. Thus, the one-loop coefficient is parametrically suppressed by $\,m_{\mu}^{2}/(\alpha_{em}m_{t}^{2})$ compared to the two-loop one. A more complete analysis of two-loop contributions to the dipole coefficients has been carried out in ref. Altmannshofer et al. (2025), including the contributions from the scalar quartic couplings, see figure 1 (right). Such contributions could be significant, particularly in our scenario where we have large scalar mass splittings. We use the Python code provided in ref. Altmannshofer et al. (2025) to compute the coefficients $c_{L}$ and $c_{R}$ numerically.

Experimentally, the most stringent constraint on $\mu\to e\gamma$ decay is provided by the MEG II experiment, with BR$(\mu\to e\gamma)\leq 1.5\times 10^{-13}$ at 90% CL Afanaciev and others (2025). It should be noted that the MEG II experiment is still running, collecting more data. It is projected that by 2026, the MEG II experiment should be able to push the bound down to BR$(\mu\to e\gamma)\leq 6.0\times 10^{-14}$ at 90% CL.

The $\mu\to e$ conversion rate in atomic nuclei is computed by matching the quark-level effective operators in eq. (27) onto the nucleon-level effective operators. The matching is done by evaluating the nucleon matrix elements

[TABLE]

where $n^{(q,N)}$ is the number of valence quark $q$ inside the nucleon $N$ and $f^{(q,N)}$ is the form factor. For light quarks, $q=u,d,s$, the form factors are determined through lattice calculations. The heavy quarks ($q=c,b,t$) form factors are related to those of the light quarks through the trace anomaly Shifman et al. (1978); Jungman et al. (1996). The values of quark form factors are listed in table 1.

Equipped with the quark form factors, one can readily compute the $\mu\to e$ conversion rate. It is given by

[TABLE]

where $D$, $V^{p}$ and $S^{p,n}$ are the overlap integrals. Their numerical values, for various atomic nuclei, have been tabulated in ref. Kitano et al. (2002). The effective couplings $g_{LV}^{q}$ and $g_{LS}^{q}$ are given by

[TABLE]

The effective coupling $g_{RV}^{q}$ and $g_{RS}^{q}$ can be obtained from their left-handed counterparts by the replacement $Y_{e\mu}\to Y_{\mu e}$. Note that we only include photon-penguin diagrams in vector couplings and tree-level scalar exchange in scalar ones. In principle, there are other contributions from $Z$-penguin and box diagrams, but they are deemed insignificant. The $Z$-exchange diagrams, for example, are suppressed by $m_{\mu}^{2}/m_{Z}^{2}$, while the box diagrams are suppressed by light quark masses. In addition, the box contributions can also result in scalar operators. However, that requires a chirality flip from muon propagator, inducing an additional $m_{\mu}/v$ suppression.

The experimental search $\mu\to e$ conversion in gold nuclei provides the strongest constraint. The SINDRUM II collaboration has placed an upper limit $\Gamma(\mu\to e$ conv.)/$\Gamma$(captured) $<7\times 10^{-13}$ at 90% CL Bertl and others (2006). The corresponding muon capture rate in the gold nucleus is $13.07\times 10^{6}\text{ s}^{-1}$ Suzuki et al. (1987). Finally, the overlap integrals for the gold nucleus are given by $D=0.189$, $V^{p}=0.0974$, $S^{p}=0.0614$ and $S^{n}=0.0918$.

As has already been mentioned, the $\mu\to e\gamma$ and $\mu\to e$ conversion constraints depend on the quartic couplings $\lambda_{3}$, $\lambda_{4}$ and $\lambda_{5}$ through the dipole coefficients $c_{L}$ and $c_{R}$. The couplings $\lambda_{4}$ and $\lambda_{5}$ are determined from the scalar masses $m_{A}$, $m_{H}$ and $m_{H^{\pm}}$, leaving $\lambda_{3}$ as an independent parameter. However, in the event that $h\to AA$ is open, the ATLAS $h\to AA$ search dictates $\lambda_{3}\simeq 2(m_{H^{\pm}}^{2}-m_{A}^{2})/v^{2}$, see section 5.2.

In a scenario where $m_{A}\leq 100$ GeV and $m_{H}=m_{H^{\pm}}\gtrsim 300$ GeV, we find that the $\mu\to e\gamma$ and $\mu\to e$ conversion constraints depend only on whether the decay channel $h\to AA$ is open or closed, and not on the actual value of $m_{A}$. In this region of parameter space, $\mu\to e\gamma$ provides a tighter constraint than $\mu\to e$ conversion. Figure 2 shows the conservative bounds on $s_{2\alpha}\sqrt{Y_{e\mu}^{2}+Y_{\mu e}^{2}}$ as a function of $m_{H}$ for the cases where the decay channel $h\to AA$ is open (solid lines) and closed (dashed lines). In obtaining these bounds in the scenario where $h\to AA$ is closed, the coupling $\lambda_{3}$ is taken to be at the minimum value allowed by Eq. (14), so that the constraint from $\mu\to e\gamma$ is at its weakest. Allowing $\lambda_{3}$ to take on a larger value improves the $\mu\to e\gamma$ limit by at most 25% (4%) for $m_{H}$ = 300 (600) GeV.

4 Current collider constraints on the model

In our setup, the doublet $\Phi_{2}$ couples to leptons but not quarks, see eq. (26). As a result, only the heavy scalar $H$ can be singly produced, provided that it mixes with the SM-like Higgs boson. In our analysis, we will focus on the LFV decays induced by coupling $Y_{e\mu}$ and $Y_{\mu e}$. The CMS collaboration has conducted a dedicated search for the heavy Higgs boson decaying into an electron-muon pair in the mass range 110 GeV $<m_{H}<$ 160 GeV Hayrapetyan and others (2023). This search provides the most stringent bound on the LFV decay of $H$. For heavier $m_{H}$, the LFV constraints can be obtained from searches for a resonant production of the sneutrino decaying to $e\mu$. The mass window of 160–200 GeV is covered by the D0 search Abazov and others (2010). The corresponding CMS search covers the mass from 200 GeV up to several TeV Tumasyan and others (2023). We reinterpret these analyses by assuming that $H$ and the sneutrino have comparable signal acceptance. Since the $H$ production cross-section is proportional to the mixing angle squared, the resulting limits are presented in terms of $s^{2}_{\alpha}\times$ BR$(H\rightarrow e^{\pm}\mu^{\mp})$, see figure 7(b). Note that the LFV constraints of the D0 search are weak, with $s^{2}_{\alpha}\times$ BR$(H\rightarrow e^{\pm}\mu^{\mp})<\mathcal{O}\left(10^{-2}\right)$. We further note that for $m_{H}>m_{A}+m_{Z}$, the decay $H\rightarrow AZ$ dominates due to the smallness of the LFV Yukawas, suppressing the branching ratio BR$(H\rightarrow e^{\pm}\mu^{\mp})$.

Additional constraints arise from the pair production of $H$ and $A$ via an $s$-channel $Z$ boson, followed by $H\to AZ$ and $A\to e^{\pm}\mu^{\mp}$. The resulting final state contains four leptons, coming from the $AA$ pair, and the $Z$ boson. The $Z$ boson will then decay to a pair of leptons, quarks or missing energy. This process is constrained by the CMS multilepton analysis Tumasyan and others (2022c), in particular the $4\ell$G and $4\ell$H signal regions. Both signal regions require at least four leptons with $p_{T}>10$ GeV. If more leptons are present, the four highest $p_{T}$ leptons are selected. In the $4\ell$G signal region, the four selected leptons are grouped into pairs of opposite-sign same-flavor (OSSF) leptons. The two pairs are indicated by $Z_{1}$ and $Z_{2}$, where the invariant mass of $Z_{1}$ is the closest to the mass of the $Z$ boson. The region is further categorized according to the invariant mass of $Z_{2}$, and the transverse mass variable ($M_{T2}$) that involves $Z_{1}$, $Z_{2}$ and the missing transverse energy. Meanwhile, the signal region $4\ell$H requires exactly one OSSF pair, denote $Z_{1}$. The signal region is characterized by the invariant mass of $Z_{1}$ and the separation $\Delta R$ between the non-OSSF leptons pair.