New Physics Models Facing Lepton Flavor Violating Higgs Decays at the Percent Level

We speculate about the possible interpretations of the recently observed excess in the $h \to \tau \mu$ decay. We derive a robust lower bound on the Higgs boson coupling strength to a tau and a muon, even in presence of the most general new physics affecting other Higgs properties. Then we reevaluate complementary indirect constraints coming from low energy observables as well as from theoretical considerations. In particular, the tentative signal should lead to $\tau \to \mu \gamma$ at rates which could be observed at Belle II. In turn we show that, barring fine-tuned cancellations, the effect can only be accommodated within models with an extended scalar sector. These general conclusions are demonstrated using a number of explicit new physics models. Finally we show how, given the $h \to \tau \mu$ signal, the current and future searches for $\mu \to e \gamma$ and $\mu \to e$ nuclear conversions unambiguously constrain the allowed rates for $h \to \tau e$.


I. INTRODUCTION
The discovery of the Higgs boson at the LHC [1,2] imbues the standard model (SM) of particle physics with completeness and self-consistency. Nonetheless, its failure to account for nonvanishing neutrino masses is one of the main motivations for considering physics beyond the SM. Incidentally, the accidental SM symmetries that prevent neutrinos from acquiring mass also completely suppress lepton flavor violating (LFV) processes. The observation of the former thus provides ample motivation for a rich experimental program to search for the latter. The CMS collaboration has recently reported a slight excess with a significance of 2.4 σ in the search for LFV decay h → τ µ [3]. The best fit for the branching ratio of the Higgs boson to τ µ final state (summed over τ − µ + and τ + µ − ), assuming SM Higgs production, is found to be B(h → τ µ) = 0.84 +0. 39 −0.37 % .
It is thus an imperative to either confirm or reject validity of this tantalizing hint with more data by both ATLAS and CMS experiments. At the same time, it is instructive to revisit expectations for this observable within various new physics (NP) scenarios and in particular re-evaluate the feasibility of obtaining such a large signal in light of severe indirect constraints on LFV Higgs interactions coming from low energy probes.
Without loss of generality, one can parameterize the mass terms and Higgs boson couplings of charged leptons after electroweak symmetry breaking (EWSB) as where the ellipsis denotes non-renormalizable interactions involving more than one Higgs boson and i = e, µ, τ . In the SM, the Higgs couplings are diagonal and given by y ij = (m i /v)δ ij , where v = 246 GeV . On the other hand, non-zero y τ µ and/or y µτ induce h → τ µ decays with a branching ratio of Assuming that the total Higgs boson decay width (Γ h ) is given by its SM value enlarged only by the contribution from h → τ µ itself, i.e., Γ h = Γ SM h /[1 − B(h → τ µ)], where Γ SM h (m h = 125 GeV) = 4.07 MeV [11], the measurement in Eq. (1) can be interpreted as a two-sided bound on the |y τ µ | 2 +|y µτ | 2 combination of couplings. These limits read (see also the left-hand side panel of Fig. 1) 0.0019(0.0008) < |y τ µ | 2 + |y µτ | 2 < 0.0032(0.0036) at 68% (95%) C.L. . (4) In general, the experimentally measured h → τ µ event yield depends not only on the values of y τ µ and y µτ , but also on other Higgs couplings contributing both to its total decay width Γ h as well as its production cross-section (σ h ). In particular, a given signal can be reproduced for larger (smaller) values of |y τ µ | and |y µτ | by enhancing (suppressing) Γ h and/or suppressing (enhancing) σ h . Since both Γ h and σ h affect other currently measured Higgs observables, their individual effects can be disentangled by performing a global fit to all Higgs production and decay event yields at the LHC. The details of this procedure can be found in Appendix A, while the resulting best fit χ 2 values as functions of |y τ µ | 2 + |y µτ | 2 are shown in the right-hand side panel of Fig. 1.
In particular, σ h is well determined by the measurements of both inclusive and separate exclusive Higgs production channels in several different decay modes. Similarly, Γ h is bounded from below by the observations (or at least indications [12]) of the dominant SM Higgs decay modes (h → bb, h → W W * , h → τ + τ − , etc.). Consequently the lower range of allowed |y τ µ | and |y µτ | values does not change much compared to Eq. (4) when considering the global fit. On the other hand, the fact that the total Higgs decay width is currently only weakly bounded from above [13,14] allows significantly larger |y τ µ | and |y µτ | couplings to reproduce the same observed signal in the general case. 1 Numerically, we find 0.0017(0.0007) < |y τ µ | 2 + |y µτ | 2 < 0.0036(0.0047) at 68% (95%) C.L. .
The rest of the paper is devoted to interpreting these ranges in terms of hypothetical new physics (NP) effects. In Sec. II we review the model independent considerations of flavor violating Higgs couplings using an effective field theory approach. Sec. III is devoted to a case study of the two Higgs doublet model (THDM) with a generic Yukawa structure (the so-called type-III THDM) that can account for the observed anomaly at tree level and be in agreement with present low energy flavor constraints. Scenario that relies on vector-like fermions to explain the data is presented in Sec. IV. Possibility to explain the anomaly via one-loop physics effects is discussed in Sec. V, where we investigate phenomenology of scalar leptoquarks to demonstrate difficulties with this 1 It has been argued recently [15,16], that the current total Higgs decay width measurements using off-peak data [13,14] introduce some model-specific assumptions into the fit. We have checked that removing the total Higgs width constraint from the fit, no upper bound on |y τ µ | 2 + |y µτ | 2 can be set.  type of approach. In particular, Sec. V D gives a very special solution that can accommodate the anomaly through an interplay between two different one-loop effects, and satisfy constraints imposed by low energy flavor physics experiments at the price of fine-tuning. We briefly conclude in Sec. VI.

