Higgs Lepton Flavor Violating Decays in Two Higgs Doublet Models

The discovery of a non-zero rate for a lepton flavor violating decay mode of the Higgs boson would definitely be an indication of New Physics. We review the prospects for such signal in Two Higgs Doublet Models, in particular for Higgs boson decays into $\tau \mu$ final states. We will show that this scenario contains all the necessary ingredients to provide large flavor violating rates and still be compatible with the stringent limits from direct searches and low-energy flavor experiments.


INTRODUCTION
The discovery of the Higgs boson in 2012 at the LHC by the ATLAS and CMS collaborations constitutes a historical milestone for particle physics and another brilliant triumph for the Standard Model (SM). With this long-awaited completion of the SM particle spectrum, it stands as one of the most successful theories ever built, providing precise predictions for a wide range of particle physics phenomena, in good agreement with a large amount of experimental results at low and high energies.
Despite this remarkable success, many fundamental questions remain unanswered in the SM. The list is long and contains experimental observations that cannot be addressed in the SM and theoretical issues that cannot be fully understood in its context. It includes the origin of neutrino masses, the nature of dark matter, the conservation of CP in the strong interactions or the reason for the replication of fermion generations, to mention a few. These open problems clearly call for an extension of the SM with new states, presumably present at high energies, and/or new dynamics.
New Physics (NP) may manifest in the form of Higgs boson properties different from those predicted by the SM. For this reason, it is crucial to look for deviations in the Higgs couplings to fermions and gauge bosons and the Higgs boson total decay width or the existence of new Higgs decay channels. In particular, new degrees of freedom coupling to the SM leptons and the Higgs boson could induce non-zero lepton flavor violating (LFV) Higgs boson decays, such as h → ℓ i ℓ j with i = j, indeed a common feature in many models with extended scalar or lepton sectors. The observation of these processes, strictly forbidden in the SM, would provide a clear hint of NP at work.
Higgs lepton flavor violating (HLFV) signatures face many indirect constraints, since most of the NP scenarios that lead to them modify other Higgs properties as well, some already experimentally determined to lie close to the SM prediction. Moreover, HLFV signatures typically come along with other LFV processes, such as the ℓ i → ℓ j γ radiative decays. While model-independent studies [1][2][3][4][5][6] have shown that large HLFV rates are in principle compatible with the existing experimental constraints, this is not generally the case in specific models. In fact, most models predict HLFV rates below the current LHC sensitivity. In contrast, the Two Higgs Doublet Model (2HDM) [7,8] has been shown to be able to accommodate large h → ℓ i ℓ j branching ratios [9], clearly within the reach of the ATLAS and CMS detectors.
The rest of the manuscript is structured as follows. In section 2 the current experimental bounds for the rates of several LFV processes, including h → ℓ i ℓ j , are briefly discussed. Section 3 contains general model-independent considerations, while section 4 argues that a type-III 2HDM may accommodate sizable h → ℓ i ℓ j branching ratios and makes this point explicitly after introducing our 2HDM notation and conventions. Finally, we comment on some 2HDM scenarios that generate neutrino masses and may lead to large HLFV rates in section 5 and conclude in section 6.

