Revisiting Generalized Two Higgs Doublet Model in the Light of Muon Anomaly and Lepton Flavor Violating Decays at HL-LHC

One of the main motivations to look beyond the SM is the discrepancy between the theoretical prediction and observation of anomalous magnetic moment of muon. Alleviating this tension between theory and experiment and satisfying the bounds from lepton flavor violation data simultaneously is a challenging task. In this paper, in view of generalised Two Higgs Doublet Model, we explore muon anomaly and lepton flavor violation along with constraints coming from B-physics, theoretical constraints, electroweak observables and collider data that can restrict the model parameter space significantly. We find that within the framework of generalised Two Higgs Doublet Model it is possible to obtain regions allowed by all constraints, that can provide an explanation for the observed muon anomaly and at the same time predicts interesting signatures of lepton flavor violation. Furthermore, we consider the flavor violating decay of low-mass CP-odd scalar to probe the allowed parameter space at future runs of the LHC. With simple cut-based analysis we show that part of that parameter space can be probed with significance $>5 \sigma$. We also provide Artificial Neural Network analysis which definitely improves our cut-based results significantly.

Our work is organised as follows. In section 2 we present a brief outline of the model. In section 3 we discuss the muon anomaly and its impact on our model parameter space. We then move to section 4 where we explore the allowed parameter space taking into account the theoretical constraints, constraints coming from low evergy observables and from the collider. We present a cut-based as well as neural-network-based collider analysis in section 5. We summarize our results and conclude in section 6.

The Generalised Two Higgs Doublet Model
In this section we briefly discuss the model in consideration. We follow the convention as in [17] 1 . Two scalar doublets Φ 1 and Φ 2 with hypercharge Y = 1 2 are present in this model. The most general scalar potential can be written as: (2.1) In general, M 2 12 , λ 5 , λ 6 and λ 7 can be complex while the rest of the parameters are real. However, in this work we assume CP is conserved, therefore M 2 12 , λ 5 , λ 6 and λ 7 are taken to be real. We should mention here in the absence of Z 2 symmetry (Φ 1 → Φ 1 , Φ 2 → −Φ 2 ) λ 6 and λ 7 can take non-zero values in general.
The two scalar doublets of the model can be expanded as
A key parameter of the model is tan β = v 2 v 1 . Charged Goldstones(G ± ) and physical charged Higgs state (H ± ) are produced as a linear combination of φ ± 1 and φ ± 2 . On the other hand, the same linear combination of Im[Φ 0 1 ] and Im[Φ 0 2 ] gives rise to neutral CPodd Goldstone (G 0 ) and physical CP-odd state (A). The mixing is shown in the following equations.
(2.5) 1 For general 2HDM review one should look into Ref [37] Diagonalizing the mass matrix for the CP-even neutral states we get the mass eigenstates h and H, where the states in the mass basis and in the flavor basis are related by the following rotation.
Where either h or H is assumed to behave like the Higgs of Standard Model with mass 125 GeV.
Having discussed the Higgs potential of generalized 2HDM, we proceed towards the Yukawa sector of the model. It is well-known that to avoid tree-level flavor-changing neutral current, a Z 2 symmetry is imposed on the general Yukawa Lagrangian of 2HDM. The doublets Φ 1 , Φ 2 and the fermion fields behave either even or odd under this Z 2 symmetry and depending on this behavior four common types of 2HDM, namely Type I, Type II, Lepton-specific or Type X and Flipped or Type Y are formed. In Type I, up and down type quarks and leptons couple to Φ 2 . In Type II, up-type quarks couple to Φ 2 and down-type quarks and leptons couple to Φ 1 . In Type X model up and down type quarks couple to Φ 2 and leptons couple to Φ 1 . In Type Y 2HDM, up type quarks and leptons couple to Φ 2 and the down-type quarks couple to Φ 1 .
Unlike these aforementioned types of 2HDM, in the generalized 2HDM, no Z 2 symmetry is imposed on the Lagrangian and therefore this model produces tree-level FCNC. In this case both the doublets Φ 1 and Φ 2 couple to all the leptons and quarks. In the generalized 2HDM, the Yukawa couplings to the quarks and leptons can be written as: are the Yukawa matrices whose flavor indices have been suppressed andΦ i = iσ 2 Φ * i . Expanding this equation in terms of the VEVs and physical fields, we can get the fermion mass matrix.
Without assuming any particular relation between the matrices Y 1 and Y 2 it is not possible to diagonalize the two of them simultaneously, which leads to tree-level scalar mediated FCNC. To diagonalize the fermion mass matrices by bi-unitary transformation we need two unitary matrices U f L and U f R . We adopt the convention of [38] and consider the Yukawa Lagrangian as a perturbation of Type X model in terms of FCNC couplings. Therefore we diagonalize Y u 2 , Y d 2 and Y 1 matrices where Y u 1 , Y d 1 and Y 2 remain non-diagonal leading to tree-level FCNC in the Yukawa sector. After diagonalization the Yukawa Lagrangian involving the neutral scalars can be written as follows.
, s β−α = sin(β−α) and t β = tan β. It can be easily checked that the FCNC processes occur due to the presence of non-zero Σ f matrices. When Σ f = 0 the Yukawa couplings in Eq. 2.9 reduce to the couplings in the Type X 2HDM. Following the convention of [39] Σ f matrices are parametrized in terms of dimensionless free parameters χ f s in the following manner.
In general χ f ij may not be equal to χ f ji , but for simplicity we assume χ f ij = χ f ji in our analysis.
As the rotation matrix for charged scalars is the same as that for pseudoscalars which can be seen from Eq. (2.4 and 2.5), therefore, the Yukawa couplings of the charged Higgs boson are similar to those of the CP-odd scalar and can be written as Where the sum over flavor indices is indicated, V ≡ U u L U d † L is the Cabibbo-Kobayashi-Maskawa (CKM) matrix, and P R,L = (1 ± γ 5 )/2 are the chiral projection operators. The expressions for ξ f matrices are given below.
Here also we can see that non-zero Σ matrices and equivalently non-zero ξ matrices will introduce non-trivial coupling structure even in the charged Higgs interaction with the quarks and leptons. One can check that in the absence of these matrices the couplings reduce to couplings in the Type X 2HDM.