II. MODEL INDEPENDENT CONSIDERATIONS
In the following we assume the SM contains all the relevant degrees of freedom at energies O(few 100) GeV. 2 In particular, any additional heavy field can be integrated out, so that we can apply effective field theory (EFT) methods. If the electroweak symmetry is realized linearly, then current data already strongly suggest that h for the most part should be a neutral CP even component of a linear combination of Y = 1/2 hypercharge weak doublets H α = (h + α , v α +x α h+ . . .) T , where the ellipsis denotes possible projections to other heavier scalar mass eigenstates and α |x α | 2 < ∼ 1/2. In general, EWSB can receive contributions from several sources (for example several Higgs multiplet condensates). However, electroweak precision tests severely constrain contributions from higher representations, so that also α v 2 In the minimal extension of the SM with Dirac neutrinos accounting for the observed neutrino masses, a nonvanishing h → τ µ rate is induced at the one-loop level. However, this contribution is negligibly small since the relevant decay amplitude is suppressed by two powers of neutrino Yukawa couplings in addition to the tau Yukawa.
In the decoupling limit we are considering, the contributions from the heavier Higgses are power suppressed and can be completely described in terms of higher dimensional operators. The leptonic Yukawa sector of such a theory, including the leading (dimension six) non-renormalizable operators suppressed by the EFT cut-off scale Λ, can be written as where L i and E i correspond to the leptonic weak doublets (L i = ( i L , ν i ) T ) and singlets (E i = i R ), respectively, and we explicitly show the flavor indices i, j = 1, 2, 3. Any additional dimension six operators coupling the Higgses to the leptons can be shown to be either redundant or not to contribute to the effective Lagrangian (2) after EWSB [17]. We can thus identify where withv α = v α /v and the unitary matrices V L,R diagonalize the leptonic mass term as We can note immediately the two possible sources of non-vanishing y τ µ and/or y µτ . In a theory with multiple Higgs doublets, the first term in the square brackets in Eq. (8) We observe that percent level B(h → τ µ) can be accommodated within the EFT approach (containing a single light Higgs doublet) with a sizable mass gap. At the same time, new degrees Since λ and λ contribute to both as well as m the observed hierarchical structure of the charged lepton masses can be used to define the natural ranges of ij . In particular, without delicate cancellations between the various contributions in Eq. (9), the observed hierarchy between the muon and tau lepton masses implies [18,19] This constraint is compared to the CMS preferred range for |y τ µ | and |y µτ | from Eq. (4) in Fig. 2.
We observe that the two indications are not in sharp conflict and thus the observation of B(h → τ µ) at the percent level can in principle be explained by natural NP although a mild hierarchy between y τ µ and y µτ is preferred in this case.
Another class of model-independent constraints comes from the fact that the dimension six operators contained in L Y in Eq. (6) will in general mix with other operators contained in the effective SM. In particular, restricting the discussion to a single Higgs setup with H ≡ H 1 , λ will mix under charged lepton Yukawa renormalization into (H † H) 3 affecting, for example, Higgs boson pair production at the LHC. However, as shown in Refs. [20,21], such mixing is suppressed by three powers of the charged lepton Yukawa matrix and is thus (baring fine-tuned cancellations) completely negligible in practice. More phenomenologically relevant are the UV finite one-and two-loop contributions to the operatorsLH(σ · B)E andLτ a H(σ · W a )E, where σ µν = i[γ µ , γ ν ]/2, B µν and W a µν are the hypercharge and weak isospin field strengths, respectively, and τ a are the Pauli matrices. These operators can mediate radiative flavor violating lepton decays as well as contribute to leptonic anomalous electric and magnetic dipole moments, all of which are already tightly constrained by experiment. The most stringent constraint comes from the τ → µγ decay [17], mediated by the effective Lagrangian where Q L,Rγ = (e/8π 2 )m τ (μσ αβ P L,R τ )F αβ , P L,R = (1 ∓ γ 5 )/2 and F αβ is the electromagnetic field strength tensor. The Wilson coefficients c L,R receive comparable one-and two-loop contributions. In the experimentally justified approximation y µµ y τ τ (both assumed real) and m µ m τ m h , they are given by [17,22,23] c where y tt is the top quark Yukawa with the SM value of y tt =m t /v = 0.67 for a MS top mass of given by is shown in Fig. 3 (diagonal dashed orange line), assuming SM values of all Higgs boson couplings except y τ µ and y µτ . In the same plot, the CMS preferred range of B(h → τ µ) in Eq. (1) is displayed by the horizontal blue band, while the current (B(τ → µγ) < 4.4 × 10 −8 @ 90% C.L.) [24] and projected future (B(τ → µγ) < 3 × 10 −9 @ 90% C.L.) [25] indirect constraints are shaded in light and dark gray vertical bands, respectively. We observe that within the EFT approach, the CMS signal is well compatible with the non-observation of τ → µγ at the B factories and will marginally remain so even at Belle II. This is in contrast with the situation in most explicit NP models generating non-zero y τ µ or y µτ , as we demonstrate in Secs. III-V . Of course many more possibilities exist at the loop level (e.g. [26,27]). In practice however, all such models predicting sizable LFV Higgs interactions necessarily suffer from a severe fine-tuning problem, if they are to simultaneously avoid the stringent τ → µγ constraint. We demonstrate this on an explicit example in Sec. V D. A general heuristic argument however goes as follows. The The accidental cancelation required in the matching procedure for the radiative decay in order to accommodate both results is thus expected to be approximately one part in 10 3 at the amplitude level. As we demonstrate in Secs. III, IV and V, these expectations are confirmed in explicit models.
A. h → τ µ vs. h → τ e A positive experimental indication for the h → τ µ decay when combined with existing stringent experimental limits on µ − e LFV processes can be used to constrain LFV τ − e processes. In particular, the product of the h → τ e and h → τ µ branching fractions is constrained from above by the rates of µ → eγ, and µ − e conversion on nuclei. This is due to the fact that tree level Higgs decays to τ µ (τ e) depend on y µτ,τ µ (y eτ,τ e ) while the same sets of couplings contribute at the loop level to µ → eγ and µ − e conversion via diagrams with a virtual τ .
The contributions to the µ → eγ process stemming from a virtual τ are m τ enhanced with respect to diagrams with intermediate µ or e states [17] leading to where we have neglected the effects of the light lepton masses while m τ dependence is retained in branching fraction is thus sensitive to a distinct combination of the LFV Yukawas: On the other hand, µ−e conversion on nuclei is most sensitive to vector current effective operators (ēγ ν P L,R µ) (qγ ν q) which probe different combinations of LFV Yukawas. The largest contribution to the vector coupling of the proton follows from Eq. (A17) of Ref. [17] and reads, in the limit of massless light leptons, The above effective coupling is evaluated in the q 2 → 0 limit, where q is the momentum exchanged in the conversion process. As before the right-handed current coefficient g RV is obtained from the above result by replacing y ij with y * ji . Note that due to strong experimental limits on B(µ → eγ) we can at present safely neglect the dipole operator Q L,R contributions to µ − e conversion [17].
The complementary information on the LFV couplings extracted from the rates of µ → eγ and µ − e conversion is sufficient to completely correlate the Higgs LFV h → τ µ and h → τ e decays: where in the second line we have approximated Γ h Γ SM h . The best experimental limit on B(µ → e) Au < 7 × 10 −13 (at 90% C.L.) was achieved by the SINDRUM II Collaboration [29] while the best upper bound on B(µ → eγ) < 5.7 × 10 −13 (at 90% C.L.) was recently determined by the MEG Collaboration [30]. Note that with the current experimental data the sum on the right-hand to both µ → eγ and especially µ − e conversion in nuclei is expected to improve significantly.
In particular the DeeMe [31] experiment aims at a sensitivity of the order of 10 −14 while the Mu2e [32] experiment could improve existing SINDRUM II bounds by four orders of magnitude.
Given B(h → τ µ) at the percent level, such a measurement will indirectly probe B(h → τ e) at the order of magnitude 10 −5 .

