Predictions for Muon Electric and Magnetic Dipole Moments from $h \rightarrow \mu^+ \mu^-$ in Two-Higgs-Doublet Models with New Leptons

We calculate chirally enhanced corrections to the muon's electric and magnetic dipole moments in two-Higgs-doublet models extended by vector-like leptons, and we explore a sharp correlation between $h \rightarrow \mu^+ \mu^-$ and the muon's dipole moments in these models. Among many detailed predictions, for a model with new leptons with the same quantum numbers as standard model leptons, we find that $0.38 \lesssim \tan \beta \lesssim 21$ necessarily requires a muon electric dipole moment to be observed at near-future experiments, assuming $h \rightarrow \mu^+ \mu^-$ is measured within $1\%$ of the standard model prediction for the current central value of the measured muon magnetic moment. In all studied models, the predicted values of the electric dipole moment can reach up to current experimental limits. Moreover, we show that in some models there can be two sources of chiral enhancement, parametrizing the correlation between $h \rightarrow \mu^+ \mu^-$ and the dipole moments by a complex number. This leads to sign-preferred predictions for the electric dipole moment.


I. INTRODUCTION
Recently, measurements of fundamental properties of the muon have gained renewed interest in serving as harbingers of new physics. Measurements of the standard model (SM) Higgs boson decay into muons allow for a deviation by more than a factor of 2 compared to the SM prediction [1], while the most recent measurement [2] of the muon's anomalous magnetic moment, ∆a µ , suggests that the experimental world average is discrepant by more than five standard deviations from the SM prediction of 2020 [3]. Current theoretical and experimental efforts aim to understand this result [4,5]. Additionally, the muon electric dipole moment has been a popular tool in studying the possible effects of large CP violation from new physics [6].
In this paper, we calculate contributions to the electric and magnetic dipole moments of the muon in a two-Higgs-doublet model (2HDM), assuming a Z 2 symmetry which enforces type-II couplings to SM leptons, extended with vector-like leptons. We explore representations of new leptons up to and including SU (2) triplets that were previously studied in SM extensions [7,8]. Such models (and other possibilities, see for example Refs. [9][10][11][12]) can generate simultaneous corrections to the muon mass and dipole moments as δm µ ≃ λ 3 v 3 /M 2 , ∆a µ ≃ m µ vRe[λ 3 ]/16π 2 M 2 , and d µ ≃ evIm[λ 3 ]/32π 2 M 2 , where λ 3 and M 2 are the products of the couplings and mass scales of new particles. The corrections to dipole moments enjoy a chiral enhancement of λv/m µ compared to the naive estimate, ∆a µ ≃ m 2 µ Re[λ 2 ]/16π 2 M 2 and d µ ≃ em µ Im[λ 2 ]/32π 2 M 2 , and allow for the largest scales of new physics to explain the anomaly. Chirally enhanced corrections to the muon dipole moments from vector-like leptons may also be connected to new physics relevant for models of dark matter [13][14][15][16][17][18] and the Cabibbo anomaly [19,20]. 1 In a 2HDM with type-II couplings (2HDM-II), these contributions can be even further enhanced by factors of tan 2 β [25,26] [where tan β is the ratio of the vacuum expectation values (VEV) of the two-Higgs doublets]. For the models not explored in [25,26], we find yet another source of chiral enhancement. For all of the models that we consider, we present 1 Additional explanations for ∆a µ consistent with searches for new Higgses can be accomplished exclusively in 2HDMs, requiring light scalar and pseudoscalar masses below the electroweak scale and large couplings to leptons, such as in type-X and flavor-aligned 2HDMs [21,22], and even in the general 2HDM where flavor-changing currents are permitted [23]. A muon-specific 2HDM can also accomplish this [24]. complete expressions of both dipole moments in the mass eigenstate basis. These formulas are completely generic and can be applied to any model with new leptons and an extended Higgs sector. These could be applied, for instance, to other types of 2HDMs [27][28][29][30][31][32][33][34] or more generic extensions of the SM with vector-like leptons and new scalar particles [33,[35][36][37][38][39].
In addition to calculations in the mass eigenstate basis, we also present results using a SM effective field theory (SMEFT)-like 2HDM effective field theory working at dimension six, where the SU (2) L × U (1) Y symmetry is linearly realized and both Higgs doublets are kept light while vector-like leptons are integrated out. These results reproduce the leadingorder v 2 /M 2 corrections in the mass eigenstate basis. We find that working in the so-called Higgs basis [40][41][42][43] greatly simplifies the matching calculation, leading to strikingly simple parametrizations for the muon dipole moments.
Although chirally enhanced dipole moments can allow for extremely heavy scales of new physics currently out of reach at the Large Hadron Collider (LHC) and possible future colliders, it was recently shown that these corrections are highly correlated to the h → µ + µ − decay rate [9,[44][45][46][47][48][49][50], although this was explored earlier for lighter new leptons [7,8].
In particular, the dimension-six operator generating the muon dipole operators is directly proportional to the one that modifies the muon mass by a real factor, k. This allows for a simultaneous parametrization of three observables−∆a µ , d µ , and h → µ + µ − −through a single model-dependent factor, providing a novel way to test high scales in the theory that are directly inaccessible to current and future colliders. In a 2HDM-II extended with vector-like leptons, the proportionality factor can now be complex due to the interference of multiple sources of chiral enhancement. This was first pointed out in a model-independent fashion in [47,48]. The complex factor leads to asymmetric and sign-preferred predictions of d µ , previously not investigated in the literature. We find a similar behavior can also occur when including subleading loop corrections to the muon mass. The correlation of muon observables can be used to obtain novel bounds on the masses and parameter space of 2HDMs with vector-like leptons inferred from current and future limits of h → µ + µ − and projected measurements of d µ .
The paper is organized as follows. In Sec. II we introduce models of vector-like leptons in the context of 2HDM-II allowing for complex Yukawa couplings and discuss our conventions for the scalar potential. In Sec. III we present general formulas for contributions to ∆a µ and d µ in the mass eigenstate basis. In Sec. IV we present calculations of the dipole moments in the 2HDM-II effective field theory. In Sec. V we explore the connection between ∆a µ , d µ , and h → µ + µ − , present detailed results, and extend discussions from previous works.
We conclude in Sec. VI. We also include lengthy appendices providing details of the 2HDM scalar potential, representations for new leptons and our conventions for calculations in the mass eigenstate basis, approximate formulas for couplings of leptons, and further details of the effective field theory calculations.

A. Model setup and parameters
We consider a 2HDM extended with charged vector-like lepton doublets L L,R and singlets E L,R , as well as neutral vector-like singlets N L,R . We assume a Z 2 symmetry, enforcing couplings of the SM leptons to the Higgs doublets as in a type-II 2HDM. This model was first considered in [51] and studied extensively in [25,26] in order to explain the anomalous magnetic moment of the muon. The most general Lagrangian assuming only mixing of second-generation leptons with new leptons is described by Quantum numbers of the fields are given in Table I. Couplings of H d and H u to other SM fermions are like those in the typical type-II 2HDM. Additional models with different representations of vectorlike leptons are described in Appendix C. An alternate version of these models exists motivated by the minimal supersymmetric SM extended with vectorlike leptons [25]. Because the superpotential is holomorphic, terms involving H † d and H † u are forbidden (but similar terms appear through H u and H d , respectively). In principle, each of the eight couplings and three vector-like masses can be complex. Field redefinitions can be chosen such that several phases of the parameters are unphysical. The Lagrangian in Eq. (1) admits four additional physical phases (two in the charged and two in the neutral lepton sectors) that cannot be rotated away. For example, we could take the mass parameters M L,E,N and all Yukawa couplings to be real except for λ, λ, κ, and κ (although we do not make this assumption). We explore the impact of these complex couplings later in Sec. III.  The doublet components are defined as After electroweak symmetry breaking (EWSB), the neutral components of the Higgs doublets GeV and their ratio is parametrized by v u /v d = tan β. The diagonalization procedure for rotating these fields to the physical basis is outlined in Appendix C and also provided for other representations.
Integrating out the heavy lepton fields reduces the above Lagrangian to where The parameter m LE µ would be the muon mass if y µ were zero. The effective field theory Lagrangian additionally contains dimension-six operators such as C µHulL µ R H d H † u H u ; however, after EWSB, they vanish at tree level and do not contribute to the muon mass. In Sec. IV we explore non-zero contributions of these types of operators to the muon mass at loop level via quartic interactions in the scalar sector.
After EWSB, the physical muon mass m µ is given by and the muon Yukawa coupling is Since (λ h µµ ) SM = m µ /v, the decay rate of h → µ + µ − compared to its SM rate is which is currently limited by R h→µ + µ − ≤ 2.2 [1].

B. Softly broken 2HDM potential
The most general scalar potential of the two-Higgs doublets, H d and H u , consistent with our conventions 2 while allowing soft Z 2 -breaking terms is given by where the explicit "·" contracts SU (2) doublets through the antisymmetric ϵ ij , e.g., The exact Z 2 -symmetric potential is modified by the additional free parameter m 2 12 which is required to softly break the symmetry. For further details, see Appendix A.
After diagonalization of the scalar fields to the mass eigenstate basis, the quartic couplings 2 In 2HDM models with a Z 2 symmetry, the doublets are usually defined as [52]. Matching to our definitions of H d,u requires rotating the λ 1 , λ 2 , λ 3 , λ 4 , and λ 5 can be written as sin β cos β , sin β cos β , If we consider m 2 h , this regime also enforces that α → β − π/2 (i.e. the alignment limit), in which the light eigenstate h behaves as the SM Higgs boson [52]. In this limit, the above couplings highly simplify to Notice that these couplings are not independent from one another in the decoupling limit, but rather λ 2 = λ 1 / tan 4 β and λ 3 = λ 1 / tan 2 β. In addition, the stability conditions of Eq. (A1) all reduce to the same lower bound: Now, this range is further restricted by imposing perturbative unitarity constraints on the above couplings. The condition for unitarity is expressed in the form of partial-wave scattering from W ± and Z bosons in the high-energy limit, which are realized by the Goldstone boson equivalence theorem through their longitudinal (scalar) modes in 2 → 2 scattering [53][54][55]. Applying this analysis to the softly broken 2HDM-II (see [56][57][58] for details while using the condition for unitarity in [59]), partial-wave unitarity and perturbativity require Applying these unitarity constraints on each coupling, we find that the heavy Higgs mass is bounded above: sin β cos β , sin β cos β , sin β cos β .
For an exact Z 2 symmetry, these unitarity conditions lead to upper limits on each individual Higgs mass discussed in [60,61], similar to the Lee-Quigg-Thacker upper mass limit obtained in a single-Higgs-doublet model [53,54]. However, with a softly broken Z 2 symmetry, m 2 H ≃ m 2 A ≃ m 2 H ± can be arbitrarily large when appropriately choosing m 2 12 for a given tan β. In the limit m 2 12 ≫ 8πv 2 , the bounds in Eq. (12) all become equal. For sufficiently heavy m 2 H ≃ m 2 A ≃ m 2 H ± , the entire range of tan β is allowed by direct searches [62,63].

III. CONTRIBUTIONS TO THE MUON DIPOLE MOMENTS
We turn our attention to calculating the anomalous magnetic (∆a µ ) and electric (d µ ) dipole moments of the muon. They can be extracted via the effective Lagrangian (where we use e > 0 and the time-like metric) which are both calculated via Fig. 1. The contributions to (g − 2) µ discussed in [25] are not specific to any type of 2HDM and, in a similar manner, the contributions to d µ presented here are general for any model with an extended Higgs sector. For completeness, we list contributions to ∆a µ in addition to d µ . Couplings to mass eigenstates are defined identically as therein, but for complex Lagrangian parameters studied here and are collected in Appendix B. Approximate formulas for couplings relevant to (g − 2) µ and d µ are given in Appendix C.
In the mass eigenstate basis, charged or neutral lepton eigenstates f a couple to the Z boson through   and generate contributions to ∆a µ and d µ as where x a Z = m 2 ea /M 2 Z parametrizes the following loop functions F Z (x) and G Z (x): In a similar way, couplings of charged and neutral leptons to the W ± boson are given by The corresponding contributions to ∆a µ and d µ are where For neutral Higgs scalars ϕ = h, H, A, their couplings to charged leptons are defined as where their contributions to ∆a µ and d µ are given by with x a ϕ = m 2 ea /m 2 ϕ and The couplings that describe interactions between charged and neutral leptons to the charged The contribution from the charged Higgs to ∆a µ and d µ are Models involving doubly charged fermions coupling to the W ± boson and a singly charged fermion are described by Contributions from the W ± boson to ∆a µ and d µ are given by as defined above. Finally, the couplings that describe interactions between singly and doubly charged fermions to the charged Higgs boson H ± are The charged Higgs contributes to ∆a µ and d µ by /m 2 H ± and with the functions defined above, The mass operators with C where µHu ) tan 2 β, including loop corrections. Note that there are additional operators involving H 2 ; however, they do not contribute to the mass or dipole moments. After EWSB, the Wilson coefficients C µB and C µW combine to generate the muon's dipole moments in Eq. (13) and we find where C µγ = cos θ W C µB − sin θ W C µW and θ W is the weak mixing angle. We now present the contribution to ∆a µ and d µ for the five representations of vector-like leptons listed in [7], calculated in the unbroken SU (2) L × U (1) Y theory from diagrams in Fig. 2. Allowing for only down-type couplings in the limit M L,E ≫ v d × (λ L , λ E , λ, λ), ∆a µ and d µ become Notice how each representation experiences a tan 2 β enhancement in addition to the standard chiral enhancement. The factors L,E ) represent contributions from light (SM) and new scalar fields, respectively, parametrized by x   Figure 3: Caption M → 1, with λ * = 0 (second column) and λ * =λ (third column). The Z 2 charge assignments for L and E are as they are in Table I. and are collected in Table II.
Contributions involving theλ coupling arise from diagrams where each heavy fermion propagator contains its mass insertion. However, diagrams where each propagator contains its momentum term are responsible for generating contributions from λ. In the latter situation, diagrams involving λ in the top row of Fig When we consider the other limit, m h ≲ M 2 ≪ M , all terms with λ vanish and Q 2 → Q 1 ≡ Q. For all representations given in Table II, Q reduces to 1, 5, or 9. Note that the heavy Higgses are not required to be near the EW scale; this is already a good approximation Scenarios involving the mixing between new heavy leptons and the muon via only the SM Higgs were first studied in the mass eigenstate basis in [7,8,65] and revisited in [48,64] in the SMEFT landscape. To provide insight into these scenarios, although unphysical, one given in Appendix C, expanding to ∼ O(v 2 /M 2 L,E ).   can take the limit tan β → 0 in Eq. (47) which isolates loop contributions only from the SM doublet H 1 . We find agreement with results in [64] in every representation. However, we disagree with [7], notably, in representations including doubly charged vector-like fields.

6-Dimensional Contributions to C µH
The quartic couplings λ 1 , λ 2 , λ 3 , λ 4 , and λ 5 in the scalar potential are relevant contributions at the one-loop level to the Wilson coefficient C µH 1 . In principle, the combination of tan β, M 2 , and large couplings can generate sufficiently large contributions to overcome loop-suppressed effects, competing with the tree-level contribution C (tree) Turning our attention to the main representation of the paper, the contribution to C µH 1 assuming only down-type couplings and the same common mass M L,E = M for new leptons up to one loop in Fig. 3 is where x (1,2) M parametrizes the loop functions If we assume once more that x M → ∞ while using Λ 3 = 2λ 1 / tan 2 β in Appendix A, the muon mass is additionally modified by the above Wilson coefficient: In principle, for the representations listed in [7], one can construct similar diagrams to Fig. 3 and calculate its full one-loop correction to the muon mass from the same dimensionsix operators. However, we expect their contributions to be of similar order in all cases and it suffices to present results in the 2 −1/2 ⊕ 1 −1 representation. We will discuss the impact of these loop effects in the following section.

V. RESULTS
The current experimental world average of the muon's anomalous magnetic moment deviates from the SM prediction by 5.1 σ [2, 3], and the current upper bound on the muon's electric dipole moment is [66] |d µ | ≤ 1.9 × 10 −20 e · cm.
The Fermilab Muon g − 2 Collaboration estimates an improvement in d µ up to a level of [67] |d µ | < 10 −21 e · cm, whereas the Paul Scherrer Institute (PSI) intends to host a new experiment involving the frozen-spin technique [68], and projects a measurement of Precision electroweak measurements of Z and W ± bosons constrain possible modifications of couplings to the muon at the 0.1% level which, in the limit of small mixing, translate to the following bounds on λ L and λ E [25] specifically in the 2 −1/2 ⊕ 1 −1 model: respectively, at 95% C.L. Additionally, we impose a bound on the mass parameters M L and M E , such that M L > 800 GeV and M E > 200 GeV to generically satisfy constraints placed by searches for new leptons [69][70][71].

A. 2HDM-II ellipse of dipole moments
It was first pointed out in [47,48] that the connection between the three observables h → µ + µ − , ∆a µ , and d µ is a byproduct of a correlation between the Wilson coefficients that modify the muon mass (C µH 1 ) and generate the muon dipole moments (C µγ ) through a real model-dependent factor k. In the 2HDM-II, this relation is This connection parametrizes the new physics contributions to the dipole moments and modification of h → µ + µ − in a given model on an ellipse defined by This is derived from rewriting the definition of R h→µ + µ − in Eq. (7) in terms of the moments in Eq. (46). The k factor is defined for a class of models in which C µH 1 is generated at either tree or loop level [47,48]. Although k may include radiative corrections, we will not consider their effects until Sec. V C.
From Eq. (47), the k factor for the 2HDM-II extended with vector-like leptons is which is generalized for different representations as summarized in Table II. As follows from Eq. (59) as well as mentioned in [47], only models with k ≲ 19 or 483 ≲ k ≲ 750 are consistent with ∆a µ within 1σ and |d µ | = 0 assuming R h→µ + µ − = 1 ± 10%. In other words, models with 1.31 ≲ (Q 1 + Q 2 tan 2 β) ≲ 33 or (Q 1 + Q 2 tan 2 β) ≲ 0.84 necessarily require a nonzero value for |d µ | with the expected precision level of R h→µ + µ − at the LHC.
Further improvement of the precision of R h→µ + µ − to 1% would extend these ranges to 1.28 ≲ When heavy Higgses are significantly lighter than vector-like leptons, we saw in the previous section that Q 1 = Q 2 ≡ Q = 1, 5, or 9 for all five representations. In this limit, the k factors become k = 64π 2 /Q(1+tan 2 β). In Fig. 4, we see that as Q increases for fixed tan β, the center of the ellipse shifts to larger vales of ∆a µ and d µ . The experimental value of ∆a µ restricts the parameter space of models to specific ranges of d µ and R h→µ + µ − . For example, for tan β = 1, models with Q = 9, 5, and 1 predict |d µ | between ∼ (6 − 12) × 10 −22 e · cm, This is just above the current experimental limit.
We show the full range of R h→µ + µ − consistent with experiment and predicted values of |d µ | for the central value of ∆a µ as tan β varies in Fig. 6 for models with Q = 1. The solid and dashed lines represent the smaller and larger of the two solutions for tan β [or k through Eq. (60)] that predict the same value of |d µ | and R h→µ + µ − . When 0.41 ≲ tan β ≲ 6.4, R h→µ + µ − = 1 ± 10% necessarily requires values of d µ that can be seen at PSI. As tan β increases, the predicted range of |d µ | extends from zero to current experimental limits.
Further assuming R h→µ + µ − = 1 ± 1% extends the range of tan β that can be fully tested at PSI to 0.38 ≲ tan β ≲ 21. Note that the blue hatched region is where the top Yukawa coupling becomes nonperturbative.
In the left plot of Fig. 7, we plot curves explaining ∆a µ within 1σ for a common scale of new physics versus the size of couplings for different values of tan β = 1, 2, 5, 20, and 50 as well as the SM+VL scenario for Q = 1. As tan β increases, the slope of the contours increases and further limits the allowed range of the predicted mass spectrum. In the right panel, the red and blue curves highlight the behavior of, respectively, R h→µ + µ − = 0.9 and  M )], meaning that solutions where k < 0 for R h→µ + µ − > 1 are never possible (as opposed to some specific FFSand SSF-type scenarios considered in [48]). For example, when d µ → 0 for the dashed red (blue) curves, future LHC precision measurements of λ h µµ at the 10% (1%) level will rule out models with tan β ≲ 6.4 (21), restricting the mass spectrum to M ≲ 18 (10) where the axes of the ellipse in the ∆a µ − d µ plane are shifted by the phase angle ϕ k .
The center of the ellipse is now located at (2m 2 µ cos ϕ k /|k|v 2 , em µ sin ϕ k /|k|v 2 ), which is no longer symmetric along |d µ | = 0. The asymmetry along with a large complex phase suggests that certain models may shift the entire ellipse to a region where a nonzero |d µ | is always predicted. This occurs when the edge of semiminor axis intersects with |d µ | = 0 and is satisfied whenever R h→µ + µ − ≤ sin 2 ϕ k .
In Fig. 8, we demonstrate this behavior for a variety of phases of k = 50e iϕ k for ϕ k = where the λ coupling is present alongside λ, such as the 2 −3/2 ⊕ 1 −1 model, and similar arguments can apply to 2 −3/2 ⊕ 3 −1 , 2 −1/2 ⊕ 3 0 models as well. In the limit that M ≃ M 2 for arbitrary λ, k = 64π 2 /(5 + tan 2 β(17 + λ * /λ)/6). In order to generate |k| = 50 with phases between ϕ k = −π/4 and −π/2 for tan β = 2, the magnitude of |λ * /λ| ≃ 17.4 − 31.0 and its phase should lie between ϕ λ * /λ ≃ −0.88 + π and −0.66 + π rad, respectively. Notice that the π phase shift is required for Re[λ * /λ] to cancel with the other real terms in order to generate a large magnitude of k, which prefers a larger tan β for a more comparable size of couplings (or vice versa as presented here). Im[λ * /λ] is responsible for generating the phase of k. We note that a wide variety of scenarios can accommodate these large ratios of couplings while respecting perturbativity and EW precision constraints needed to explain ∆a µ . However, for ∆a µ , because of the partial cancellation from the real part of λ * /λ, a large complex phase in the combination λ L λ E λ is required. 5

