Probing Lepton Flavor Violation at the 13 TeV LHC

We investigate the bounds on tau-mu lepton flavor violation (LFV). Our main focus is on the collider constrains on tau-mu LFV. We use the Type-III Two-Higgs-Doublet-Model (2HDM) as a set up for our study. While the LFV branching fraction of the 125 GeV is well constrained by current LHC searches, the heavier neutral states could have a large branching fraction to tau and muon. We estimate the LHC reach for the 13 TeV center of mass energy with 300 $\text{fb}^{-1}$ luminosity for a neutral boson decaying into a tau and a muon. We identify parts of the LFV parameter space where the searches for heavy scalar and pseudoscalar decaying into a tau and a muon are more sensitive than the similar search for the 125 GeV boson.


Introduction
The upgraded Large Hadron Collider (LHC) opens up the possibilities to explore a higher energy scale where new physics may lie. The lepton flavor volation (LFV) is an interesting possible new physics that might show up during this next run of the LHC. In Ref. [1,2], various LFV decay channels of the 125 GeV scalar h were explored. The authors found that LHC constraints on the decay h → τ µ can be superior to the bounds from low energy experiments such as τ → µγ and τ → 3µ. This decay has been probed at the LHC run-1 [3][4][5] and early run-2 [6]. The ATLAS and CMS experiments constrained the LFV branching fraction to be BR h→τ µ < 1.43% and BR h→τ µ < 1.2% respectively. Additionally there is a 2.4σ hint of LFV branching fraction from run-1 CMS with BR h→τ µ = 0.84 +0. 39 −0.37 % [4], which is marginally compatible with the previous constraints.
The Two Higgs Doublet Model (2HDM) is one possible extension of the standard model (SM). In this model, the particle content of the SM is enlarged by an introduction of an additional scalar doublet. The extra doublet brings with it a host of interesting LHC phenomenology. Firstly, there are additional new particles that might be observed at the LHC. Assuming the model is CP conserving, these new particles are neutral heavy scalar (H), neutral pseudoscalar (A) and a pair of charged scalars. Secondly, the new scalar doublet introduces new Yukawa couplings to the fermions which could give rise to LFV couplings. The LFV couplings, in turn, lead to LFV decays of the neutral scalars and pseudoscalar 1 , in particular the decay into tau and muon. Hence this model can be a simple UV completion of LFV decay considered in Ref. [1,2]. Many other models can incorporate the LFV decay, for example, the minimal supersymmetric Standard Model and its extension [7,8], models with an electroweak triplet [9,10], models with vector-like heavy lepton [11,12], composite Higgs model [13], and the Little Higgs models [14,15]. As has already been mentioned, the LFV h → τ µ has been searched for at the LHC and the branching ratio is tightly constrained. However, the corresponding LFV branching ratio of the heavy neutral particles can, in principle, be large. Thus it is possible that the larger branching ratio compensates the smaller cross section of heavier scalars, making it possible to probe these decays at the LHC. These LFV decay searches could explore more parameter space compared to the current h → τ µ search alone [16][17][18]. While no specific LHC estimates have been given in refs. [16,17], ref. [18] recasts the CMS h → τ µ search to incorporate the heavy CP-even scalar H. By considering an optimistic scenario, in which the heavy scalars couple to the top quark and is produced copiously through gluon-fusion, the authors of ref. [18] found that at the 8 TeV LHC the inclusion of H → τ µ search excludes more parameter space of 2HDM.
In this paper we examine the bounds on LFV decays at the 13 TeV LHC with a luminosity of 300 fb −1 . Moreover, we fully explore the parameter space of the 2HDM by including the pseudoscalar A that has not been taken into account in the previous works. In addition, we consider various possibilities of heavy resonances couplings to SM particles in the context of type-III 2HDM. Hence we cover both the optimistic scenarios, where the heavy resonances production cross-section and their LFV branching fractions are large, and the pessimistic scenario (small production cross-sections and small LFV branching ratios) for discovering the LFV via the heavy resonance searches.
The paper is structured as follow. In Sec. 2 we review the Higgs phenomenology in 2HDM relevant for our LFV analysis. We then discuss the current constraints on LFV in the Higgs sector, including both the direct and indirect constraints, in Sec. 3. In Sec. 4 we study the collider constraints on the Higgs LFV decays. We also give the projected bound for the 13 TeV LHC with 300 fb −1 luminosity. We then conclude in Sec. 5.