III. TYPE-III TWO HIGGS DOUBLET MODEL
Several well motivated scenarios of physics beyond the SM (e.g. supersymmetric extensions or models addressing the strong CP problem via a PQ symmetry) require the addition of another SU (2) L doublet scalar to the minimal SM Higgs sector. The most extensively studied version of the THDM is the type-II model in which the absence of tree level flavor changing neutral currents (FCNCs) is imposed by the requirement that only one of the doublets couples to up-type quarks, while the other one couples to down-type quarks and charged leptons. FCNCs in the quark sector are relegated to the one-loop level via the interactions involving charged scalars, whose couplings to quarks follow the CKM structure of the SM. In this model LFV is absent as in the SM due to massless neutrinos and because only a single λ α in Eq. (6) is non-vanishing.
In order to accommodate LFV Higgs interactions one thus has to depart from this version and consider the THDM with a generic Yukawa structure (type-III THDM). In the following we restrict our discussion to the parameter space with a MSSM-like Higgs potential, for which the following relations hold [33] tan where 0 < β < π/2 is the angle of the rotation matrix that diagonalizes the mass squared matrices of charged scalars and pseudoscalars and −π/2 < α < 0 is an analogous angle for the neutral  [34,35]. Note however that generically in the MSSM, LFV radiative lepton decay contributions will be generated at the same order as discussed in Sec. II, leading to tighter constraints, c.f. [36]. We thus do not attempt to interpret the type-III THDM model parameters in terms of the underlying MSSM parameters, but consider them as essentially uncorrelated, only subject to experimental bounds and theoretical constraints such as perturbativity conditions.
The relevant part of the Yukawa Lagrangian in the charged lepton mass eigenstate basis is [33] where the Yukawa couplings can be written as The terms f i parametrize the off-diagonal charged lepton Yukawa couplings and are free parameters of the model. Finally, the coefficients x k q for H k = (H, h, A) are given by [33] x k u = (− sin α, − cos α, i cos β) , In addition to tree level exchanges of neutral Higgses we have loop contributions to FCNCs due to charged Higgs, whose interactions with leptons are parameterized by the following Yukawa couplings: A. h → τ µ Using preceding relations, at the tree-level the type-III THDM contributes to the effective Higgs couplings as and we can translate Eq. (4) into a 1 σ two-sided bound on the parameters For reference we also give the expression for the branching fraction We show in Fig. 5 the points in the plane spanned by (| τ µ | 2 + | µτ | 2 ) 1/2 and m A for which B(h → τ µ) = 0.84% when all other Higgs decay rates are held SM-like. Interestingly, the decay exhibits small dependence on β for large values of tan β. Indeed, an expansion around the decoupling limit For light pseudo-scalar masses (m A ), mixing with the heavy scalar Higgs could affect the couplings of the light Higgs to the gauge bosons. However, the overall modification factor of hW W and hZZ couplings relative to the SM values is given by sin(β − α) and starts to differ from 1 only at order O(ξ 2 ): We have checked that even for m A = 150 GeV the hW W coupling deviates less that 2 % from its