C. One-loop corrections to k
In the Higgs basis, additional contributions from quartic couplings can appear from oneloop corrections to the mass operator C µH 1 [Eq. (48)]. The full one-loop corrected k factor for our main model of interest, new leptons with quantum numbers identical to SM leptons, Notice that C µH 1 now also involves λ at loop level, which in general can generate a nonzero phase for complex Lagrangian parameters. while the contribution is nearly ∼ 16% everywhere for λ 1 = 1. We note that the behavior is similar in both plots because, irrespective of sign, there is a partial cancellation between λλ and (λλ) * terms, leaving terms such as |λ| 2 , |λ| 2 unaffected by their relative sign. In fact,    Table II, the same pattern of contributions occurs for both when λ * = λ and produce integer values.
The additional corrections modify the ellipse equation through an additional phase in the k factor. In fact, for k ∈ C, a striking feature occurs: situations exist where the prediction for d µ is asymmetric and sign-preferred, shifting the ellipse to a region where a portion of the parameter space can be ruled out while being consistent with ∆a µ . This behavior appears for models whose contribution to C µγ also comes from an additional source for generally complex λ * , whenever ϕ λ * ̸ = ϕ λ + πn for n = 0, 1, · · · and M 2 ≃ M . Other sources of complex k can occur via subleading corrections to the muon mass, generating up to 40% of the contribution compared to the tree-level piece when all couplings are near their perturbativity limits in nearly the entire range of tan β. However, this will reduce the magnitude of k by the same amount, causing the center of the ellipse to shift to larger ∆a µ while predicting a large nonzero d µ within the range of R h→µ + µ − = 1 ± 10% [47], increasing the discovery potential of these effects in near-future experiments.
The mass and couplings of the muon remain among the last vestiges of indirect hints of new physics beyond the SM. The need to precisely measure these quantities is currently driving new experimental efforts [67,68] and inspiring the next generation of particle colliders [72][73][74][75][76][77]. The complementarity of future measurements of the muon dipole moments and h → µ + µ − will provide an important road map to understanding high scales of new physics whose presence may already be leaving clues in experiments. λ 1 , λ 2 > 0, λ 1 λ 2 + λ 3 > 0, and After EWSB when the neutral components of the doublets acquire a VEV, ⟨H 0 d ⟩ = v d and ⟨H 0 u ⟩ = v u , the rotation angle α diagonalizes the CP -even scalar fields h and H to the physical basis and yields From the original basis, we can rotate via the angle β to the Higgs basis where the SM and additional scalar fields become separated: We see that in the alignment limit when β − α = π/2, the two doublets simplify to where H 1 contains only the SM degrees of freedom and H 2 contains the additional Higgs fields in this new basis. From Eqs. (A7) and (A8), one notices that ⟨H 0 1 ⟩ = v and ⟨H 0 2 ⟩ = 0, such that H 1 becomes the SM Higgs doublet in the alignment limit. Inverting Eq. (A6), the where we used the fact that iσ 2 = ϵ and ϵ 12 = +1. Inserting these definitions into the 2HDM-II scalar potential (8), we find in this basis If we utilize the extrema conditions of the potential in Eq.