Type-III 2HDM
2HDM is one of the most studied extensions of the SM, for a review see Ref. [19]. There are many realizations of the model in the literature. They can be characterized by the structure of Yukawa couplings. In Type-I, Type-II and Type-X 2HDM, each fermion type (up-type quarks, down-type quarks and leptons) is coupled to only one scalar doublet. Thus there is no flavor violating Yukawa couplings of the neutral scalar boson in this case [20,21] 2 . In type-III 2HMD, however, both of the scalar doublets couple to all the fermions. As a result, there is a possibility of flavor violation in the neutral scalars Yukawa couplings. Therefore, we will focus on the type-III 2HDM in this work.

Conventions and notations
In this section we set up our conventions and notations. The two SU (2) L scalar doublets are taken to have hypercharge Y = 1/2. With this convention, the electric charge generator is Q = σ 3 2 + Y where σ i are the Pauli matrices. The scalar potential consistent with SU (2) In general, M 2 12 , λ 5 , λ 6 and λ 7 are complex while the rest of the parameters are real. However, since in this work we are interested in the simplified scenario where CP is a good symmetry, we will take M 2 12 , λ 5 , λ 6 and λ 7 to be real. We leave the CP violating scenario for a future work.
In this paper we use the Higgs basis [23] where the vacuum expectation value (VEV) resides only in Φ 1 . In the Higgs basis, the fields Φ 1 and Φ 2 can be expanded as where v is the VEV, G ± and G 0 are the would be Goldstone bosons, H ± is the charged Higgs, A is the neutral CP-odd Higgs, φ 1 and φ 2 are the neutral CP-even Higgs. Minimizing the potential leads to the relations The fields H ± and A are mass eigenstates with masses The fields φ 1 and φ 2 are in general not mass eigenstates. They are related to the mass eigenstates h and H by where mixing angle α is given by (2.6) The scalar masses are We take the h to be the 125 GeV scalar resonance discovered at the LHC. The Yukawa sector in the Type-III 2HDM is given by where m f are fermion masses, Y f are the Yukawa coupling matrices, V is the CMK matrix and the indices i, j, k run over fermion families. The scalar doubletΦ is defined asΦ = iσ 2 Φ * . The fermion doublets are taken to be Note that the fermion fields L(R) , u L(R) and d L(R) in Eq. (2.8) and (2.9) are in the mass eigenbasis. After electroweak symmetry breaking, the Yukawa couplings in the physical basis read where P R,L = (1 ± γ 5 )/2 are the chiral projections, f runs over fermion types ( , U and D) and (2.11) Notice that the flavor violating Higgs Yukawa couplings are encoded in the Yukawa matrices Y , Y U and Y D .
The Yukawa matrices Y f 's affect the production cross-sections and decay rates of the Higgs bosons. Since we're interested in the lepton flavor violating (LFV) decay of the neutral scalar φ → τ ± µ ∓ , we will restrict our attention to the neutral scalars sector.

Production cross-sections of neutral Higgs bosons
The production cross-sections of the neutral Higgs bosons, φ, in each channel can be most conveniently described in terms of the would be SM Higgs boson cross-sections. They are [19] where q = U, D, δ h = cos α, δ H = sin α, δ A = 0, i runs over generation index, τ x = 4m 2 x /m 2 φ and the loop functions A H 1/2 (τ ) and A A 1/2 (τ ) are given in Eq. (A.1). In our analysis below, we use the would be SM Higgs boson cross-sections provided by the LHC Higgs Cross Section Working Group [24].
The pseudoscalar A can only be produced via the Yukawa interactions. Thus its production depends most sensitively on Yukawa coupling matrix Y U and Y D . Fig. 1  For the neutral scalar H, its production cross-section also depends on the mixing angle α. Fig. 2 Y U cc =0.1 Figure 1: The production cross-section of the pseudoscalar A as a function of its mass. The gluon-fusion channel is shown in red and the ttA channel is shown in orange. On the left pane Y tt U = 0.1 while on the right pane Y cc U = 0.1. Other Yukawa couplings are taken to be 0.