B. Constraints from τ → µγ
The coefficients c L,R in Eq. (12) are generated via one-loop diagrams with charged or neutral Higgses. In the limit m µ = 0 (for the sake of consistency also y H k µµ = µµ = 0), these coefficients are [33] c (1-loop) The suppression of the above result by two small Yukawa couplings is relaxed at the two-loop level, where Barr-Zee type diagrams require only a single LFV Yukawa coupling along with y tt .
Therefore the contributions of the neutral Higgses at the two-loop level can dominate over the one-loop result. Here we employ the Barr-Zee contributions due to neutral Higgses calculated for µ → eγ in Ref. [37] and adapted to τ → µγ in Appendix B1. The 2-loop contributions to B(τ → µγ) involve exactly the same combination of elements as B(h → τ µ) in Eq. (28).
For the charged Higgs contributions to the best of our knowledge no complete calculation exists in the literature. A partial result with the charged Higgs contribution inducing the H 0 k γγ * or H 0 k γZ * vertices has been presented in Ref. [35]. However, the diagrams with an effective H ± W ∓ γ coupling via loop of t and b are still missing. In the following we use the known partial results to estimate the order of magnitude of the constraint. The missing contributions could potentially further strengthen the τ → µγ bound but are not expected to substantially weaken it.
In order to explore the nontrivial correlation between the B(h → τ µ) and B(τ → µγ) we sample the parameter space of τ µ , µτ , τ τ for a specific choice of tan β and m A . The ranges allowed for τ µ,µτ are required to fulfill the naturalness criterium of Eq. (11) and that each of the effective couplings is within the perturbative range The 1-loop amplitude of the τ → µγ process further depends on the diagonal y h τ τ Yukawa coupling which is experimentally constrained by the searches for h → τ τ decays within ATLAS and CMS experiments. In fact both collaborations report strong evidence for the existence of this mode with the signal strengths (normalized to SM expectations) µ τ τ = 1.43 +0. 43 −0.37 and µ τ τ = 0.78 ± 0.27, observed by ATLAS [38] and CMS [39], respectively. In the following we use a naïve average of the two experimental results: In the decoupling limit (c.f. Eq. (30)) and with y h tt SM-like, we can assume SM-like Higgs production cross sections and model the departure of µ τ τ from 1 by an appropriate shift of y h τ τ away from − √ 2m τ /v (i.e. non-vanishing τ τ ).
A scenario with SM-like y h τ τ coupling corresponding to µ τ τ = 1 for fixed tan β = 10 and m A = 0.3 TeV is presented in Fig. 3 by a black narrow stripe. This scenario easily passes both experimental constraints. A perfect one-to-one correspondence between the two observables is spoiled by the asymmetry between the left-and right-handed τ → µγ amplitudes at one-loop, see Eq. (31). On the other hand, for masses m A significantly larger than 500 GeV it is not possible to reconcile both predictions with the corresponding experimental values.
Allowing τ τ to vary within the experimental bounds imposed by the h → τ τ decays, i.e.
Eq. (34), one obtains predictions for h → τ µ and τ → µγ as represented by the green band in Fig. 3. In particular, this additional freedom in τ → µγ breaks, to some extent, the strict correlation with the h → τ µ rate. The remaining correlation signals the dominance of the 2-loop contributions in B(τ → µγ) which are independent of τ τ .
The interaction vertices that induce the LFV decay of Higgs boson to τ µ pairs also contribute to the LFV decay of tau lepton to three muons, τ − → µ − µ + µ − . The current best experimental upper limit on its branching fraction has been reported by the Belle Collaboration B(τ − → µ − µ + µ − ) ≤ 2.1 × 10 −8 at the 90% C.L. [40]. This process proceeds through the tree level exchange of the neutral scalars h, H and the pseudoscalar A. Since all the relevant amplitudes are suppressed by flavor-diagonal muon Yukawa, this constraint is not competitive with τ → µγ in terms of sensitivity to τ µ and µτ .
Here we repeat the derivation of the bound on h → τ e described in Sec. II A in the context of the type-III THDM model. In this case, both µ → eγ and µ − e conversion receive contributions not from a single Higgs but from three neutral Higgses while the charged Higgs contributions are relatively suppressed by light lepton masses. The effective operator coefficients presented in Sec. II A now have to be summed over the three Higgses with proper mass and effective coupling replacements. Ultimately, the observables depend on t β ≡ tan β, m A , and the LFV parameters ij , where dependence on the first two parameters can be absorbed into an overall coefficient: Tree level LFV Higgs decays h → τ µ and h → τ e are still properly described in the effective framework of Eq. (26). In contrast with what has been found in the effective framework (c.f. Eq. (19)), the sensitivity to ij is now similar in both low energy observables. Also, the coefficients B µ→eγ,µe 0 now have a strong dependence on tan β and m A . Finally, Higgs decays are suppressed by a projection factor so that These features suppress the h → τ e process very stringently in the type-III THDM model: The above bound is monotonically decreasing with tan β and m A as can be seen on Fig. 4. Note that at the projected sensitivity of the Mu2e experiment, the right hand side of Eq. (36) would start to be dominated by the µ → eγ constraint, so an improvement expected from the planned MEG upgrade (MEGII) [41] will be essential in this context, as can be seen by comparing the two bottom-most dashed contours in the plot.