EXPERIMENTAL STATUS
A remarkable effort has been devoted to the search for LFV in processes involving charged leptons, resulting in impressive bounds in some channels. This has been particularly well motivated after the discovery of neutrino flavor oscillations, which imply that charged lepton flavor violating processes must exist, although perhaps with low rates.
In what concerns HLFV, the ATLAS and CMS collaborations have searched for h → ℓ i ℓ j , setting limits for the corresponding branching ratios in the 10 −5 − 10 −3 ballpark, as shown in Table 1 1 . The CMS limit on BR(h → µe) has been obtained using √ s = 8 TeV data, whereas the rest of ATLAS and CMS limits have been updated including √ s = 13 TeV data. Dedicated strategies can in principle improve these limits substantially with future LHC data [3], in particular in the 14 TeV HL-LHC phase [63]. Furthermore, HLFV can also be searched for at e + e − colliders, which offer a very clean environment, perfectly suited for the exploration of the Higgs boson properties. As shown 1 We ∼ 10 −9 [99] µ − , Ti → e − , Ti 4.3 × 10 −12 [103] ∼ 10 −18 [104] µ − , Au → e − , Au 7 × 10 −13 [105] µ − , Al → e − , Al 10 −15 − 10 −18 [106] µ − , SiC → e − , SiC 10 −14 [107] by several analyses [63,66,86,87], the planned future e + e − facilities (CEPC, FCC-ee, and ILC) can probe HLFV branching ratios as low as 10 −5 − 10 −4 , improving on the current LHC limits by about one order of magnitude for channels involving the τ lepton. The NP degrees of freedom and interactions leading to h → ℓ i ℓ j also generate other LFV processes, such as ℓ i → ℓ j γ . Since these are subject to much stronger experimental bounds, they tend to be crucial constraints in most specific models and must be considered in any phenomenological study on HLFV decays. Table 2 collects the current bounds and future sensitivities for several LFV processes of interest. Muon LFV observables have the best experimental limits due to the existence of high-intensity muon beams, while the branching ratios for tau LFV decays are bound to be below ∼ 10 −8 . The most constraining processes in many models are the ℓ i → ℓ j γ radiative decays. The MEG collaboration has established the strong limit BR(µ → eγ ) < 4.2 · 10 −13 , a bound that will be improved by about an order of magnitude in the MEG-II phase [97]. Regarding the ℓ i → ℓ j ℓ k ℓ k 3-body decays, the µ → eee branching ratio sensitivity is expected to be improved by four orders of magnitude by the Mu3e experiment [101]. Finally, the most spectacular progress in the search for LFV are expected in µ − e conversion experiments, which are also expected to improve the current limits for different nuclei by several orders of magnitude [104,106,107]. See Calibbi and Signorelli [108] for an experimental and theoretical review of the current situation in charged lepton flavor violation experiments.