Decays of neutral Higgs bosons
The decays of the neutral Higgs bosons, φ, can be expressed in terms of the would be SM Higgs decays. They are where δ h = cos α, δ H = sin α, δ A = 0, N c is the number of color, Q f is the electric charge of fermion f and the loop functions A H 1 (τ ) is defined in Eq. (A.1). For the φ → γγ decay, there is also a contribution from the charged Higgs loop. However, this contribution is small thus we will drop it from our analysis. In our analysis below, we use the would be SM Higgs boson branching ratios provided by the LHC Higgs Cross Section Working Group [24].
The off-diagonal elements of the Yukawa coupling matrices in Eq. (2.10)   can be written as where in the last line we make the approximation m φ m f . Since our main interest in this work is on the LFV decays involving tau and muon, therefore the only nonzero off-diagonal elements of Y 's that we consider are Y τ µ and Y µτ .
In addition to the decay channels listed in Eq. (2.14), the pseudoscalar A can decay to a  Figure 3: The branching ratios of the pseudoscalar A as a function of its mass. In both plots, Y τ µ = Y µτ = 0.1 and sin α = 0.1. The A → τ µ channel is shown in blue, tt channel in orange, hZ channel in magenta and gg (cc) channel in red. On the left pane, Y tt U = 0.1 while on the right pane Y cc U = 0.1. Other Yukawa couplings are taken to be 0. scalar and a Z boson. The partial decay width for A → hZ is given by The Γ A HZ can be obtained by making a replacement sin α → cos α and m h → m H . If the mass m H m A , as is the case when an approximate SO(3) symmetry is imposed on the scalar sector [25], the decay channel A → HZ is closed. Fig. 3 shows the branching ratio of A as a function of its mass for the cases Y tt U = 0. is kinematically open. We can parametrize this decay channel by introducing the coupling λ φφH such that L ⊃ 1 2 λ φφH vHφ 2 for neutral resonance φ. For the charged Higgs, the coupling is defined without a 1/2 factor. Then, the partial decay width H → φφ For the case of H → H + H − , there is an extra factor of 2 dues to H + and H − being distinct particles. We note the coupling λ φφH depends strongly on the scalar potential. Fig. 4 Figure 4: The branching ratios of the heavy neutral scalar H as a function of its mass. The decay into τ ± µ ∓ is shown in blue, the ttA in orange, the digluon in red, the W ± W ∓ in green, the ZZ in magenta, the hh in gray, the bb in brown, the cc in cyan and the τ + τ − in purple. In all plots, Y τ µ = Y µτ = 0.01, sin α = 0.1 and λ hhH = 1. On the upper left, Y tt U = 0.1, on the upper right Y tt U = 0, on the bottom left Y cc U = 0.1 and on the bottom right Y τ τ = 0.01. Other Yukawa couplings are taken to be 0.

Lepton flavor violation in the Higgs sector
The off-diagonal elements in the Yukawa matrice Y f 's induce flavor violating decays of the Higgs bosons. In this paper, for simplicity, we will assume that flavor violations reside only in the lepton sector. Moreover, we will focus our attention on the τ − µ LFV. Thus the only non-zero off-diagonal entries of Y that we consider are Y τ µ and Y µτ .
The couplings Y τ µ and Y µτ can be probed directly at the LHC or indirectly via low energy precision measurements. Here we summarize the current constraints on Y τ µ and Y µτ .

Direct constraints on Y τ µ and Y µτ
The coupling Y τ µ and Y µτ lead to flavor violating decays φ → τ µ for the neutral Higgs bosons. In particular, for the 125 GeV resonance, h, the partial width for such a decay is This decay mode has been searched for at the LHC by the CMS collaboration. Currently, CMS has placed an upper bound on the branching ratio at Br(h → τ µ) < 1.20% [6]. This bound translates to the bound on the mixing angle α and the coupling Y τ µ and Y µτ as shown by the solid orange line in Fig. 5. Notice that the branching ratio h → τ µ is tightly constrained by the CMS limit. However, the branching ratio H → τ µ (or A → τ µ) can still be large. Figure 6: The one-loop diagrams contributing to τ → µγ decays induced by the Higgs bosons with flavor violating Yukawa couplings Y τ µ and Y µτ .

Indirect constraints on Y τ µ and Y µτ
Extensive analyses of flavor constraints on the Yukawa structure of the Type-III 2HDM have been carried out in Ref. [26]. Here we'll focus on the flavor observables relevant for tau-muon LFV.

τ → µγ
The flavor violating couplings Y τ µ and Y µτ lead to a rare decay τ → µγ. Their contributions can be computed by first matching to the effective operators [1] where F αβ is the U (1) EM field strength tensor. In terms of the Wilson coefficients c L and c R , the decay rate for τ → µγ can be written as The experimental bound on the branching ratio is Br(τ → µγ) < 4.4 × 10 −8 [27]. The Wilson coefficients c L,R get contributions from both one-loop and two-loop diagrams. At the one-loop level, the contributions arise from diagrams shown in Fig. 6. Their contributions are [28] where we have assumed |y τ τ ,φ | |y µµ ,φ |. Notice that c L and c R are related by Y ij ↔ Y ji * . The coefficients c L and c R also get corrections from two-loop processes which can be as large as the one-loop contributions. The expressions for the two-loop contributions, ∆c L and ∆c R , are given in App. B. In the case that Y τ µ = Y µτ are the only non-vanishing Yukawa couplings, the bound from τ → µγ is shown in the solid black line in the left pane of Fig. 5.
The Y τ µ and Y µτ contribution to the muon magnetic dipole moment.