IV. VECTOR-LIKE LEPTONS
In this section we discuss the possibility of generating sizable LFV Higgs couplings by mixing the SM chiral leptons with new vector-like leptons (VLLs). Such states can appear in the low energy spectra of grand unified theories [42][43][44][45][46] and are predicted to exist in composite Higgs scenarios with partial compositeness [47][48][49][50]. In this setting the chiral leptons obtain additional couplings to the Higgs by mixing with the heavy vector-like leptons.
We first consider the addition of VLLs in a single vectorial representation of the SM gauge In this case one can prove [51] that LFV Higgs couplings are directly related to LFV Z-boson couplings where g 0.65 and c W ≡ cos θ W 0.88 are the weak isospin gauge coupling and the cosine of the weak mixing angle, respectively. In particular one obtains a relation Severe constraints on X τ µ,µτ , Y τ µ,µτ < ∼ 10 −3 from searches for τ → 3µ [52] then preclude any sizable contributions to B(h → τ µ) . 3 De-correlating LFV Higgs and Z boson couplings requires the introduction of VLLs in both SM gauge representations (Ψ E , Ψ L ) mixing with chiral leptons Integrating out the vector-like states one finds the LFV Higgs and Z boson couplings can now be completely de-correlated in the limit where the direct couplings of chiral fermions to the Higgs vanish, and the SM leptons obtain their masses exclusively through mixing with VLLs in Eq. (39).
In particular the lepton flavor off-diagonal Higgs Yukawas in Eq. (7) are then of the form [54] All the LFV phenomena are driven by this contribution and one can derive a one-to-one correlation between tree-level Higgs decay h → τ µ and the 1-loop radiative decay τ → µγ [54] B(h → τ µ) The 1-loop suppression factor of τ → µγ is not small enough to evade the experimental upper bound and at the same time accommodate the h → τ µ at the percent level, as illustrated by a pink-dotted line in Fig. 3.
This type of diagrams requires a helicity flip on one of the fermion lines. One therefore expects suppressed amplitude when all the fermions are substantially lighter than v. It turns out that size The physical Higgs can couple to the scalar ∆ or to the top quark as shown in the first two diagrams in Fig. 6. While the strength of the coupling relevant for the latter process is fixed by the top Yukawa, the former process depends on an unknown hLQLQ coupling, λv, that originates from the marginal "Higgs portal" operator, Here ∆ is the scalar LQ and H is the SM Higgs doublet.
The Yukawa couplings of ∆ 1 are given by the following Lagrangian where Q i = (u i L , d i L ) T and U i = u i R are the quark weak doublets and up-type singlets, respectively. We explicitly show flavor indices i, j = 1, 2, 3, and SU (2) indices a, b = 1, 2, with 12 = 1. Also, here y L ij and y R ij are elements of arbitrary complex 3×3 Yukawa coupling matrices. After expanding the SU (2) indices, we obtain where V CKM and V PMNS represent Cabibbo-Kobayashi-Maskawa and Pontecorvo-Maki-Nakagawa-Sakata mixing matrices, respectively. All fields in Eq. (44) are specified in the mass eigenstate basis. The Wilson coefficients of the h → τ µ effective Lagrangian (3) are then where N c = 3 is the number of colors. The relevant loop function further depends on the portal coupling λ, given in terms of Passarino-Veltman functions B 0 and C 0 . The h → τ µ decay width in the ∆ 1 LQ scenario is then given by The state ∆ 1 also contributes to the τ → µγ through leptoquark Yukawas present in h → τ µ.
The Wilson coefficients of effective Lagrangian introduced in Eq. (12) are with x t = m 2 t /m 2 ∆ 1 , in agreement with [55]. The τ → µγ branching fraction in presence of ∆ 1 depends on identical Yukawa combination as the h → τ µ decay rate: In the left-hand side panel of Fig. 7 we show the dependence of B(h → τ µ) on the leptoquark Yukawa couplings assuming the total Higgs decay width to be SM-like. Here, we turn off the Higgs portal coupling. We set the mass of ∆ 1 LQ to 650 GeV, consistent with current direct search limits at the LHC [56,57]. In the same panel of Fig. 7, we show the CMS measurement reported in [3]. The B(τ → µγ) dependence on the leptoquark Yukawa couplings for ∆ 1 case is shown in the right-hand side panel of Fig. 7 (dashed line), again for m ∆ 1 = 650 GeV.