Exploration of Muon (g − 2)
The muon anomalous magnetic moment is one of the biggest triumphs of quantum field theory. A precise measurement of this helps one to comprehend the higher order corrections contributing to it. Moreover, it indicates an existence of new physics because of the longstanding discrepancy between SM prediction and experimental observation [3]. Possibly the ongoing E989 experiment at Fermilab [4] and future E34 experiment at J-PARC [5] will shed new light on this discrepancy between the theory and experiment.
In classical quantum mechanics the value of g µ (gyromagnetic ratio for µ) is 2. It receives correction from loop effects in quantum field theory. These corrections are parametrized in terms of a µ = gµ−2 2 . In SM it receives contribution via QED, electroweak and hadronic loops. A great deal of effort has been put forth to determine the SM prediction to an unprecedented level of accuracy. SM contributions up to three orders in the electromagnetic constant, has been calculated by [40]. Taking into account pure QED, electroweak and hadronic contribution, the predicted value for muon anomaly in the SM is given by [41] a SM µ = 116591810(43) × 10 −11 (3.1) The most recent bound comes from BNL(2006) data [42] which gives a exp µ = 116592091(54)(33) × 10 −11 (3.2) The difference between the theory and experiment denotes a 3.7σ discrepancy which can certainly leave us a room to search for NP scenarios.
In this work we consider one loop as well as two loop Bar-Zee type contribution to ∆a µ in generalized 2HDM. It is shown in earlier works [43,44] that the two-loop Bar-Zee diagrams can bring sizeable contribution for a large region of parameter space. We present the diagrams for one loop contribution to ∆a µ in Fig. 1 and two loop Bar-Zee diagrams in Fig. 2. We mention here that the two loop Bar-Zee contributions dominate over the oneloop contributions. Although the two loop diagrams have a loop suppression factor but also have an enhancement factor of M 2 m 2 µ , where M is the mass of the heavy particle running in the loop namely t, b, H ± , W ± as can be seen from Fig. 2. This enhancement factor usually dominates over the loop suppression. The Bar-Zee contribution from an internal photon and heavy fermion or H ± running in the loop has been studied in great detail in the past and these diagrams give rise to major contribution to ∆a µ . The contribution from diagrams where a Z boson replaces the internal photon is negligible due to coupling as well as mass suppression. Also, the diagram involving W ± in the loop, will have negligible contribution because of suppression in the coupling between W ± bosons and non-standard CP-even Higgs in the alignment limit. We have also considered the Bar-Zee diagrams where charged Higgs replaces the neutral Higgs and W ± substitutes the internal γ in Fig. 3. It has been shown in [44] that their contribution can be sizeable in some regions of the parameter space. In Fig. 4, we also consider the diagrams with internal charged Higgs or W ± where the grey circle represents the same loops as in Fig. 3, excluding the fermion loops for the W ± diagram, because that will be a pure SM contribution. We compute ∆a µ taking into account all the aforementioned diagrams following [43,44]. Next we scan the parameter space of our model and plot allowed region in the m A − tan β plane in Fig. 5. One can check that a low mass pseudoscalar with an enhanced coupling to the τ leptons will give significant contribution to ∆a µ (see Fig. 2(top left)). In our model the coupling of pseudoscalar with a pair of τ leptons is proportional to tan β. Therefore, low m A and large tan β region is favored in the light of g µ − 2 data. While scanning the parameter space we have used the 2σ upper and lower bound on the experimentally observed central value of ∆a µ (Eq. 3.3).