Muon magnetic dipole moment
The Y τ µ and Y µτ coupling also contribute to a magnetic dipole moment for the muon, a µ ≡ (g µ − 2)/2, as shown in Fig. 7. Their contribution can be easily translated from the result in Ref. [1] δa where φ = h, H and A. Note in the above expression we have dropped terms suppressed by m µ /m τ and m τ /m φ . The discrepancy between SM prediction and the measured value is [27,29] This can be used to bound the LFV contribution to a µ , δa µ ≤ 4.47 × 10 −9 at 95% C.L..

Muon electric dipole moment
In effective theory, the muon electric dipole moment is described by where φ = h, H and A and we have dropped the term suppressed by m µ /m τ and m τ /m φ . Since in this work we are interested in the simple case where Y τ µ and Y µτ are real, the muon electric dipole moment constraint does not apply to our scenario. Figure 8: The diagrams contributing to τ → µ + − decays. The diagram on the right arises from the effective dipole operators in Eq. (3.2).
The flavor violating couplings Y τ µ and Y µτ also lead to a decay τ → 3µ as well as τ → µe + e − as shown in Fig. 8. Their contributions can be easily computed by matching onto effective operators. The relevant effective operators for τ → 3µ are the dipole operators in Eq. (3.2) and the 4-fermion operators where {x, y} = {L, R} and we assume LFV resides only in the tau-mu couplings. For the case of a real Yukawa matrix Y and dropping terms suppressed by m µ /v, we get (3.11) The expressions for c e xy can be obtained by a replacement Y µµ → Y ee . In terms of these Wilson coefficients, the doubly differential partial width for where m 2 ij = (p i + p j ) 2 . Note in the above expression we set m µ = 0. Similarly, the doubly differential partial decay width for Numerically, they are As can be seen from the Wilson coefficients in Eq. (3.11), Γ τ →3µ and Γ µe + e − depends on the coupling Y µµ and Y ee . These couplings are constrained by the direct h → µµ [31] and h → ee [32] search with sin α Y µµ 2.13 × 10 −3 and sin α Y ee 1.76 × 10 −3 . Current experimental bounds on the τ → 3µ and τ → µe + e − decays are [27] Br(τ → 3µ) = 2.1 × 10 −8 , (3. 16) In the limit that Y ee = Y µµ = 0, these bounds translate to the limit on the LFV couplings Y τ µ and Y µτ as shown in Fig. 5.

τ → µM
In this subsection we discuss low energy constraints from the decay of τ → µM , where M represents a light meson. Tau decays into a muon and a pseudoscalar meson (π, η, η ) constrains the coupling of the pseudoscalar, A, to light quarks, in addition to the LFV coupling (Y τ µ and Y µτ ). We follow Ref. [30] to calculate the branching fraction of LFV decays τ → µM . The expression for the decay width are given in Appendix C. The decay of τ → µπ constrains the coupling the pseudoscalar to the up and down quarks. We can compare these bounds against the bounds obtained from the kinematics of the light Higgs (h) productions. These are y uu U,h < 0.011 and y dd D,h < 0.013 [33]. The coupling of the pseudoscalar to the strange quarks are best constrained by the decay τ → µη. 3 This bound can also be compared against bounds from the LHC light Higgs precision measurement: y ss D,h < 0.029 [34]. Figure 9 shows the bounds on Y qq Q from the LFV τ → µπ and τ → µη decays for Y τ µ = 0.01. In the figure, we also show the bounds from the light Higgs precision measurements for sin α = 0.1.
The light and heavy scalar coupling to light quarks can be probed using the decay τ → µπ + π − and τ → µρ. We follow Ref. [30] to calculate the partial decay width of the tau.  Figure 9: The solid blues line in these plots show the bounds for Y uu U , Y dd D and Y ss D from τ → µP as a function of m A for Y τ µ = 0.01. Here P represents the pseudoscalars mesons: π or η. The region above the lines are excluded by the searches. The shaded regions are excluded from the kinematics of light Higgs productions for sin α = 0.1. In calculating these bounds, we take only one of the Y 's to be nonzero, while the others are zero.
We have found that the bounds from τ → µπ + π − are slightly stronger than the one from τ → µρ. Hence we only show the τ → µπ + π − results in this paper. In addition to the scalar couplings to the light quarks, the bounds also depends on sin α, Y τ µ and m H . Figure  10 shows the bounds for Y τ µ = 0.01 and sin α = 0.1. The constraints from the light Higgs precision measurements are also shown in the plots. (c) Bounds on Y ss D from τ → µπ + π − Figure 10: The solid blue lines in these plots show the bounds for Y uu U , Y dd D and Y dd D from τ → µπ + π − as a function of m H for Y τ µ = 0.01 and sin α = 0.1. The region above the lines are excluded by the searches. The shaded regions are excluded from the light Higgs precision measurements. In calculating these bounds, we take only one of the Y 's to be nonzero, while the others are zero.