B. Impact of the Higgs portal and its constraints
Next we study the impact of a non-zero Higgs portal coupling in the scenario with ∆ 1 . As an example, for m ∆ 1 = 650 GeV the loop function dependence on the portal coupling is CMS 1σ 10 -5 10 -4 0.001 0.010 0.100 1 10 -11 τ → µγ constraints. The ∆ 1 (∆ 2 ) case is rendered in dashed (solid) line. The following transformation needs to be applied when going from the ∆ 1 to the ∆ 2 case: y L ij → y L ji .
Thus, a positive large λ could in principle relax the leptoquark Yukawa couplings and yield sizable h → τ µ rates without violating the τ → µγ constraint. However, the Higgs portal coupling also induces corrections to the h → γγ decay and to gluon-gluon fusion (ggF) induced Higgs production with the leptoquark running in the triangular loop. The modified ggF production cross section normalized to its SM value is Here, N ∆ i is the number of ∆ i components in the weak multiplet ∆. C(r ∆ ) is the index of color representation r ∆ of ∆ and for the triplet (C(3) = 1/2). We consider heavy enough colored scalars such that the loop function is in the decoupling limit. Similarly, the modified h → γγ decay width, normalized to its SM value, is given by The sum in Eq.
Fit to the latest available LHC Higgs data including the CMS h → τ µ measurement.
If besides the λ and leptoquark Yukawa couplings, we also allow other Higgs couplings to vary (see Appendix A) we instead get the dark and light grey regions at 1 σ and 2 σ levels, respectively.
C. The ∆ 2 = (3, 2, 7/6) case The Yukawa couplings of the ∆ 2 leptoquark to SM fermions are where we explicitly show the flavor indices i, j = 1, 2, 3, and SU (2) indices a, b = 1, 2. y L and y R in Eq. (53) are arbitrary complex 3 × 3 Yukawa matrices. In the mass eigenstate basis we have where a superscript on ∆ 2 denotes the electric charge of a given SU (2) where the loop function g 1 has been introduced in the previous section. The h → τ µ decay rate is then Allowed values for the Higgs portal coupling λ can be inferred from a global fit to the Higgs data as has been done for the portal coupling of the (3, 1, −1/3) state. (See Fig. 8.) We do not attempt to repeat the same procedure for the state (3, 2, 7/6) since we do not expect a drastic change in the allowed range of λ.
The Wilson coefficients for τ → µγ are again proportional to the couplings responsible for h → τ µ: with x t = m 2 t /m 2 ∆ 2 and agree with the formulas presented in [55]. Finally, the corresponding branching fraction is given by In the left-hand (right-hand) side panel of Fig. 7 we show the B(h → τ µ) (B(τ → µγ)) dependence on the leptoquark Yukawa couplings for the ∆ 2 LQ case with the Higgs portal coupling turned off and taking m ∆ 2 = 650 GeV.
Also in this leptoquark scenario the bound on B(τ → µγ) excludes sizable B(h → τ µ) due to the strict correlation between the two observables. See the orange stripe in Fig. 3, where the portal coupling is restricted to |λ| < 1.