Constraints on the model
We have seen from the discussion in the previous section that the major contribution to the anomalous magnetic moment of muon comes from the low mass pseudoscalar contribution at moderate to large tan β. However, in the presence of the non-diagonal Yukawa matrices it is inevitable that the similar contributions will also generate FCNC processes. This flavor changing processes include loop induced µ → eγ, τ → eγ, τ → µγ, µ → 3e and µ − e conversion, all of which put a strong constraint on the flavor changing couplings as well as the (pseudo)scalar masses and tan β. It is evident that low mass pseudoscalar and large tan β regions will be disfavored from the limits coming from the low energy measurements, which seems to be in tension with the requirement of muon (g − 2). We study these constraints carefully in the upcoming subsection and explore the regions of parameter space consistent with the limits from the non-observation of low energy flavor violating processes as well as the experimental observation of (g µ − 2). Moreover, these flavor changing vertices also give rice to FCNC in the (pseudo)scalar mediated tree-level decays in the leptonic final states. Our main objective of this work is to probe this region of parameter space in the collider. We mention here that to get sufficient event rate at the collider, we focus on the low mass range of the decaying (pseudo)scalar. We now proceed to discuss various constraints on our model which further guide us to choose our benchmarks for the upcoming direct search analysis at the collider.

Limits from low energy measurements
In the SM, lepton flavor is conserved since neutrinos are massless. In neutrino oscillation [6,7], LFV has been observed in the neutrino sector. However, till date LFV has not been observed in the charged lepton sector. Therefore, lepton flavor violation can be treated as one of the important tools to search for new physics. Many new physics models can accommodate LFV processes. Since, no such signal has been observed yet, there are strong limits on these LFV processes [11].
The recent bound on BR(µ → eγ) < 4.2×10 −13 comes from MEG experiment [12]. The other important constraint in LFV will come from µ → 3e which is a natural consequence of µ → eγ decay when the resulting photon converts to e + e − pair. Apart from that, µ − e conversion in nuclei can also be an important signature of LFV. The relation between µ → eγ and other possible LFV constraints, namely, BR(µ → 3e) and µ−e conversion(CR)  are given by [10,13] Figure 5. The allowed region in m A − tan β plane from g µ − 2 data at 2σ.
If we try to translate the limits according to the relations in Eq. 4.1, we find that the limit on CR(µ T i → e T i) has to be < 2.1 × 10 −15 for it to be of same strength as the limit from µ → eγ, while from experiment [14] this upper limit is 6.1 × 10 −13 [9]. Similar argument holds for BR(µ → 3e) branching ratio. From Eq. 4.1, BR(µ → 3e) should be < 2.62 × 10 −15 , for this upper limit to be of same strength as BR(µ → eγ) while the experimental upper bound is 1.0 × 10 −12 [8]. Therefore, it is evident that these two constraints are relatively weak. Hence, for our analysis we take into account the strongest bound which is coming from µ → eγ .
Similar to the µ − e sector, there are strong constraints on (τ → eγ) and (τ → µγ) branching ratio. Current Bound on BR(τ → eγ) < 3.3 × 10 −8 [11] and BR(τ → µγ) < 4.4×10 −8 [11] puts a strong constraint on y τ e and y µτ respectively. One should also take into account the limits on the BR(τ → 3e) < 2.7 × 10 −8 [15] and BR(τ → 3µ) < 2.1 × 10 −8 [15]. However, compared to BR(τ → eγ), the limit on BR(τ → 3e) is weaker due to an addition suppression of factor α. Same is true for the limit on BR(τ → 3µ) which puts much weaker limit compared to BR(τ → µγ). We calculate the LFV processes in one-loop as suggested by Ref. [13]. The presence of flavor non-diagonal Yukawa matrices in the generalized 2HDM gives rise to flavor violating coupling between scalars and leptons at the tree-level. This in turn enables the LFV decays shown in Fig. 6 at one loop. It is worth mentioning that the loop contribution from the neutral scalars dominates over the contribution from the charged Higgs loop. Among the neutral scalars primarily the pseudoscalar loop contributes to this process i → j γ owing to low pseudoscalar mass.
It is interesting to note that BR(τ → eγ) and BR(τ → µγ) constrain the couplings y τ e and y µτ respectively. The reason behind this is the following. The major contribution to the corresponding amplitudes come from the τ loop (see Fig. 6 (left)), where y τ e and y µτ appear with m τ . Therefore these terms dominate over the other terms which are accompanied by m µ , m e or are product of two LFV couplings. Hence the upper limit on the aforementioned branching ratios constrains particularly y τ e and y µτ . However, the situation is different in case of BR(µ → eγ). In this case the τ loop which has the highest contribution in terms of loop integral comes with a coefficient y µτ × y τ e . Therefore its contribution can be comparable with the e or µ loop, with coefficients y µe multiplied with m e or m µ . Therefore BR(µ → eγ) is not solely controlled by y µe .
We have seen that BR(τ → eγ) constrains y τ e < 10 −4 and BR(τ → µγ ) constrains y τ µ < 10 −4 . However, for y µe the situation is not so straightforward. Unlike τ → eγ and τ → µγ, the decay µ → eγ does not primarily constrain y µe coupling. With suitable choice of y µτ and y τ e couplings it is therefore possible to satisfy BR(µ → eγ) with a wide range of y µe coupling. However, to satisfy all the three conditions simultaneously, where y µτ and y τ e couplings are already strongly constrained by BR(τ → eγ) and BR(τ → µγ), the coupling y µe also gets strong upper bound (< 10 −5 ).
It is important to note that the branching ratios we just described also depend strongly on the scalar masses and tan β, along with the flavor changing couplings. In Fig. 7-10 we have plotted the regions allowed by LFV constraints in m A −tan β plane for specific choices of flavor changing Yukawa couplings. In Fig. 7-10, we have also superimposed the region allowed by (g µ − 2) data on the region allowed by low energy LFV data.
It can be clear from Fig 7,8 that the two limits tend to favor non-overlapping regions, unless the LFV Yukawa couplings are below certain values. In Fig. 9 and 10, we present our choice of LFV Yukawa couplings for which we get the an overlapping region that is allowed by (g µ − 2) as well as low energy LFV constraints. We specifically concentrate on the scenario depicted in Fig. 10, because the values of flavor violating Yukawa couplings in this case produce adequate event rate which can be probed at the HL-LHC. Therefore this region is of primary interest to us from the collider point of view.