Appendix B: Models and mass eigenstate basis
First, we introduce models and define their mass eigenstate basis in each of the five representations. We additionally present couplings of lepton eigenstates to SM bosons Z, W ± , h and new bosons H, A, H ± . These couplings are defined with a common notation and will be further characterized by diagonalization and Yukawa matrices in each individual representation.
After EWSB, the Lagrangian of Eq. (1) gives the charged lepton mass matrix: and the neutral lepton mass matrix, We construct the mass eigenstate basis in the following steps. First, we find U L,R such that In this procedure, U † L M U R becomes diagonal where the diagonal entries are the physical masses up to possible phases. Applying this procedure to the mass matrices M e and M ν above, we have where m µ , m e 4 , m e 5 , m ν 4 , and m ν 5 are the physical masses of the leptons. Further, we define so that and The diagonalization matrices U L and U R represent rotations to the mass eigenstate basis e L → U e Lê L , e R → U e Rê R and ν L → U ν Lν L , ν R → U ν Rν R . This procedure will be also followed for other representations of leptons.
In Appendix C we present approximate formulas for the diagonalization matrices and resulting couplings. Those formulas are valid only when following the above procedure.
In this representation, the doublet L and triplet E are defined as and couple to H d and SM leptons via the following Lagrangian: After EWSB, the mass matrices are Additionally, there is a doubly charged mass eigenstate E −− with mass M E . Hence, The doublet with singly and doubly charged components is defined as The Lagrangian for this representation becomes We remind the reader that the explicit · corresponds to contracting SU (2) indices by ϵ ij , such . Alternatively, we could continue to contract SU (2) indices with δ ij by introducing ϵH † d ≡ H d . The difference between the two notations is an overall minus sign for the term λ H d L L E R .
After EWSB, the the mass matrices are which is diagonalized by U e L and U e R as in Eq. (B7), but with the replacement of λ → −λ. There is one massless neutral eigenstate ν µ and one doubly charged eigenstate L −− with mass M L . Hence, The doublet L containing both a singly and doubly charged field as well as the triplet E are defined as which couple to the SM leptons in the following way: After EWSB, the mass matrices are (B21) Finally, this representation is described by the following doublet L and triplet E: These fields interact with the SM leptons in the following way: The mass matrices after EWSB are and