MODEL-INDEPENDENT CONSIDERATIONS
In order to explore HLFV in a model-independent way, it proves convenient to adopt an approach based on Effectively Field Theory (EFT). This is particularly well motivated due to lack of NP signals at the LHC, which arguably implies that any new particles responsible for HLFV would lie clearly above the electroweak scale. We will now continue along the lines of Herrero-Garcia et al. [72] 2 .
In addition to canonical kinetic terms, the SM Lagrangian contains the Yukawa terms for the leptons where ℓ = ν e L ∼ 1, 2, − 1 2 , e ≡ e R ∼ (1, 1, −1) and denote the SM lepton SU(2) L doublets and singlets and Higgs doublet, respectively, and we give the representation under the SM gauge group SU(3) c × SU(2) L × U(1) Y . Ŵ e is a 3 × 3 complex matrix that can be taken to be diagonal without loss of generality. Therefore, the three lepton flavors are exactly conserved in the SM, which then possesses a G f = U(1) e × U(1) µ × U(1) τ global flavor symmetry 3 .
The flavor symmetry G f gets generally broken in the presence of NP. This can be generically parametrized by means of nonrenormalizable gauge-invariant operators of dimension d > 4 that encode the LFV effects induced by unknown heavy states, Here is the scale of NP, of the order of the masses of the states whose decoupling induces the dimension-d operator Q i , and C i the associated Wilson coefficient. There are many of such LFV operators. However, the only dimension-six operator giving rise to Higgs LFV is Q eϕ , defined as This operator was first highlighted in the context of HLFV in Harnik et al. [3], later denoted the Yukawa operator in Herrero-Garcia et al. [72] and is one of the operators in the Warsaw 2 See also the comprehensive reference Doršner et al. [33] for a similar reasoning in a multi-Higgs EFT that further strengthens the case for potentially large HLFV effects in the 2HDM. 3 The global flavor symmetry of the SM is known to be broken since the experimental observation of neutrino flavor oscillations. Therefore, the NP behind the generation of neutrino masses must necessarily violate G f and induce LFV processes such as ℓ i → ℓ j γ and h → ℓ i ℓ j . The resulting rates in some specific models are too low to be observed in any foreseeable experiment. For instance, in the SM minimally extended with Dirac neutrino masses one expects tiny LFV branching ratios, as low as BR(µ → eγ ) ∼ 10 −55 [109] or BR(h → τ µ) ∼ 10 −56 [72]. However, this is not a generic expectation, since these rates get hugely enhanced in most NP scenarios. We refer to section 5 for more details about the connection between HLFV and neutrino masses.
basis of the Standard Model Effective Field Theory [110] 4 . Any additional gauge-invariant dimension-six operator leading to HLFV can be shown to be redundant, and therefore reducible to Q eϕ by using equations of motion, Fierz transformations or other field redefinitions [3]. Therefore, all HLFV effects are encoded (at least in scenarios leading to NP contributions of dimension six) by Q eϕ . After electroweak symmetry breaking, the SM Yukawa term in Equation (1) and the NP contribution encoded in Q eϕ add up to where is the 3 × 3 charged lepton mass matrix and are the Higgs boson couplings to a pair of charged leptons. We have used the decomposition ϕ 0 = 1 the Higgs boson couplings to charged leptons read While the first term in Equation (8) is proportional to the charged lepton masses, the second one can in general contain off-diagonal entries and induce HLFV processes. In particular, this piece leads to the h → ℓ i ℓ j decays, with i = j. One finds where Ŵ h (≃ 4 MeV in the SM) is the total Higgs boson decay width.
So far, we only discussed the Q eϕ operator, which induces the HLFV decays h → ℓ i ℓ j we are interested in. However, in a complete ultraviolet theory other operators will be generated as well. In particular, operators that give rise to other LFV processes, such as ℓ i → ℓ j γ , with much stronger experimental bounds. At dimension six, two gauge-invariant operators of this type exist [110], where τ I , with I = 1, 2, 3, are the Pauli matrices. At low energies, these two operators are matched to the dimension-five photonic dipole operator 5 which is directly responsible for the ℓ i → ℓ j γ radiative LFV decays. Defining the contribution of O eγ to the low-energy effective Lagrangian as where L eγ is its Wilson coefficient, the resulting branching ratios are with m i and Ŵ i the mass and total decay width of the charged lepton ℓ i , respectively. Any theory that induces Q eϕ will also generate O eγ , since these operators transform in the same way under flavor and chiral symmetries [33] and mix under renormalization group evolution [113]. Therefore, one cannot simply get rid of the latter. However, different NP scenarios predict a different balance between the C eϕ and L eγ Wilson coefficients, and this is what determines the magnitude of the allowed HLFV effects in a specific model. Let us consider an example to illustrate this connection: a model inducing predominantly (Q eB ) 12 at the highenergy scale . In this case the operator Q eϕ 12 gets induced due to renormalization group running while O eγ 12 is obtained after matching at the electroweak scale. Since the HLFV operator is induced by operator mixing effects, the resulting coefficient is suppressed and one expects the relation BR(h → µe) ≃ 10 −14 log 2 (m h / ) BR(µ → eγ ), which clearly precludes the observation of the HLFV decay. More generally, in models with L eγ ∼ C eϕ or L eγ > C eϕ , as in the example we just considered, the strong constraints derived from the non-observation of ℓ i → ℓ j γ would imply tiny HLFV rates. In contrast, models predicting L eγ ≪ C eϕ may accommodate sizable HLFV effects. As we will see in section 4.1, the 2HDM is one of such models.