Theoretical constraints
The constraints from the requirements of vacuum stability and perturbativity have been studied in detail in earlier works [45,46]. It has been pointed out that large separation between the m A and m ± H is disfavored from the theoretical considerations of vacuum stability and perturbativity. Since we are interested in the low mass pseudoscalars from the requirements of (g µ − 2), it is imperative to check the upper limit on m ± H compatible with low m A . The vacuum stability and perturbativity conditions put bounds on the λ parameters and thereby correlate the masses of different neutral and charged scalars. The vacuum stability condition requires [46,47] The resulting squared-masses for the CP-odd and charged Higgs states are given by [46] It is clear from 4.4, the mass-square difference m ± H 2 − m 2 A is proportional to λ 5 − λ 4 , which should be less than λ 3 + √ λ 1 λ 2 . Along with the vacuum stability criteria, the requirement of perturbativity, ie. all the quartic couplings C H i H j H k H l < 4π puts an upper limit on m ± H − m A . In the scenario when m h = 125 GeV, the parameter space allowed by stability and perturbativity constraints are indicated in [45].
In Fig. 11 we show the parameter space allowed by stability, unitarity and perturbativity constraints. We show only the low m A region, we will be interested in this region in our collider analysis. In Fig. 11(left) we show the upper limit of m H ± as a function of m A as pointed out in the foregoing discussion. We can see that m H ± < 170 − 180 GeV is allowed for low m A . In Fig. 11(right) we show that constraints in the tan β − m A plane. We see that although very large tan β is allowed from perturbativity considerations, low to moderate tan β values are much more favored compared to the high values.
It is important to notice that the low m A region of parameter space which we are interested in, yield a substantial branching fraction for h → AA decay, where h is the 125 GeV Higgs and m A < m h 2 . The experimental upper limit on this branching ratio is rather strong [48], where a stringent limit comes from the search for (pp → h → AA) process in the µ + µ − τ + τ − final state. The only way such a small branching ratio can be achieved is when the coupling g hAA is extremely small. This in turn imposes a relation between m 2 12 , tan β and m A [49]. However m 2 12 is a parameter a crucial parameter which ensures perturbativity. It is required for perturbativity that m 2 12 ∼ m 2 H tan β . It is shown in [49], in the case where 125 GeV Higgs is the lightest Higgs boson, and m H > 125 GeV, it is possible to obey the perturbativity constraints as well as the upper limit on BR(h → AA) for low tan β < 10 and the mass gap m H − m h is not very large. Although this region is phenomenologically viable, the (g µ − 2) requirements(see Fig. 5) impose that m A should also be very small, ie m A < 10 GeV. This statement is only valid in the 'right-sign' region of 2HDM where Higgs coupling with the fermions and gauge bosons are of same sign. In the so-called 'wrong-sign' region where Higgs coupling to the fermions and gauge bosons Figure 10. The magenta, green and cyan regions are the allowed range for µ → eγ, τ → eγ and τ → µγ respectively. The blue band is the allowed 2σ allowed range for muon anomaly. The overlapping regions satisfy both constraints. The flavor changing couplings are taken to be y µe = 2 × 10 −7 , y τ e = 5 × 10 −5 , y µτ = 5 × 10 −5 . are of opposite sign, gives rise to entirely different allowed region and phenomenological signatures. We have not explored this scenario in the current work and will be pursuing in a future study. The other possibility is to consider the case when the heavier CP even Higgs is SM-like, ie m H = 125 GeV. However in this case the LEP limit implies either m A or m h can be < m H 2 . We consider the low mass pseudoscalar, and therefore m h > m H 2 . Here also, like the previous case the limit on BR(h → AA) will indicate extremely small value of the coupling g HAA whose expression is given as follows: In this case there is more freedom compared to the previous case, in terms of the allowed parameter space. One can have a pseudoscalar mass > 10 GeV with moderate tan β, with suitable value of m 2 12 and m h , while satisfying perturbativity condition and the small BR(H → AA) simultaneously. This point onwards, we will explore this particular scenario, ie. for our work m H = 125 GeV.