Analysis
In this section we discuss the LHC search for the LFV decay of the scalar and pseudoscalar resonances. We follow closely the analysis of CMS run-1 [4]. In the CMS analysis, two decay channels of tau are discussed: tau decays into electron and neutrinos (denoted by µτ e channel) and the hadronically decay tau (denoted by µτ h channel). The background for the µτ h channel is dominated by j(W → µν), which the jet is misidentified as τ h . Since we cannot simulate the jet misidentifications accurately, we will not consider the µτ h channel in our work. This exclusion weakens our expected bound, hence our estimated bound is a conservative one. The µτ e channel requires exactly one muon with p T > 25 GeV and |η| < 2.1, and one opposite charge electron with p T > 10 GeV and |η| < 2.3. The events are categorized according to the number of jets in the event. The jets are required to have p T > 30 GeV and |η| < 4.7. The cuts for each jet categories are defined in Table 1. In order to distinguish the signal from the background, the collinear approximation is used. The collinear mass of the system is defined as where m eµ is the invariant mass of the electron and the muon and (4. 2) The main backgrounds for this search are Z → τ τ , W W , and ZZ. 4 The signal and the background event samples are generated using Madgraph 5 [35] followed by parton shower, hadronization and matching simulations in PYTHIA 6 [36]. Delphes 3 [37] was used to simulate the detector environment. We include pileup in our simulation with 21 pileups per bunch for 8 TeV simulation and 40 pileups per bunch for 13 TeV run [38]. The signal models are generated using Feynrules [39]. The comparisons between our simulation and the 8 TeV CMS simulation [4] are shown in Fig. 11, while the m col distribution for the signal is given in Fig. 12.
For the estimated 13 TeV reach, we define the signal regions to be m φ − ∆ < m col < m φ + ∆, where m φ is the mass of the resonance we are interested in. For each m φ , the value    of ∆ is varied to get the best estimated bound. The 95% C.L. bound is found by solving for N sig that satisfies 5 where N back is the number of estimated background events and the probability P follows the Poisson distribution, P (k, λ) = λ k e −λ k! . N sig is related to the signal cross section, σ, and the branching fraction BR φ→τ µ , by N sig = σ BR φ→τ µ L , where L is the luminosity and is the acceptance and efficiency of the detector estimated from the simulation. In our analysis  below, we focus on the gluon-fusion production channel because it is the dominant production mechanism for both the heavy scalar and the pseudoscalar as can be seen in Fig. 1 and 2. The estimated bounds for the 13 TeV LHC with 300 fb −1 luminosity are given in Fig. 13. In deriving these bounds, we only consider the 0-jet category as it contains signal predominantly produced by gluon-fusion [4].

Results
The constraints on LFV from the φ → τ µ search (φ = H, A) can roughly be categorized into two cases according to its production cross-section. If the production cross-section is large, the bounds are expected to be strong. On the contrary, if the cross-section is small, the bounds are expected to be weak. For comparison, we also give the estimated bounds from h → τ µ search.