HLFV IN THE 2HDM
We now concentrate on the 2HDM. First, in section 4.1 we particularize the previous model-independent discussion to the case of a 2HDM in order to motivate this model as the perfect scenario to obtain large HLFV rates. Section 4.2 will introduce our notation and conventions for the 2HDM. Finally, we will concentrate on τ µ flavor violation, discuss τ → µγ in the 2HDM in section 4.3 and show some selected phenomenological results on h → τ µ in the 2HDM in section 4.4.

EFT Motivation
A NP model with a second scalar SU(2) L doublet would induce the Q eϕ operator, as shown in Figure 1. This is one of the possible topologies contributing to the HLFV operator, as discussed in great detail in Herrero-Garcia et al. [72]. The scalar must be an SU(2) L doublet for the diagram to be gauge invariant, and would be identified with a second Higgs doublet in a 2HDM. Moreover, it should be noticed that this topology requires both scalar doublets to couple to leptons. The coupling is clearly shown, and ϕ was already assumed to couple to leptons, see Equation (1). For this reason, the 2HDM behind the generation of this topology would be a type-III 2HDM, in which both Higgs doublets are allowed to couple to leptons in a general way 6 . We note that Q eϕ is generated at tree-level in this scenario, thus enhancing the Wilson coefficient C eϕ . We expect where λ and y are the quartic scalar and Yukawa couplings involved in the topology shown in Figure 1 and m the scalar mass term. Next, we consider the generation of O eγ . This operator can be generated by attaching one of the external ϕ lines in Figure 1 to a charged lepton line and then adding the photon. There are, however, three effects that suppress the generation of this operator in the 2HDM: • The O eγ operator gets induced at the 1-loop level.
• Closing the loop by attaching the ϕ line to the charged lepton line introduces a charged lepton mass insertion. • The O eγ operator requires a chirality flip, which in a SM extension with only new scalar fields (such as the 2HDM) implies an additional charged lepton mass insertion.
These considerations allow us to estimate This estimate is known to be very poor due to the presence of several additional contributions in a complete 2HDM, as we will show in section 4.4. Nevertheless, it serves as a good motivation to consider this scenario as potentially promising in what concerns HLFV, since the required hierarchy between C eϕ and L eγ is naturally obtained.

2HDM: Model Basics
In the following, we consider the general 2HDM, usually referred to as type-III 2HDM, and denote the two Higgs doublets as ϕ 1 and ϕ 2 . In contrast to other variants of the 2HDM, no distinction between the two Higgs doublets is introduced. This has two consequences for our discussion. First, both ϕ 1 and ϕ 2 are allowed to couple to all fermion species, and in particular to leptons, a fundamental ingredient for the generation of HLFV effects, see section 4.1. And second, one can perform arbitrary U(2) basis transformations in {ϕ 1 , ϕ 2 } space, without any impact on physical observables. This freedom can be used to go to a specific basis in which only one Higgs doublet acquires a VEV, the so-called Higgs basis [117][118][119]. In this basis, the scalar potential of the model is given by 7 and the Higgs doublets can be decomposed as 7 We follow the conventions of Davidson and Haber [120], with a slightly different notation.
Here φ 0 1 , φ 0 2 and A are neutral scalars, H + is a charged scalar and G + and G 0 are Goldstone bosons. Assuming CP conservation, the CP-even states φ 0 1 and φ 0 2 do not mix with the CP-odd state A, which is a physical mass eigenstate. φ 0 1 and φ 0 2 are related to the mass eigenstates h and H (with m h < m H ) as where s β−α ≡ sin(β − α), c β−α ≡ cos(β − α) and β − α is a physical mixing angle. The lightest CP-even state, h, is identified with the Higgs boson discovered at the LHC. With these definitions at hand, one can derive several relations between the potential parameters and the physical Higgs masses [120], Let us now discuss the Yukawa interactions of the model. Again, we use the freedom to choose specific weak bases. In the Higgs basis for the scalar doublets and the mass basis for the fermions, the Yukawa Lagrangian can be written as +q ρ u u ϕ 2 +q ρ d d ϕ 2 +l ρ e e ϕ 2 + h.c. . (24) Here we denote ϕ a = i σ 2 ϕ * a , with a = 1, 2, and define the diagonal matrices M u = diag(m u , m c , m t ) and M d = diag(m d , m s , m b ). K is the CKM matrix and ρ u,d,e are general 3×3 complex matrices in flavor space, which in the following will be assumed to be Hermitian for simplicity. Using Equations (17) and (18), and expanding in SU(2) L indices, the leptonic part of Equation (24) can be rewritten as L ρ e e R A +ν L U † ρ e e R H + + h.c. , (25) where U is the PMNS matrix. This expression allows us to extract the couplings of the neutral scalars of the model to leptons. By following the same steps with the quarks, one finds the Frontiers in Physics | www.frontiersin.org general expressions where f = u, d, e. He have introduced s f = +1 for down-type quarks and charged leptons and s f = −1 for up-type quarks. Equation (26) must be compared to the general expression in Equation (8). Again, the first term is diagonal, whereas the second may contain off-diagonal entries. We therefore conclude that the ρ e matrix is the source of the HLFV processes that we are about to discuss. Finally, the neutral scalar couplings to a pair of gauge bosons are fully dictated by the gauge symmetry. One has and the couplings to a pair of Z-bosons follow the same proporcionalities.