Electroweak constraints
In this section we analyze the impact of constraints arising from electroweak precision measurements, especially the oblique parameters [50,51] on our model. The experimental collaboration Gfitter group [52] has published a contour in the plane of S and T parameter taking into account the correlation between them. The status of 2HDM in the light of the most recent global electroweak data has been presented in [53]. We mention here that we have used the elliptic contour which has been computed with U as a free parameter. This choice leaves us with a less constrained parameter space compared to the cases U = 0.
We have calculated the oblique parameters in our model and obtained the allowed region of parameter space at 3σ. Here also we concentrate on low m A region and considered the case m H = 125 GeV, ie. the second lightest CP-even Higgs is SM-like. m H ± has been varied from 90 GeV -200 GeV. In Fig. 12, We present the allowed region in the plane of m A and (m h − m H ± ). It is evident from the figure that as the pseudoscalar mass decreases, the mass difference m h − m H ± should also decrease to obey the constraints from oblique parameters.

Constraints from B-physics
From the discussions of Section 2, it is clear that in presence of flavor changing couplings in the Yukawa sector the charged Higgs couplings to quarks and leptons are also modified. With new free parameters in the Lagrangian, interesting phenomenologies are likely to show up in rare decay processes that were suppressed in the SM. One such possibility in the rare processes involving B−meson. The free parameters of the model gets constrained by the experimental upper bounds on such rare FCNC processes. It is clear from Eq. 2.10 that the FCNC within the first two generations are naturally suppressed by the small quark masses, while a larger freedom is allowed for the FCNC in the top and bottom quark sector. In our analysis also we have taken only λ tt and λ bb to be non-zero where λ tt and λ bb are the htt and hbb coupling strengths respectively, considering h to be the non-SM like CP-even Higgs.
The rare processes involving B-mesons primarily include the decay B → X s γ, B s → µ + µ − , B ± → τ ± ν τ , B 0 −B 0 mixing whose strength is determined by the mass difference ∆M B between the two states. The most stringent constraint comes from the B → X s γ decay. The impact of these constraints in terms of specific types of 2HDM have been studied in great detail in earlier works [54,55]. In conventional type I and type II 2HDM, the dominant additional contribution to the loop induced decay B → X s γ comes from the charged Higgs boson-top quark penguin diagrams and its contribution depends on m H ± . In type II 2HDM, this extra contribution interferes constructively with its SM counterpart and therefore the lower bound on the charged Higgs boson becomes rather high (m ± H > ∼ 600 GeV). In type I, the charged Higgs penguin diagram's contribution interferes destructively with its SM counterpart and gives negligible result at large tan β. Therefore no strong constraint appear on the mass of the charged Higgs in type I model. The type X model has same structure as type I, as far as the interactions of Higgs with the quark sector is concerned. Therefore Type X models also do not receive any strong lower bound on m H ± . As we can think of our model as a perturbation from the type X scenario, in the absence of the extra terms in the Yukawa Lagrangian, there is no strict lower bound on the charged Higgs mass. However, even in the presence of non-zero FCNC Yukawa matrix elements, it is possible to have low enough charged Higgs mass [20,31,34] with suitable choice of λ tt and λ bb couplings. We have taken in our analysis λ tt ∼ 0.5 and λ bb ∼ 2, which allows a charged Higgs mass m H ± > ∼ 150 GeV.
Another decay process which can constrain our model parameters space is B ± → τ ± ν τ where charged Higgs enters at the tree level itself. The observed branching ratio for the process B ± u → τ ± ν τ = (1.06 ± 0.19) × 10 −4 [56]. The decay B ± c → τ ± ν τ , although has not been observed, but puts an upper limit (< 30%) [56] on the branching ratio for this decay. However, we have assumed only λ tt and λ bb are non-zero in the quark sector, we find out that these limits essentially reduces to a limit on λ bb and tan β. In [34], it has been shown that λ bb ∼ 2 is favored for large or moderate tan β. The constraint from ∆M B puts an an upper limit on λ tt as a function of the charged Higgs mass [20]. m H ± > ∼ 150 GeV is allowed for λ tt < ∼ 0.5. Therefore our choice of parameter space obeys this constraint as well.