Large cross-section case
We begin with an optimistic scenario where the production cross sections of the heavy scalar and the pseudoscalar are large with a significant LFV branching fraction. This scenario can be achieved by taking Y tt U large while keeping the other Y 's vanishing (except Y µτ and Y τ µ ). However, Y tt U is constrained by low energy experiment and hence cannot be arbitrary large. For example, the coupling Y tt U , in a combination with Y τ µ and Y µτ , can be constrained by the τ → µγ measurement, see for example, Eq. (B.1) and (B.3). Moreover, Ref. [40] analyzed constraints from b → sγ and ∆M B d and found that |Y tt U | < 1 for m H + 500 GeV. Additionally, the LHC Higgs data provide another handle on Y tt U . A combined analysis of CMS and ATLAS from the LHC run-1 data found the gluon-fusion signal strength µ ggF ≡ σ ggF (σ ggF ) SM = 1.03 +0. 16 −0.14 [41] 6 . It is predicted that this signal strength can be measured with a precision of 6% at the 14 TeV LHC with 300 fb −1 luminosity [42]. The bound from the gluon-fusion signal strength can be seen in Fig. 14. In the figure we also include the bounds on sin α from h → W W * measurements from the LHC run-1, σ V BF BR h→W W (σ V BF BR h→W W ) SM = 0.84 +0.4 −0.4 [41], together with the future estimates for the 14 TeV LHC with 300 fb −1 luminosity [42].
The cross section of the light scalar, h, is determined by y tt U,h which depends on the mixing angle α and Y tt U as shown in Eq. (2.11). Current data show that the branching fraction BR h→τ µ is smaller than a few percents, hence it does not significantly affect the total decay width. Therefore the branching fraction BR h→τ µ mostly depends on y τ µ ,h = Y τ µ sin α and y µτ ,h , see Eq. (3.1). Hence the constraints from the h → τ µ decay can be displayed in the y tt U,h -y τ µ ,h plane. Fig. 15a shows our projection of such a bound for the 13 TeV LHC with 300 fb −1 luminosity. In Fig. 15b, we take sin α = 0.05, 0.1, 0.5 as our benchmark values. For small values of sin α, the bounds on Y τ µ does not significantly depend on Y tt U . The production cross section of the heavy scalar depends on m H and α and Y tt U through y tt U,H , while its LFV branching fraction also depends on Y τ µ , Y µτ and λ hhH (if the decay H → hh is open). Given the number of free parameters, we pick a set of benchmark points to show the LFV H bounds. The estimated constraints in the Y τ µ -m H plane are shown in Fig. 16   channel of the H opens, the H LFV bound becomes weaker due to a smaller branching fraction for H → τ µ. This is especially pronounced in the case of a small Y tt U and a small mixing angle α where the production cross-section of the H is small. In the figure we also show the LFV h bounds for comparison. Note also that for each value of sin α in Fig. 16, the h LFV bound is approximately the same for the values of Y tt U under consideration. For the cases of small mixing angle, sin α = ±0.05 and ±0.1, the BR h→τ µ is small compared to BR H→τ µ . Hence the H LFV bounds are generically stronger for a wide range of m H , see Figs. 16a, 16b, 16d and 16e. Note that in Figs. 16d and 16e, there is no H LFV bound for Y tt U = 0.05 because the H production cross-section is too small. Also from the figure, it's clear that for the case Y tt U = 1 the H LFV bound is particularly strong because of a large production cross-section of the H. On the other hand, for the case sin α = ±0.5, only small value of Y tt U is viable, see Fig. 14. Thus, we only consider the case Y tt U = 0.05 for sin α = 0.5, Fig. 16c and Y tt U = 0.05 and 0.1 for sin α = −0.5, Fig. 16f. In these cases the H production cross-section is small while the h LFV bound is strong because of a large BR h→τ µ , see Fig. 15b. As a result, we estimate the LFV H searches only provide a better constraint for m H 160 GeV in the case of sin α = ±0.5. It should be noted that as the mixing angle decreases, the LFV h bound gets weaker. Hence the search for LFV H decays become more and more relevant.
Having discussed the heavy scalar LFV bounds, now we turn our attention to the pseudoscalar LFV bounds. The production cross section of A depends only on Y tt U and m A while the LFV branching ratio also depends sin α, Y τ µ and Y µτ . Similar to the case of the heavy scalar, the A LFV branching fraction decreases as new decay channels (A → hZ and A → HZ) open. The estimated LFV pseudoscalar bounds for various benchmark points are shown in Fig. 17. Note the partial decay width of A depends on sin 2 α, thus we only consider the positive values of sin α. Let's first consider the case where the decay channel A → HZ is closed (the top rows of Fig. 17). For the case of a small mixing angle, sin α = 0.05 and 0.1, the BR h→τ µ is small, thus the bounds from the A LFV search are stronger than the h LFV bounds for a wide range of m A . On the other hand, for a large mixing angle, ie. sin α = 0.5, the viable value of Y tt U is small. In this case the BR h→τ µ is large while the production crosssection of A is small. Thus the bound from the h LFV search is much stronger than the A LFV bound. In the case that the decay channel A → HZ is open, the LFV branching fraction reduces significantly. Hence for low values of Y tt U , the LFV A bounds is stronger than LFV h bounds only for m A m H + m Z . For a bigger value of Y tt U , the A production cross-section increases. As a result, the LFV A bounds can be stronger for higher values of m A . Finally, we note that in the SO(3) limit where m H = m A , the pseudoscalar bounds is more constraining than the heavy scalar bounds.

Small cross section case: mixing
For the rest of this section we'll consider the opposite limit where Y tt U is vanishing. In this limit, the production cross-section of the heavy resonances, H and A, are small. As a result, one would expect the bounds on LFV from H and A searches to get weaker. Nevertheless, the LFV bounds from the heavy resonances searches can still be more constraining than the bound from the 125 GeV search in some parts of parameter space.
We starts with the simplest case where all the Y 's are vanishing except for Y µτ and Y τ µ . In this scenario the pseudoscalar cannot be produced, hence it offers no LHC bounds. The heavy scalar, H, on the other hand, can still be produced through the neutral scalars mixing. Fig. 18 shows the 13 TeV with 300 fb −1 LHC bounds on the coupling Y µτ and Y τ µ as a functions of m H with sin α = 0.05, 0.1 and 0.5 respectively. For comparison, similar bounds from the 125 GeV search are also given. From the plot, we can see that in the case of sin α 0.1, the H LFV search for m H 160 GeV yields stronger bounds than the h LFV In our analysis so far we have taken all the Yukawa couplings to vanish except Y µτ and Y τ µ . One could argue this is an optimistic scenario because the LFV branching faction of the heavy resonances are large, see for example Figs. 3 and 4. Introducing other non-zero Yukawa couplings, in principle, would dilute the LFV branching ratio and hence weakens our LFV bounds. Below, we relax this assumption and explore the consequences of allowing other Yukawa couplings to be non-vanishing.