Couplings of leptons to Z and W ± bosons
The couplings of leptons to the Z boson come from kinetic terms in the Lagrangian, whereby the covariant derivative is For couplings of left-and right-handed singly charged leptons to the Z boson, we have whereby and The couplings to W ± boson are defined by and Additionally, for representations with doubly charged fermions, the left-and right-handed couplings are The representation-dependent ξ i and χ i factors for each set of left-and right-handed Z and W ± gauge couplings are collected in Table III.

Couplings to Higgs bosons
The softly broken 2HDM-II scalar potential is given by Eq. (8) and is diagonalized to the mass eigenstate basis through Eqs. (A2)-(A5). We also work in the alignment limit where β − α = π/2 such that couplings to the light eigenstate h are SM-like. For singly charged leptons in a given model, the Lagrangian for Yukawa couplings to neutral Higgs bosons is defined by where where Y E and Y A E are Yukawa matrices defined per representation for CP -even (h, H) and CP -odd (A) bosons, respectively. The notations of Eq. (B37) are the same for all representations except for their respective diagonalization and Yukawa matrices.
Yukawa couplings of the charged Higgs boson H ± to leptons are defined by The last two terms are relevant for representations with doubly charged leptons. The couplings are defined as The matrices Y H ± E ′ ,N ′ are the Yukawa matrices present for doubly charged fermions coupling to singly charged ones.
For the 2 −1/2 ⊕ 1 −1 model, the Yukawa matrices for neutral Higgs bosons are The Yukawa matrices for charged Higgs bosons are In the 2 −1/2 ⊕ 3 −1 representation, the Yukawa matrices for neutral Higgs bosons are The Yukawa matrices for charged Higgs bosons are Then, for the 2 −3/2 ⊕ 1 −1 , we have the following Yukawa matrices for neutral Higgs bosons: The Yukawa matrices for charged Higgs bosons are The Yukawa matrices for neutral Higgs bosons in the 2 −3/2 ⊕ 3 −1 model are The Yukawa matrices for charged Higgs bosons are For the last model we consider, 2 −1/2 ⊕ 3 0 , we find that the Yukawa matrices involving neutral Higgs bosons are The Yukawa matrices for the charged Higgs are Appendix C: Approximate formulas for diagonalization matrices and couplings We present in full detail approximate diagonalization matrices for our main representation, as defined in Eqs. (B3) and (B4), and the other four models. We also give a complete list of approximate couplings to bosons in our representation. Computing approximate couplings in other representations is straightforward and follows similar arguments as those presented for our model. and Note that exact couplings of fermions to Z, W ± , H, h, A, H, and H ± bosons in the mass eigenstate basis were defined first in Appendices A1 and A2 of [25]. Those definitions apply identically in this paper, except that one must replace the diagonalization matrices with their complex versions U e,ν L and U e,ν R . In this limit, couplings of charged fermions to the Z boson are approximated by Couplings of charged and neutral fermions to the W ± boson are approximated below: Approximate couplings of charged fermions to the light SM Higgs boson h are given below: Couplings to the CP -even heavy Higgs field H can be found by simply replacing cos β → sin β in the above couplings to h. Couplings to the CP -odd heavy Higgs field A are given below: Finally, the approximate couplings of fermions to the charged Higgs field H ± are given by In the same expansion limit as for the main model above, the diagonalization matrices become and The diagonalization matrices for this model are the same as in the 2 −1/2 ⊕ 1 −1 model, Eqs. (C1)-(C3), with the exception that λ → −λ.
In the same mass limit, the diagonalization matrices for this representation are and The diagonalization matrices for doubly charged fermions are as well as Evaluating the same mass limit as above, the diagonalization matrices are The diagonalization matrices for the neutral sector are Appendix D: Contributions to the dipole moments in the Higgs basis In this appendix, we list the full contributions of C µB and C µW in the context of SU (2) singlets, doublets, and triplets discussed in [7]. For compactness, we use the shorthand notation cos β ≡ c β and sin β ≡ s β . Loop functions can be found in Eqs. (D20)-(D26) in this appendix and their useful limits can be found in Eqs (D27)-(D33).
For the main representation of the paper, 2 −1/2 ⊕ 1 −1 , and translating the Lagrangian of Eq. (1) to the Higgs basis while omitting up-type couplings, the Wilson coefficients are and where Y L,E are the hypercharges of the new fermions and x (1,2) L,E = M 2 L,E /M 2 1,2 . In the limit where M 2 L,E ≫ M 2 1 and M 2 L,E ≃ M 2 2 while inserting their respective hypercharges, we find If we consider the limit where M 2 L,E ≫ M 2 1 and M 2 2 , we find the same expression. Next, the Wilson coefficients for the 2 −1/2 ⊕ 3 −1 representation, using the Lagrangian of Eq. (B10) in the Higgs basis, are and Taking the same mass limits x L,E → ∞ and x L,E → 1 while inserting respective hypercharges, we find Considering the other limit where x (1,2) L,E → ∞, we have For the 2 −3/2 ⊕ 1 −1 representation, using the Lagrangian of Eq. (B15) in the Higgs basis, the Wilson coefficients are L ) L , x L , x L , x L ) and L , x L , x E ) + C(x (2) L , x L , x L , x L , x E ) .
Likewise, using the Lagrangian of Eq. (B19) in the Higgs basis for the 2 −3/2 ⊕ 3 −1 representation, the Wilson coefficients are L ) L , x L , x L , x L ) L , x L , x L ) L , x L , x L ) L ) L , x L , x L , x Taking the limits x In the limit where x (1,2) L,E → ∞, we have C µγ = 5 λ L λ Eλ 64π 2 M L M E ec 2 β 1 + tan 2 β .
Finally, for the last representation 2 −1/2 ⊕ 3 0 , we translate the Lagrangian of Eq. (B23) to the Higgs basis and compute the following Wilson coefficients: L , x L , x L , x L ) L , x L , x L , x L ) L ) L , x L , x L , x L,E → ∞, we find Note that for this representation m LE µ = −2 × λ L λ Eλ v 3 d /M L M E when rewriting C µγ in terms of C µH 1 .