Constraints from direct search at the colliders
Constraints can be obtained from collider searches for the production and decay of on-shell neutral and charged Higgs bosons. There have been numerous searches in the past in this direction. The LEP experiments have looked for pair production of charged Higgs bosons in the process e + e − → γ/Z → H + H − . In this process all the couplings that appear are essentially gauge couplings, this predictions for this process therefore depends only on the charged Higgs mass m H ± . However the decay and branching fractions of H ± are indeed model dependent. But a combined search for H ± in τ ν and cs channel put a robust lower limit of 80 GeV on m H ± [58]. This limit only mildly depends on BR(H ± → τ ν).
At the LHC the charged Higgs search can be categorized in two types. For m ± H < m t , charged Higgs can be produced from the decay of top quark(t → bH ± ). This decay has been searched for in τ ν [59,60]  Collider searches for the non-standard neutral Higgs also put constraints on the parameter space of interest. Searches for non-standard Higgs bosons are performed at the LHC in various channels with SM final states. As we are specifically interested in the low pseudoscalar mass region with enhanced coupling to leptons, the limits which are crucial for our case comes from the search for low pseudoscalar produced in association with b quarks and decaying into τ τ final state [67,68]. Constraints from the search for low mass (pseudo)scalar produced in association with bb and decaying into bb [69, 70] has been taken into account. CMS has also searched for decay involving two non-standard Higgs bosons such as h/H → Z( )A(τ τ ) [71] and h/H → Z( )A(bb) [72]. However these limits become applicable for heavier CP-even Higgs > ∼ 200 GeV. Therefore these particular searches do not have any considerable affect on our parameter space.
We mention here that one should also take into account the limits coming from the direct search of the 125 GeV Higgs in various final states including τ τ [73,74],µµ [75,76]. Moreover, as the focus of our work is FCNC in the Yukawa sector, the constraints coming from flavor violating decays of 125 GeV Higgs boson also put constraints on the flavorviolating Yukawa matrix elements. The 125 GeV Higgs decaying to eµ and eτ final state have been looked for by the CMS experiments [77]. CMS also puts an upper limit on the branching ratio for 125 GeV Higgs decaying to µτ final state [78]. Undoubtedly, these limits are crucial for our study. However, as we strictly confine ourselves to alignment limit (cos(β − α) ≈ 0.999), the flavor violating decays of the 125 GeV Higgs(H in our case) will receive a suppression by a factor sin 2 (β − α) which can be seen from Eq. 2.9. Therefore in this limit the constraints coming from lepton flavor violating decays of the 125 GeV Higgs are trivially satisfied.
An important constraint comes from the direct search for 125 GeV Higgs decaying into two light pseudoscalar final states when it is kinematically allowed. The upper bound on this branching ratio translates into severe constraint on the parameter space of this model. We have discussed this in detail in a previous subsection 4.2 and have taken into account in our analysis.

Collider Searches
From our discussions in the previous sections it is clear that the existence of flavor violation in the lepton Yukawa sector gives rise to flavor-violating decays of µ and τ leptons. The presence of off-diagonal elements in the Yukawa matrices are the source of the lepton flavor violation in generalized 2HDM. The flavor violating decays of leptons are introduced at loop level via flavor violating coupling between the scalar and leptons at tree level. These flavor non-diagonal tree-level Yukawa coupling between the scalar and leptons will also give rise to interesting phenomenology at the colliders [16][17][18][19].
In this work we consider probing the CP-odd scalar A in flavor violating leptonic decay mode in generalized 2HDM at the HL-LHC. Our signal process is given as where = e, µ and τ denotes the leptonic decay of τ . The signal of our interest is + − + / E T . The SM processes that can give rise to similar final states are τ τ /ee/µµ, tt, W ± +jets, Diboson, SM Higgs [16,79]. Out of these backgrounds, the major background in our signal region is the τ τ . Due to large cross-section, tt also serves as the other major background.
For our analysis, we choose three benchmark points valid from all the experimental and theoretical constraints and quote their production cross-section in Table 1. We mention here that since the branching ratios of the pseudoscalar decaying to flavor violating final states is very small (BR(A → µτ ) ≈ BR(A → τ e) ≈ 10 −7 ), owing to the smallness of lepton flavor violating Yukawa couplings, we are compelled to choose low mass pseudoscalar which will have considerable production cross-section and therefore will be a viable candidate to search for in the collider.
We first present the cut-based analysis for this channel in the following subsection. Next we explore the possible improvement of our results with multivariate techniques using Artificial Neural Network (ANN).