Small cross section case: mixing and Y cc U
We begin by introducing non-zero quark Yukawa couplings Y U and Y D . The case of nonzero Y tt U has already been discussed in subsection 4.2.1, hence we will not consider it here. Similarly, we will not consider the off-diagonal elements of Y U and Y D since they are severely constrained by FCNC experiments. The couplings Y uu U , Y dd D and Y ss D are severely constrained, as discussed in Section 3.2.5. The coupling Y bb D is also tightly constrained from the h → bb search at the LHC. Therefore, this leaves the Y cc U as the least constrained coupling in this scenario. Ref. [43] estimates the bounds on the light Higgs coupling to the charm quarks from a global fit of the LHC Higgs data and found that y cc U,h 0.045. A non-zero value of y cc U,h also contributes to the h production cross-section through a charm loop and is also constrained by the Higgs measurements. The bounds on Y cc U and sin α from the above constraints are shown in Fig. 19.  Figure 19: The excluded parameter space in the plane of Y cc U and sin α. The dark blue and dark red regions are the excluded regions from the LHC run-1 ggh coupling and W W h coupling measurements respectively. The light red regions are the predicted exclusion regions from the W W h coupling measurements at the 14 TeV LHC with 300 fb −1 luminosity. The yellow region is the cch coupling bounds taken from [43].
As in the case of nonzero Y tt U , the LFV branching fraction of the light Higgs depends on Y τ µ , Y µτ and sin α through y τ µ ,h and y µτ ,h . In this scenario, the charm-loop contribution to the h production cross-section is enhanced by a non-zero Y cc U . However, in the case of h, the top-loop contribution dominates over the charm-loop contribution. The LFV bounds from the light Higgs search are shown in Fig. 20. These bounds do not vary greatly with Y cc U , especially for small values of sin α, since in these cases the charm-loop contribution is not significantly enhanced.
For the heavy scalar, H, its main production channel is through the gluon-fusion via the top-loop and the charm-loop. For a small mixing angle, sin α = ±0.05 and ±0.1, the toploop contribution and the charm-loop contribution are comparable, while for a large mixing angle, sin α = ±0.5, the top-loop contribution dominates over the charm-loop. Hence the H cross-section does not vary significantly over the range of Y cc U compatible with the LHC Higgs measurements (Fig. 19)  For the pseudoscalar, its main production channel is through the gluon-fusion via a charm-loop. Thus the A cross-section grows with |Y cc U | 2 . Its LFV branching ratio, on the other hand, decreases with |Y cc U | 2 because of the increasing A → cc partial width. However, the gain in production cross-section outweighs the drop in BR A→τ µ . Hence, the A LFV bounds get more constraining for larger values of Y cc U . Fig. 22 shows the estimated bounds for LFV pseudoscalar search for sin α = 0.05, 0.1 and 0.5. In obtaining these bounds, we assume the decay channel A → HZ is closed. If A → HZ is open, the bound would be even less constraining. For comparison, we also show the corresponding LFV h bounds. For each value of sin α, the Y cc U values are chosen so that they are compatible with the LHC Higgs data (Fig. 19). From the plot, we can see that, unless the value of sin α is small, the LFV h search is more constraining than the LFV pseudoscalar search. Even for sin α = 0.05, the LFV A search is only more constraining for m A 170 GeV. Given this estimated results, we do not consider the case which the A → HZ decay channel is open. Finally we note that, contrary to the nonzero Y tt U case, in this case the heavy scalar bounds is stronger than the pseudoscalar bounds in the SO(3) limit where m A m H .