τ → µγ in the 2HDM
Given the strong experimental bounds on µe flavor violating processes, we will concentrate on τ LFV. In particular, we will consider τ µ LFV, and therefore discuss h → τ µ and the related τ → µγ . The h → τ µ HLFV decay is the main focus of this manuscript and we show some phemenological results in section 4.4. However, in order to assess the observability of this process, one must take into account the strong constraint coming from BR(τ → µγ ), which we now proceed to evaluate in the 2HDM. The τ → µγ radiative decay is induced by the dipole operator defined by Equations (11) and (12). It proves convenient to define the form factors A L and A R as Since we assume the matrix ρ to be Hermitian, |g hτ µ | = |g hµτ |, which implies |A L | = |A R | ≡ |A|. We just need to determine the most relevant contributions to the form factor A in the 2HDM. It is well known that in the 2HDM 2-loop contributions to ℓ i → ℓ j γ may easily dominate over 1-loop ones [121]. The reason is easy to understand. Dipole transitions require a chirality flip. In a 1-loop diagram with a virtual scalar in the loop, two chirality flips take place in the Yukawa vertices, and therefore one more is required in the fermion propagator, giving a total of three. This largely suppresses the loop amplitude, which explains why 2-loop diagrams with only one chirality flip can be dominant even if one pays the extra loop suppression factor of 1/(16π 2 ). In particular, 2-loop Barr-Zee diagrams [122] can easily dominate if the involved scalar fields have large couplings to the virtual fermions or bosons running in the loops. Taking all these ingredients into account, the authors of Davidson and Grenier [19] identified three main contributions to τ → µγ in the type-III 2HDM: • 1-loop diagrams with neutral Higgs bosons and charged leptons in the loop • 2-loop Barr-Zee diagrams with an internal photon and a third generation quark • 2-loop Barr-Zee diagrams with an internal photon and a W-boson These contributions were computed in Chang et al. [123] and are shown in Figure 2. One can write where the different contributions were computed by Davidson and Grenier [19] and are given explicitly in Appendix. Armed with these expressions we are ready to explore the HLFV phenomenology of the type-III 2HDM.

