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 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 refs. [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 -1 -
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 refs. [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 section 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 section 3. In section 4 1 The charged scalar also has a LFV decay into ν , where ν is the neutrino with a different flavor than . However since none of the LHC experiments can detect neutrino flavors, this LFV decay will be indistinguishable from the flavor conserving decay ν . 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 section 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 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 Ref. [22] shows that a small perturbation in the lepton Yukawa structure of the Type-X 2HDM could lead to an observable h → τ µ decay.

JHEP05(2017)055
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 φ 1 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  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.
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 . Figure 1 shows the production cross-section of A, σ A , at the 13 TeV LHC as a function of its mass, m A , for the cases Y tt U = 0.  For the neutral scalar H, its production cross-section also depends on the mixing angle α. Figure 2 shows the production cross-section of H, σ H , at the 13 TeV LHC as a function of its mass, m H , for the cases Y tt U = 0.1 (top row, relevant for the analysis in section 4.2.1) and Y cc U = 0.05 (bottom row, relevant for the analysis in section 4.2.3) with the mixing angle sin α = 0.1 and 0.5.

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) lead to flavor violating decays of the neutral Higgs bosons. The partial decay width into final states f if j 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 scalar and a Z boson. The partial decay width for A → hZ is given by  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, 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. Figure 3 shows the branching ratio of A as a function of its mass for the cases Y tt U = 0.1 (relevant for the analysis in section 4.2.1) and Y cc U = 0.1 (relevant for the analysis in section 4.2.3). The heavy neutral scalar H can also decay into a pair of lighter Higgs bosons, φ = h, A and H + H − if it 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 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. Figure 4 show the branching ratios of H as a function of its mass for 4 different benchmark scenarios. Each benchmark is relevant for our analysis in section 4.2. In all these plots, we assume the decay H → AA and H → H + H − are kinematically closed.

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 figure 5. Notice that the branching ratio h → τ µ is tightly Figure 6. The one-loop diagrams contributing to τ → µγ decays induced by the Higgs bosons with flavor violating Yukawa couplings Y τ µ and Y µτ .
constrained by the CMS limit. However, the branching ratio H → τ µ (or A → τ µ) can still be large.

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 figure 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 appendix 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 figure 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 figure 7. Their contributions 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.. From the form of the Yukawa couplings in eq. (2.11), one can see that the pseudoscalar contribution and the scalar contributions have opposite sign. Thus in the case where m H and m A are nearly degenerate and a small mixing angle, their contributions to δa µ cancel each other. Therefore the bound from the muon magnetic dipole moment is expected to be weak, see for example the left pane of figure 5. On the other hand, if m H and m A are different, the bound from a µ could be strong as can be seen from the right pane of figure 5.

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).

τ → 3µ and
The flavor violating couplings Y τ µ and Y µτ also lead to a decay τ → 3µ as well as τ → µe + e − as shown in figure 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 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 τ − (p) → µ − (p 1 )e + (p 2 )e − (p 3 ) is 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 figure 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. 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.

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  . 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. 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 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 figure 11, while the m col distribution for the signal is given in figure 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 N back k=0 P (k; N back + N sig ) < 0.05, (4.3) 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 figure 1 and 2.
The estimated bounds for the 13 TeV LHC with 300 fb −1 luminosity are given in figure 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. 5 We ignore systematic uncertainties in this work.

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 figure 14. In the figure we also include the bounds on sin α from h → W W * measurements from the LHC run-1, [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. Figure 15a shows our projection of such a bound for the 13 TeV LHC with 300 fb −1 luminosity. In figure 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 figure 16 for sin α = ±0.05, ±0.1, ±0.5 and Y tt U = 0.05, 0.1, 1. Note that when a new decay 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 figure 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 figures 16a, 16b, 16d and 16e. Note that in figures 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 figure 14. Thus, we only consider the case Y tt U = 0.05 for sin α = 0.5, figure 16c and Y tt U = 0.05 and 0.1 for sin α = −0.5, figure 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 figure 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 figure 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 figure 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 cross-section 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. Figure 18  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 figures 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 figure 19.
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,      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 . Figure 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 ( figure 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 τ τ . 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.
The estimated bounds from the h and the H LFV searches are given in figure 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.

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 figure 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 JHEP05(2017)055 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 section 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.
x(2x−1) 2B (B.14) The 2-loop contributions to c R can be obtained from ∆c L with the replacement y ji * f,φ → y ij f,φ .
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.