Cut-based Analysis
The signal events are generated in Madgraph5@NLO [80] implementing the model file in FeynRules [81]. We generate both signal and SM backgrounds events at leading order (LO) in Madgraph5@NLO [80] using the NNPDF3.0 parton distributions [82]. The parton showering and hadronization are done using the built-in Pythia [83] within Madgraph. The showered events are then passed through Delphes(v3) [84] for the detector simulation where the jets are constructed using the anti-K T jet algorithm with minimum jet formation radius ∆R = 0.5. The isolated leptons are considered to be separated from the jets and other leptons by ∆R i 0.4, i = j, .
To generate our signal and background events, we employ the following pre-selection cuts.
The b-jets are tagged with the p T -dependent b-tag efficiency following the criteria of Ref. [85] which has an average 75% tagging efficiency of the b-jets with 50 GeV < p T < 200 GeV and 1% mis-tagging efficiency for light jets.
Additionally, we propose the following selection cuts on certain kinematic observables to disentangle the signal from the SM backgrounds that would enhance the signal significance. We describe those observable in the following.  • p T of the leptons: In Fig. 13, we present the transverse momentum p T of the leading and sub-leading leptons. For the signal, the leptons coming from the decay of a low mass pseudoscalar, tend to have low p T . Since the distributions are mostly overlapping for both signal and τ τ background, it is very difficult to put any hard p T cut. However, to affirm that our signal has 2 isolated leptons, we reject any third lepton with p T ( ) > 10 GeV. Moreover, since our signal is hadronically quiet, we put a jet-veto of with p T (j) > 20 GeV. We also reject any b -jet with p T (b) > 20 GeV. This helps us reduce the tt semileptonic and W ± + jets background. • Selecting / E T : For signal the neutrino is coming from the leptonic decay of τ . The τ comes from the decay of low mass pseudoscalar. Therefore for signal the / E T tends to be small. For τ τ background, the neutrinos are produced almost back to back. So, the lower / E T bins are populated both for signal and τ τ background. On the other hand, top decay being a three-body decay, the / E t produced in tt event peaks at a larger value. We do not put any additional cut on / E T . However for sake of completeness, we present the distribution in Fig. 14(left).
• Invariant mass of the di-lepton pair: In Fig. 14(right) we show the invariant mass of the dilepton pair. In case of signal, the leptons are coming from the decay of a low mass pseudoscalar, therefore its distribution peaks at a lower value, compared to τ τ and tt background. This variable plays a crucial role to discriminate between signal and background. We use this variable for our cut-based analysis.
• The collinear mass: The collinear mass is defined as follows: with the visible momentum fraction of the τ decay products being, and M vis is the visible mass of the τ − system. The variable M collinear essentially reconstructs the mass of the pseudoscalar. From Fig. 15 (left) it is clear that M collinear yields a clear distinction between the signal and the irreducible τ τ background. A suitable cut should be imposed on this variable to reduce the τ τ background.
• The transverse mass variable: The transverse mass is defined as where ∆φ − / E T denotes the azimuthal angle between the leading lepton and / E T . From Fig. 15 (right) it is clear that a cut on M T variable helps us reduce the tt background.
• Angle between the lepton: The angle between two leptons is strictly correlated to the invariant mass. Since for the signal the invariant mass of the dilepton pair is small, the azimuthal angle between the two leptons ∆Φ( + , − ) appears at lower value compared to the τ τ background where the leptons are produced back to back and as a result ∆Φ( + , − ) peaks around π. It is clear from Fig. 16 that a suitable cut on this observable will enhance the signal-background separation. With optimized cuts on the aforementioned variables (listed in Table 2), we get the signal significance as presented in Table 2. The significance [86] is calculated using the following formula. S = 2[(S + B)ln(1 + S B ) − S] where S and B denotes the number of signal and background after applying all the cuts. We mention here that, we multiply the signal and background cross-sections with respective k-factors to take into account the next-to-leading-order (NLO) effects. For signal, we use the k-factor of 2 [87] while for tt and τ τ background, we take the k-factor to be 1.6 [88] and 1.15 [89] respectively.