HLFV Phenomenology
Following Aristizabal Sierra and Vicente [9], we show now some phenomenological results on HLFV in the type-III 2HDM. We refer to Crivellin et al. [ [89], and Sher and Thrasher [124] for additional HLFV phenomenological studies in the 2HDM. First, we must make an observation about the type-III 2HDM. As already explained, in this version of the 2HDM one can apply rotations in Higgs space that modify the Higgs VEVs. For this reason, the usual ratio of VEVs tan β is not uniquely defined. Given that we are mainly interested in tau flavor violation, we define [120] tan We note that tan β τ is the physical ratio between the tau Yukawa coupling and √ 2m τ /v, which would correspond to the usual tan β in the Type-II 2HDM.
Let us now discuss our parameter choices. The results presented here are based on a random scan of the parameter space, taking the parameter ranges, 200 GeV < m H < 1000 GeV , 400 GeV < m A < 1000 GeV , These are based on the following considerations and experimental constraints: • It proves convenient to use as input the scalar masses, rather than the scalar potential parameters. These should nevertheless be computed to make sure that they never exceed the perturbative limit of 4π. The 1 , 2 , 4 , 5 and 6 parameters can be computed by means of Equations (19)- (23). The remaining parameters do not have any impact on the observables studied here, but they can be chosen by demanding that the scalar potential is bounded from below [125][126][127].
• A small m H ± − m A mass difference is required by electroweak precision data. In particular, larger values could potentially lead to a T oblique parameter value outside the 1 σ range determined by Baak et al. [128], T ∈ [−0.03, 0.19]. • The lower limit on sin(β − α) is motivated by the fact that the observed Higgs boson has SM-like couplings to fermions and gauge bosons. Lower values of sin(β −α) could potentially induce deviations, see Equations (26) and (29), in tension with the experimental results. In our analysis we consider the CMS measurements presented in CMS Collaboration [129] and require that the signal strengths for h → ττ , bb, WW * , ZZ * , defined as µ = (σ × BR) / (σ × BR) SM , are within the CMS 1 σ ranges [129]. For the determination of the Higgs production cross-section we assume gluon fusion. • The lower limit on m H ± is motivated by flavor physics (mainly B physics, see for instance [130]).
Finally, our scan also fixes m h = 125.1 GeV [131]. In what concerns the Yukawa matrices, and in order to reduce the number of free parameters, we will consider specific textures for the ρ matrices. Inspired by the Cheng-Sher ansatz [132], we express ρ e as By construction, κ τ τ = 1. However, the other entries of the κ matrix are free parameters. In particular, κ τ µ = κ * µτ is the relevant parameter giving rise to the h → τ µ and τ → µγ decays. In our random scan we will take 0.1 < |κ τ µ | < 3.0. For the quark ρ matrices we assume the usual Type-II textures The parameters are fixed as explained in the text. The vertical lines correspond to the current bound BR(τ → µγ ) < 4.4 × 10 −8 [98] and the expected Belle-II sensitivity, estimated to be ∼ 10 −9 [99]. Finally, the horizontal line indicates the limit by the CMS collaboration BR(h → τ µ) < 0.0025 [95].
This ansatz is particularly convenient since it ensures compatibility with the (already constraining) experimental bounds on the Higgs boson couplings to quarks. Furthermore, it can be regarded as a departure beyond the popular Type-II 2HDM, with the only deviation in the τ µ coupling [133] 8 . After these preliminaries, we are ready to show some results on h → τ µ. Figure 3 shows BR(h → τ µ) as a function of BR(τ → µγ ). The vertical lines shown in this figure correspond to the current experimental bound BR(τ → µγ ) < 4.4 × 10 −8 [98] and the expected sensitivity of the Belle-II experiment, of about ∼ 10 −9 [99]. The horizontal line corresponds to the limit BR(h → τ µ) < 0.0025, set by the CMS collaboration [95]. As can be seen from this figure, the correlation between these two observables is not exact, and this can be traced back to the different contributions to τ → µγ , which might even cancel in some cases. We find that, in general, the dominant contribution to the τ → µγ amplitude comes from 2-loop Barr-Zee diagrams with internal W bosons, although the other contributions typically play a relevant role as well.
The main qualitative message that one can extract from Figure 3 is that the type-III 2HDM can induce h → τ µ rates close to the current bound while being in agreement with all experimental constraints. All LFV observables increase with κ τ µ , and in some regions of the parameter space they can be close to their current experimental limits, explicitly shown in Figure 3. These regions are characterized by tan β τ 2, sin(β − α) ∼ 0.9 and κ τ µ 0.1. Finally, some additional comments are in order. A Higgs doublet with µτ LFV couplings can also address the longstanding muon g − 2 anomaly. This was studied in relation to the LFV decay h → τ µ in Aristizabal Sierra and Vicente [9], Davidson and Grenier [19], Omura et al. [34], Liu et al. [47], Benbrik et al. [53], Altmannshofer et al. [70], Wang et al. [76], and Crivellin et al. [130] and very recently in Iguro et al. [134] and Wang and Zhang [135]. It could also be linked to the popular b → s [31,48] or b → c [42] anomalies observed in B-meson decays or be a crucial ingredient for lepton-flavored electroweak baryogenesis [80]. In fact, the type-III 2HDM with generic Yukawa couplings has a very rich flavor phenomenology, see Crivellin et al. [130]. The analogous quark flavor violating decay h → bs was studied in Crivellin et al. [136]. It is also remarkable that the 2HDM with a BGL symmetry [137][138][139] can also lead to large h → τ µ branching ratios, strongly correlated with other flavor observables, as shown in Botella et al. [46]. Here we concentrated on the 2HDM. For HLFV studies in other multi-Higgs doublet models see the interesting works [140][141][142][143].

2HDMs, NEUTRINO MASSES AND HLFV
Neutrino flavor oscillations constitute the only existing experimental proof of LFV. Since these are sourced by non-zero neutrino masses and mixings, it is interesting to discuss their possible connection to HLFV [72,79]. First, one should notice that while the existence of non-zero neutrino masses and mixings implies the violation of lepton flavor, the observation of lepton flavor violating processes does not require neutrinos to be massive. In fact, there are many examples of the latter, the 2HDM being the simplest one. Indeed, in the model presented in section 4.2, the most general 2HDM, neutrino masses are exactly zero but processes such as ℓ i → ℓ j γ or h → ℓ i ℓ j are perfectly possible. Similarly, neutrino masses vanish in the Minimal Supersymmetric Standard Model, but LFV processes are indeed induced if the slepton soft mass contain off-diagonal entries. Other examples are also known, see for instance [144].
There are, however, many neutrino mass models that require the introduction of a second Higgs doublet, and these may potentially lead to an interesting connection between the generation of neutrino masses and HLFV. One of the most popular examples of this link is the Zee model [145], a setup that induces neutrino masses at the 1-loop level 9 . This model can actually be regarded as an extension of the general 2HDM with the addition of a singly-charged scalar, k ∼ (1, 1, 1) . (38) If both lepton doublets couple to leptons, as in the type-III 2HDM, the simultaneous presence of the Yukawa term fl c ℓ k and the trilinear scalar potential term µ ϕ 1 ϕ 2 k † breaks lepton number in two units, thus inducing Majorana neutrino masses. Therefore, the model contains all the ingredients to induce neutrino masses and observable HLFV rates. Interestingly, these two consequences of the Zee model are connected in a non-trivial way. Assuming for simplicity f eµ = 0, neglecting the electron mass compared to the muon and tau masses, and keeping the term proportional to the muon mass in the (3,3) element to get three massive neutrinos, the neutrino mass matrix is given by Herrero-García et al. [83] where A is a dimensionless combination of model parameters, containing the corresponding loop function, s 2δ ∝ µ quantifies the mixing in the charged scalar sector and we denote s β τ = sin β τ . The ρ e matrix was introduced in Equation (24). We see that in order to accommodate the measured leptonic mixing angles (see for instance the global fit [147]) both ρ τ µ e and ρ τ e e must be different from zero. Therefore, the Zee model leads to correlations between the leptonic mixing angles and the h → τ µ and h → τ e rates. These can be used to set a lower limit on the HLFV rates [83]. For instance, one finds for normal (inverted) neutrino mass ordering. We refer to Herrero-García et al. [83], where a detailed exploration of the parameter space of the Zee model is performed, concluding that the model can accommodate large HLFV rates, even exceeding the current bounds. Similar findings were recently found in Nomura et al. [91], where a Zee model supplemented with a flavor-dependent U(1) symmetry was considered.
There are other neutrino mass models including a second Higgs doublet. For instance, in left-right symmetric models [148][149][150][151][152] one usually introduces a scalar field that is a doublet of both SU(2) L and SU(2) R . This bidoublet can be denoted by and decomposed as The scalar representation has two gauge invariant Yukawa couplings to leptons and can be regarded at energies below the SU(2) R breaking scale as a pair of SU(2) L doublets. However, these setups cannot be identified with a type-III 2HDM since the left-right symmetry require the two lepton Yukawa matrices to be Hermitian, thus strongly restricting the allowed parameter space. Furthermore, current limits on quark flavor violation require the second CP-even mass eigenstate, H, to have a large mass, m H 25 TeV [153], suppressing all HLFV effects. Therefore, the minimal left-right models would have to be enlarged with additional scalar fields in order to be able to provide large HLFV rates [72].

SUMMARY AND DISCUSSION
In this mini-review we have discussed Higgs lepton flavor violating decays, such as h → ℓ i ℓ j , in the context of the general 2HDM. After motivating this scenario with some modelindependent considerations, we have explicitly shown that the 2HDM can indeed allow for large HLFV rates while being in perfect agreement with the experimental constraints at low and high energies. A possible connection to the mechanism behind the generation of neutrino masses is also discussed.
The 2HDM must face many stringent constraints in order to generate large HLFV rates. The first tension comes from existing measurements of the Higgs boson couplings to fermions and gauge bosons, which already place bounds on the parameter that controls the mixing between the two CP-even scalars in the model, sin(β − α). This must necessarily deviate from 1 in order to allow for non-standard Higgs decays, but not too much in order to be in agreement with the constraints on the Higgs boson couplings. The dipole transitions ℓ i → ℓ j γ also set very strong limits on the LFV parameters, and these cannot be avoided by any flavor symmetry. However, the associated operators turn out to be suppressed in the 2HDM compared to the operators leading to HLFV, thus increasing the chances to have observable effects of the latter. Taking all these constraints into account, as well as those from electroweak precision data or direct LHC searches, we find that the 2HDM can accommodate HLFV rates arbitrarily close to the current limits, hence making them a very attractive way to search for NP.
The determination of the properties of the recently discovered Higgs boson is among the current priorities for the particle physics community. The increasingly precise measurements of the Higgs couplings and decay rates might eventually reveal a deviation from the SM expectations and hint toward the existence of NP. In particular, the search for HLFV has just started. A positive signal could shed light on questions as central as the flavor puzzle or the fundamental nature of electroweak symmetry breaking.

AUTHOR CONTRIBUTIONS
The author declares that he contributed to the design and implementation of the research, to the analysis of the results and to the writing of the manuscript.

ACKNOWLEDGMENTS
I am very grateful to my collaborators in the subjects discussed in this review. In particular, I thank Diego Aristizábal Sierra for many long and enjoyable discussions on Higgs lepton flavor violating decays. I also thank Luca Fiorini for drawing my attention to updated results on h → ℓ i ℓ j by the ATLAS collaboration. I acknowledge financial support from the Spanish grants SEV-2014-0398 and FPA2017-85216-P (AEI/FEDER, UE) and SEJI/2018/033 (Generalitat Valenciana) and the Spanish Red Consolider MultiDark FPA2017-90566-REDC.