Small cross section case: mixing and Y τ τ
In the previous subsection, we reduce the LFV branching fractions by introducing Y cc U . In that scenario, the production cross-section in the gluon-fusion channel for the scalars are slightly boosted by charm loops. If instead one introduces a non-zero Y τ τ , while keeping other Y 's zero (except Y µτ and Y τ µ ), the production cross section of the scalars does not get affected 7 . However, in this scenario the pseudoscalar does not couple to quarks hence it will not get produced at the LHC. The value of Y τ τ is constrained by the measurement of h → τ τ decay which has the signal strength µ τ τ = 1.11 +0. 24 −0.22 [41]. Since we want to consider the most pessimistic scenario here, we take the maximum allowed value of Y τ τ . The estimated bounds from the h and the H LFV searches are given in Fig. 23. In this case a large value of mixing is needed to produce enough H. For a large, but still allowed, value of sin α, the H LFV search provides better bounds than the h search only in the low mass region, m H 160 GeV. 7 Non-zero value of Y µµ and Y ee will have the same effect. However these couplings are constrained by τ → µ + − and direct h → + − searches as discussed in Section 3.2.4.  Figure 23: The estimated bounds from LFV decays of h and H in the case of nonzero Y τ τ and Y µτ for 13 TeV 300 fb −1 LHC. The value of Y τ τ is taken to be the maximum value that still satisfies a 95% C.L. of µ τ τ = 1.11 +0.24 −0.22 [41]. The region above each curve is excluded.

Conclusions and discussions
The recently upgraded LHC offers an opportunity to probe physics at a higher energy scale and a promise of new physics discoveries. One of the generic phenomena in new physics scenarios is LFV. So far, the experimental efforts to probe LFV at the LHC has been focused on the 125 GeV scalar, h, decaying into tau and muon. However, since LFV requires new physics and new physics usually comes with extra particles, it could be beneficial to search for LFV in the decays of these new particles as well. Thus in this paper we explore the possibility of utilizing these additional neutral particles in LFV search.
In our work, we focus on the CP conserving Type-III 2HDM as a concrete setup for new physics with LFV. In this scenario, the LFV decay into tau and muon originates from the Yukawa couplings Y τ µ and Y µτ , see Eq. (2.8). These couplings are constrained by both the low energy measurements and the collider experiments. However, we found that the constraints from low energy measurement is typically weaker than the collider constraints, see for example Fig. 5.
The Yukawa couplings Y τ µ and Y µτ correlate the h → τ µ decay with those of the heavy scalar (H → τ µ) and the pseudoscalar (A → τ µ). We simulate the LHC reach for 13 TeV center of mass energy with 300 fb −1 luminosity and find that the H and A search offer complimentary bounds on the LFV parameter space to the traditional h search in many of our benchmark scenarios, see Sec. 4.2. Moreover, we have considered various scenarios of the Yukawa couplings and shown that even in the pessimistic case of a small production cross-section and a small LFV branching fraction, the heavy resonance searches are still well motivated. Thus combining all three searches give the best possible bound on the LFV parameter space.
We assume in this work that CP is conserved, in particular, we take all the Yukawa couplings to be real. One could relax this assumption and consider complex Yukawas. In this case, one can study CP violation in the tau and muon decay. In fact, CP violation in h → τ µ has been studied in Ref. [44]. However, one can expand such a study to include the H and the A. We leave this for a possible future project.

B Wilson coefficients for τ → µγ
The two-loop contributions to the Wilson coefficients c L and c R had been worked out in Ref. [46]. Below we translated their results into our notation. The two-loop contributions can be split according to the particles running in the loop.
where z xy ≡ m 2 x /m 2 y , (δ h , δ H , δ A ) = (cos α, sin α, 1), and the loop functions are The 2-loop contributions to c R can be obtained from ∆c L with the replacement y ji * f,φ → y ij f,φ .

C Decay width of τ → µM
We follow Ref. [30] to calculate the decay width of τ → µM . The decay width of τ → µπ in the limit of m µ = 0 is given by where we take m τ = 1.777 GeV, m π = 134.98 MeV and f π = 130 MeV. The decay width of τ → µη in the limit of m µ = 0 is given by where we take m η = 547.86 MeV, h q η = 0.001 GeV 3 and h s η = −0.055 GeV 3 . The decay width of τ → µη in the limit of m µ = 0 is given by where we take m η = 957.78 MeV, h q η = 0.001 GeV 3 and h s η = −0.068 GeV 3 . In the equations above, we have The differential decay width of τ → µπ + π − as a function of the pion pair invariant mass squared, s = (p π + + p π − ) 2 , is given by dΓ τ →µπ + π − ds = s − 4m 2 π (m τ − s) 2 512πm τ √ s C u SL + C d SL Γ π (s) + C s SL ∆ π (s) 2 + (L → R) , (C.5) where the hadronic form factors Γ π (s) and ∆ π (s) are taken from Ref. [47]. The coefficients C q SL and C q RL are given by The decay width of τ → µπ + π − can be calculated by integrating Eq. C.5 for 4m 2 π ≤ s ≤ (m τ − m µ ) 2 . The width of τ → µρ is calculated by integrating for Eq. C.5 for 587 MeV ≤ √ s ≤ 962 MeV [48]. In the equation above, we have ignored the contributions from the fermion masses (the first terms in Eq. 2.11).