D. Fine-tuning solution
In this section we give an example of a phenomenologically viable model with fine-tuning.
We add to the SM a scalar leptoquark ∆ and a vector-like top quark partner T and turn on the following couplings On the other hand, T can cancel the top quark contribution to τ → µγ. We calculate the decay rate where h 2 (x) is defined in Eq. (57).
As an example we take m ∆ = 650 GeV and y L 33 y R * 32 = 0.8, which gives the best fit value for B(h → τ µ). In Fig. 9 we show how to avoid the τ → µγ constraint with an appropriate choice of x R * 32 x L 33 and m T . The scenario can be realized in both LQ models. In passing we note that if m T > m ∆ , the T → ∆ decays can produce spectacular signatures at the LHC, while at the same time avoiding traditional searches for vector-like top quark partners, especially if T does not mix with the top.

VI. CONCLUSIONS
Prompted by the recent experimental hint of h → τ µ events by the CMS Collaboration, we have carefully examined the implications of LFV Higgs decays at the percent level on possible extensions of the SM. In particular we have shown how a tentative B(h → τ µ) signal can be combined with other Higgs measurements to yield a robust lower bound on the effective LFV Higgs Yukawa couplings to taus and muons. Then we have reexamined the connection between LFV Higgs decays and LFV radiative decays of charged leptons, and demonstrated using EFT methods that the current CMS hint implies τ → µγ at rates, which could be observable at the Belle II experiment. In explicit models, the τ → µγ constraint is generically much more severe. In fact, an eventual observation of h → τ µ at at the LHC together with existing indirect constraints would already single out an extended SM scalar sector as a required ingredient in any natural explanation, the minimal example being THDM of type III. We have also examined purely fermionic SM extensions and models where h → τ µ is generated at loop level, only to show that without the introduction of extra Higgs doublets, reconciling all existing indirect constraints with percent level B(h → τ µ), when at all possible, requires a high degree of fine-tuning. Finally, we have shown how a positive signal of h → τ µ can be combined with experimental searches for µ → eγ decays and µ − e conversions in nuclei to yield robust bounds on B(h → τ e). In particular, considering only low energy Higgs EFT effects, the two LFV Higgs decay rates could still be comparable. On the other hand, the THDM III cannot accommodate both respective branching ratios at the percent level. Currently planned improvements in experimental searches for µ − e LFV processes will be able to probe the product B(h → τ µ)B(h → τ e) at the 10 −7 level in generic EFT and to order 10 −12 or better within the THDM III.
where index i represents a Higgs decay mode, while k denotes a production channel. Such observables could then easily be expressed in terms of new physics parameters. However, the actual experimental categories are never pure and contain events from different production mechanisms.
Furthermore, in order to fully reconstruct the total likelihood function, it is necessary to know all the correlations among the different categories, which are available only to the experimental collaborations.
It has been argued recently in Ref. [58] that the existing measurements for a given decay channel should be presented in terms of two-dimensional likelihoods, in which gluon-gluon fusion (ggF) and associated production with a top pair (ttH) are combined as one signal (µ (ggF +ttH) ), while vector boson fusion (VBF) and associated production with a gauge boson (VH) as another, (µ (V BF +V H) ). This decomposition is particularly useful since, on one hand, ttH is sub-dominant with respect to ggF and poorly explored in the present data set, while VBF and VH receive common corrections in wide class of NP models obeying the custodial invariance. Furthermore, the correlations are automatically provided by the experimental collaborations for a given decay channel.
Following the recommendation, ATLAS and CMS have combined different search categories for a given decay mode to provide separation into production mechanisms. The results are usually presented in 2D plots in (µ (ggF +ttH) , µ (V BF +V H) ) plane. In this case, we parametrize the likelihood with where the sum goes over the decay channels and the covariance matrices are given by We obtain the best-fit values (μ), variances (σ) and correlations (ρ) from the plots provided by the experimental collaborations. These are presented in Tab. I and Tab. II.
If the separation into production modes is not provided, we use the signal strength measurements from existing search categories. These in general target certain production mechanism, which, however, does not imply 100% purity. Inclusive categories are dominated by ggF (90%), while VBF-tagged categories can have 20% to 50% contamination from ggF. VH-and ttH-tagged h → µµ inclusive 1.6 ± 4.2 [64] categories are assumed to be pure. In this case, we write where ξ i represent contributions of the specified production mechanisms to the given category. We do not consider correlations here and add each measurement to total χ 2 separately, We approximate the likelihood for the total Higgs decay width measurement [73] by In Tab. I and Tab. II we finally summarize all the available LHC Higgs data used in our analyses.
The limits quoted in Eq. (4) are obtained after allowing only y τ µ , y µτ to be free parameter while fixing other Higgs couplings to their SM values. On the other hand, in the limits shown in Eq. (5) we further allow for arbitrary values of the usual SM tree level Higgs boson couplings; κ t , κ τ , κ b , κ V and κ c where κ V = κ W = κ Z as well as arbitrary new contributions to loop induced Higgs couplings; κ g and κ γ .