Improved analysis with Artificial Neural Network (ANN)
Having performed the cut-based analysis, we proceed to analyze the same channel(dileptons + / E T ) with ANN [90] at the LHC. We explore the possibility of improvement in our results. This method has been used extensively in the recent past and it has proved to improve the results of cut-based analysis multifold enabling better separation between the signal and backgrounds. Significant work has been done in the context of Higgs sector [91][92][93][94][95]. In collider search for dark matter also this methods have been proved to be extremely useful [94,96]. Therefore we employ this tool in our present analysis also where signal yield is small because of minuscule branching fraction (∼ 10 −7 ) of the flavor changing decay of the pseudoscalar and the distinction between signal and background becomes crucial. We have examined and computed the maximum signal significance for the benchmarks that we considered, achievable at the HL-LHC using these technique. The toolkit used for ANN analysis is a python-based deep-learning library Keras [97].
From our analysis in the previous section we identify the important input variable that provide substantial distinction between signal and backgrounds. We mention here that the choice of input variables play a crucial role. In Table. 3 we present all the input variables that we have used for our analysis.  For ANN analysis we have chosen a network with four hidden layers with activation curve relu at all of them. The batch-size is taken to be 1000 and the number of epochs is 100 in our case for each batch. We have used 80% of the dataset for training and 20% for validation. One possible demerit of these techniques is over-training of the data sample. In case of over-training the training sample indeed gives extremely good accuracy but the validation sample fails to achieve the same degree of accuracy as that of the training sample. However we have explicitly checked that with our choice of network parameters the algorithm does not over-train.
We find that the variables M , M collinear , M T , ∆φ and ∆R play crucial role in signal-background separation. However, there is a strong correlation between ∆R , ∆φ and M as we have already discussed in the previous subsection. We mention here that in order to obtain a better performance from the network we have applied two basic cuts, namely M < 30 GeV and M collinear < 40 GeV on signal and background events over and above the lepton selection and jet-veto. From our discussion of the cut-based analysis we know that these cuts guide us towards the so-called signal region. Therefore the network will be trained better to separate signal from background specifically in the signal region, this results in a better accuracy in the output. The accuracy we get is 99%(BP1), 98%(BP2) and 96%(BP3) which indicates very good discriminating power between signal and background. We present in Fig. 17, the Receiver Operating Characteristic (ROC) curve for the three benchmark points of our choice.
The area under curve is 0.999(BP1), 0.998(BP2) and 0.994(BP3). We plot only a part of the ROC curve which is relevant for our analysis. We scan over the ROC curve and choose suitable points which yield the maximum signal significance. In Table. 4, we present the signal significance S for all the signal benchmarks we have considered. 10.9 σ BP2 7.5 σ BP3 5.3 σ Table 4. Signal significance for the benchmark points at 14 TeV with L = 3 ab −1 with cuts+ANN.
Comparing the results of ANN in Table. 4 and the cut-based results in Table. 2 we can see that our analysis with ANN improves the results of cut-based analysis significantly.

Conclusion
The motivation behind this work is a much-desired reconciliation between the observed muon anomaly and LFV constraints. In this regard, we consider an extension of the SM with extended scalar sector, namely, generalized 2HDM. The additional non-standard scalars of this model take part in the muon anomaly and flavor non-diagonal Yukawa matrices lead to LFV processes. We show that the long-standing problem of muon anomaly and LFV constraints can be solved simultaneously over considerable range of parameter space in this model. We show such a region in Fig. 10 with flavor changing couplings fixed at y µe = 2 × 10 −7 , y τ e = 5 × 10 −5 and y µτ = 5 × 10 −5 where both muon anomaly and LFV constraints are satisfied.
We then proceed to implement theoretical constraints pertaining to the requirement of perturbativity, unitarity and vacuum stability. The flavor non-diagonal Yukawa matrices also get severely constrained by the B-physics observables. The direct searches for the SM Higgs as well as the additional scalar states in the collider put another set of bounds on the model parameter space. One such crucial direct search constraint turns out to be the search for the SM Higgs decaying to two light pseudoscalars. Our aim in this work is to search for lepton flavor violation in the scalar sector at the collider. Therefore the scalar states with low mass, prove to be the best candidate for such searches owing to their large production cross-section. We have also found that it is the light CP-odd scalar(A) of our model that helps us explain the (g µ − 2) data. The lightness of the pseudoscalar, also implies a large branching ratio of the 125 GeV Higgs into a pair of pseudoscalars when the decay is kinematically feasible. To ensure the upper bound to this branching fraction coming from collider data, along with the perturbativity requirements, one is drawn to the scenario where the observed 125 GeV Higgs is the heavier of the two CP-even states of 2HDM in the alignment limit. However, this statement is valid only in the 'right-sign' region of 2HDM which we have considered in this work. The 'wrong-sign' scenario will give rise to a different allowed region and interesting phenomenology, which we want to pursue in detail in the future.
After finding out the region allowed by all constraints, we look for flavor violating decay of CP-odd scalar (A) in the τ → + − + / E T final state, where τ decays leptonically and , = e, µ. Performing a rectangular cut-based analysis for 14 TeV LHC with 3ab −1 luminosity, we show that the significance drops from ∼ 10σ to ∼ 1σ as the mass of the scalar increases from 21 GeV to 27 GeV. We then perform an analysis using ANN and observe significant improvements of our results from the cut-based analysis. We would like to point out that although we have probed a narrow region of parameter space in terms of pseudoscalar mass, we did it in order to investigate the reach of LFV searches at 14 TeV LHC with 3 ab −1 at > ∼ 5σ significance. The results of ANN analysis in Table. 4 indicates that even higher masses of A can be probed, although with somewhat lower significance (< 5σ). Also, a further upgrade in the luminosity as well as energy frontier will enable us to probe heavier CP-odd scalar decaying into lepton flavor violating final states. [64] ATLAS collaboration, T. A. collaboration, Search for charged Higgs bosons in the τ +jets final state using 14.7 f b −1 of pp collision data recorded at of pp collision data recorded at √ s=13 TeV with the ATLAS experiment,ATLAS-CONF-2016-088, .
[66] ATLAS collaboration, T. A. collaboration, Search for charged Higgs bosons in the H ± → tb decay channel in pp collisions at √ s = 13 TeV using the ATLAS detector, ATLAS-CONF-2016-089, .