The muon magnetic moment in the 2HDM: complete two-loop result

We study the ${\text{2HDM}}$ contribution to the muon anomalous magnetic moment $a_\mu$ and present the complete two-loop result, particularly for the bosonic contribution. We focus on the Aligned ${\text{2HDM}}$, which has general Yukawa couplings and contains the type I, II, X, Y models as special cases. The result is expressed with physical parameters: three Higgs boson masses, Yukawa couplings, two mixing angles, and one quartic potential parameter. We show that the result can be split into several parts, each of which has a simple parameter dependence, and we document their general behavior. Taking into account constraints on parameters, we find that the full ${\text{2HDM}}$ contribution to $a_\mu$ can accommodate the current experimental value, and the complete two-loop bosonic contribution can amount to $(2\cdots 4)\times 10^{-10}$, more than the future experimental uncertainty.


Introduction
The measured properties of the Higgs boson discovered at the LHC [1,2] are compatible with the Standard Model (SM) [3]. However, there is room for alternative explanations of the Higgs boson and electroweak symmetry breaking in models with extended Higgs sectors. The two-Higgs doublet model (2HDM) is a particularly interesting framework to be studied. In a large part of its parameter space it is compatible with experimental data, it can originate from more fundamental theories like the MSSM, and it predicts a multitude of observable effects by which it can be studied and constrained. 1 Here we focus on the muon anomalous magnetic moment a µ = (g µ − 2)/2 in the 2HDM. This is one of the most useful precision observables to provide complementary, non-collider constraints of extensions of the SM [4][5][6]. After significant recent progress on all aspects of the SM prediction, there is a stable 3-4 σ deviation between the SM prediction and the Brookhaven measurement [7], using the indicated references for the hadronic vacuum polarization contributions. 1 Several recent studies [28][29][30][31][32][33][34] have shown that the 2HDM has viable parameter regions in which this (or at least most of this) deviation is explained. The existing studies evaluate a µ in the 2HDM using one-loop and particular two-loop diagrams, so-called Barr-Zee diagrams. Such Barr-Zee diagrams were first considered in Ref. [35] and for a µ in the 2HDM in Refs. [36][37][38][39]; the most complete calculation is presented in Ref. [30]. Here we present and document the full two-loop calculation of a µ in the 2HDM, including Barr-Zee and non-Barr-Zee diagrams.
Our calculation is motivated in two ways. Firstly, the 2HDM one-loop contributions are suppressed by two additional powers of the small muon Yukawa coupling. Thus the one-loop contributions are parametrically smaller than the two-loop contributions. In this sense our calculation completes the leading-order prediction of a 2HDM µ . Secondly, new a µ experiments are planned at Fermilab and JPARC [40,41]. These promise to reduce the experimental uncertainty significantly, in particular the Fermilab measurement plans to obtain ∆a Fermilab µ = 1.6 × 10 −10 . ( This highlights the need for reliable and accurate theory predictions also in extensions of the SM. In the electroweak SM, the full two-loop calculation has been done in Refs. [11,[42][43][44]. In other models, such as the MSSM, several classes of two-loop contributions have been evaluated [44][45][46][47][48][49][50][51]. It has been found that each class can give rise to significant corrections, and an analysis of the remaining MSSM theory uncertainty has shown that the future experimental precision can only be matched by a complete two-loop computation [52,53]. This paper is divided as follows: in Sec. 2 we review the 2HDM and introduce the phenomenological constraints adopted in our analysis. In Sec. 3 the complete renormalized 2HDM two-loop contributions to a µ is presented. Each part of the computation is documented in a series of plots and/or analytic formulas. We perform a numerical analysis of our result in Sec. 4, showing that the complete two-loop bosonic contribution can amount to (2 · · · 4) × 10 −10 , i.e. at the level of the precision of the planned Fermilab experiment. We present our conclusion in Sec. 5. Appendix A contains all analytic formulas of the renormalized bosonic two-loop contributions to a µ while in Appendix B we discuss the cancellation of M A dependence in Y A l sector.

The model and its parameters
The two-Higgs-Doublet Model (2HDM) is an extended Standard Model (SM) with two complex scalar doublets Both scalar doublets are assigned with the same hypercharge as the SM doublet. The vacuum expectation value (VEV) of the SM, v, is recovered by the relation v 2 = v 2 1 + v 2 2 . The most general form of the Higgs potential V (φ 1 , φ 2 ) depends on eleven physical parameters [54]. In this work, we consider the CP-conserving case in which all parameters are real. For simplicity, we also impose an approximate Z 2 symmetry which restricts the quartic couplings to five (denoted by λ 1,2,3,4,5 ) while the quadratic couplings are given by three parameters (m 2 11 , m 2 22 , and m 2 12 ), the last one breaking the Z 2 symmetry softly [54,55] V (φ 1 , φ 2 ) =m 2 11 Through a rotation with angle tan β ≡ t β ≡ v 2 /v 1 , we can choose new scalar In the new basis only the doublet Φ v contains the VEV and the Goldstone bosons, and the components are explicitly H ± corresponds to the charged Higgs bosons and A to the neutral CP-odd one. S 1 and S 2 are not mass eigenstates, but they are related to the CP-even neutral mass eigenstates h, H through a new mixing angle α as If β −α = π 2 , the two mass eigenstates are completely separated in each scalar doublet and the neutral CP-even Higgs boson h has just the same interactions as the SM Higgs boson, h SM . We call this the SM-limit, following Ref. [55]. The LHC data allow a small deviation η [56], which we define as In this work we present the results away from the SM-limit, where η = 0. Seven of the eight parameters, m 2 11 , m 2 22 , m 2 12 , λ 1 , · · · , λ 5 , introduced in the 2HDM potential Eq. (4) can be replaced with physical parameters such as the scalar boson masses, M h , M H , M A , M H ± , the mixing angles, β, α, the VEV, v [54,55]. The tree-level relations for the λ i can be written as where c α = cos α, s α = sin α, c β = cos β, and s β = sin β. We are still left with one more free parameter m 2 12 , or equivalently λ 1 . It is convenient to define the quantity Λ 5 , which absorbs m 2 12 or λ 1 , as The equivalent relation in terms of λ 1 can be written up to η 1 -order as 2 All the previous relations hold at tree level and might be modified at higher orders, depending on the chosen renormalization scheme for the Higgs sector parameters. Renormalization schemes for the 2HDM Higgs sector parameters have been discussed recently in Refs. [59,60]. For our purposes it will turn out that the tree-level relations are sufficient. We complete the discussion of the 2HDM by introducing the fermionic sector. The Yukawa coupling is model-dependent. In the present paper we focus on the Aligned 2HDM. The Aligned 2HDM is very general and contains the usual type I, II, X and Y models as special cases: see Table 1.
In the Aligned 2HDM it is only required that the mass matrices and the Yukawa coupling matrices in the most general Yukawa Lagrangian are 2 Λ 5 corresponds to λ 5 in the 2HDM model file of FeynArts [57].
proportional to each other with proportionality constant, ζ f [58]. The aligned Yukawa Lagrangian reads where P R,L = 1 2 (1 ± γ 5 ), and V CKM is the Cabibbo-Kobayashi-Maskawa matrix. The Yukawa coupling matrices are defined as where M f denotes the diagonal 3×3 fermion mass matrix. We have f = u, d, l and S ∈ {h, H, A, H ± }. The generation independent coefficients Y S f are specific for each model.
In the Aligned 2HDM, Y S f are dependent on (β − α) and ζ f , and we have [58] Since we focus on small deviations from the SM-limit, i.e. small η, it is useful to expand the coefficients of Eq. (18) for small η, The parameters ζ l,u,d are constrained by experimental results of other physical processes. The detailed explanation on the allowed parameter regions is given in Sec. 2.2. Types I, II, X, Y are recovered by assigning specific values of the aligned parameters ζ f as listed in Table 1.

Constraints
Following the presentation of Ref. [28], we introduce some constraints to restrict the allowed parameter region. They are mainly theoretical and electroweak (EW) constraints. As theoretical constraints, we consider the requirements of stability, and perturbativity that the scalar potential must

Theoretical constraints
The theoretical constraints faced by the 2HDM are of two different natures. The first is related to the stability of the potential, requiring that a vacuum minimum exists and that this minimum is the global minimum of the system. The second is related to perturbativity, requiring that none of the couplings exceeds a given maximal value. For the CP-conserving potential Eq. (4), all these requirements are translated into relations between the different λ i introduced on the potential as below [54,61,62]: • Stability • Global minimum • Perturbativity As Ref. [28,61], we adopt λ max = 4π. In the phenomenological analysis we employ Eqs. (9) - (13) to translate the constraints of Eqs. (20) -(22) into 7 those on the physical mass parameters. Since we do not assume the 2HDM to be necessarily a fundamental theory valid up to very high energy scales, we require the validity of the above conditions only for the tree-level parameters. For constraints from requiring conditions on running, high-scale parameters see particularly Ref. [61].

Electroweak and experimental constraints
Regarding electroweak precision data, we will include the constraints on the Peskin-Takeuchi parameters S, T and U [63,64] S = −0.03 ± 0.10, T = 0.01 ± 0.12, To implement them in our phenomenological analysis, we use 2HDMC-1.7.0 [65,66] to restrict the allowed parameter space on the masses of the scalars. We also include the model-independent constraint obtained by LEP on the mass of the charged scalar [64] M H ± ≥ 80 GeV.
Finally we introduce the constraints on the aligned parameters ζ f . As discussed in [56], in order to avoid conflict with current LHC data they should satisfy

2HDM Two-loop Contributions
The purpose of our study is to present the complete two-loop 2HDM contribution to a µ . The renormalized two-loop result a 2HDM,2 µ is the sum of the one-loop contribution a 2HDM,1 µ , two-loop bosonic and fermionic contributions and a shift from using the Fermi constant The actual renormalized two-loop contributions, a B µ and a F µ , are obtained from the sum of the appropriate two-loop and one-loop counterterm diagrams. The one-loop contribution and a ∆r-shift In the EW SM, it is sufficient to evaluate the full result only up to order m 2 µ /M 2 W and neglect higher order terms of O(m 4 µ ). In the 2HDM, however, there are potentially non-negligible terms of this order. Hence we evaluate a 2HDM,1 µ up to O(m 4 µ ), but at the two-loop level terms up to O(m 2 µ ) are sufficient. We furthermore expand the results in the small parameter η = α − β + π/2 up to the order η, and we set the mass of the Higgs boson h to the mass of the observed Higgs boson, Our calculational procedure is based on the one described in Refs. [44,45] using TwoCalc [67] for evaluating two-loop integrals and in-house routines for reduction to master integrals, large mass expansion, and analytical simplification.

One-loop contribution
The 2HDM one-loop result is expressed as [68][69][70] a 2HDM,1 where S ∈ {h, H, A, H ± }. G F is the muon decay constant. The Y S l are given in Eq. (19).
For each Higgs boson the loop-function F S is defined as The right column shows the approximations in the small x limit [28].
The numerical values and signs of the different contributions can be easily read off from the following approximation, usingx S ≡ M S /100 GeV, and neglecting terms of order η, a 2HDM,1 µ ζ l 100 9 Mass renormalization constants: Tadpole renormalization constants: δT h , δT H At this point we also remark that the EW SM one-loop result is evaluated in terms of the muon decay constant G F , G F is related to the input parameter of the on-shell renormalization scheme as As a result of this, if the on-shell scheme is used to define the counterterms for the two-loop calculation, there is an additional contribution a EW(1) µ × (−∆r); see also [11,[42][43][44].
The extra contribution for the 2HDM is then given by where ∆r 2HDM is the extra 2HDM contribution to ∆r. It is discussed in detail in [71,72]. In accordance with Ref. [72] we verified that in all the parameter space relevant for our analysis ∆r 2HDM is at most of the order of 10 −3 , and thus |a ∆r-shift µ | ≤ 2 × 10 −12 .

Counterterm contribution
The 2HDM counterterm diagrams involve the renormalization constants in Table 2. These renormalization constants are defined in the on-shell renormalization scheme [73,74]. In terms of these the electric charge renormalization constant δZ e is derived as where δZ AA is the photon field renormalization constant and δZ ZA the photon-Z mixing renormalization constant. For the mass and field renormalization constants several useful statements can be made. They are obtained from self-energy diagrams with external SM particles. In the expansion in η up to O(η 1 ) each mass and field renormalization constant can be decomposed into the SM and additional 2HDM contributions. These additional 2HDM contributions to the renormalization constants are obtained by computing the loop diagrams containing the new scalar bosons in the 2HDM. For these renormalization constants the fermionic contributions are the same in both the SM and 2HDM. Therefore the additional 2HDM contributions from these diagrams arise entirely from the bosonic parts. The tadpole renormalization constants should be treated separately. The tadpole renormalization constants are determined in such a way that the one-point Green functions of the CP-even Higgs fields vanish. In the CPconserving 2HDM there are two tadpole renormalization constants, δT h and δT H , whereas in the SM we have one tadpole renormalization constant δT SM . The contributions from gauge(g) and Goldstone(G) bosons and fermions(f ) are related to the SM counterpart by simple rescaling of couplings as However, the Higgs loops of the tadpole diagrams are proportional to the triple Higgs couplings, are t β -dependent, and do not satisfy such a simple relation.
We now turn to the counterterm diagrams. Fig. 1 shows sample diagrams. It is convenient to classify the 2HDM counterterm diagrams into three groups. The first group encompasses the SM-like counterterm diagrams without Higgs boson inside. The second group contains the counterterm diagrams with neutral physical Higgs bosons. The third group consists of the counterterm diagram with G ± -H ± mixing counterterm vertex. In the following we explain them one after the other and provide the explicit results.
• SM-like counterterms without physical Higgs bosons: The first group encompasses the SM counterterm diagrams which do not contain physical Higgs bosons. The results of SM counterterm diagrams are found in Ref. [42]. The additional 2HDM contributions from these counterterm diagrams are obtained by applying the corresponding additional 2HDM renormalization constants. This is straightforward for all diagrams except Fig. 1a. This diagram is the only one which contains the tadpole renormalization constants. In the following we explain the cancellation of the gauge and Goldstone boson contributions as well as the fermion contributions.
The additional 2HDM contribution from this diagram is the difference between the 2HDM and SM results, The counterterm vertices of the G ± -G ± propagator for the SM and 2HDM are Hence, from Eqs. (36) and (39) we find that the tadpole renormalization constants with gauge/Goldstone bosons and fermion loops drop out from the result, Eq. (37). Consequently no new 2HDM contributions are obtained from the fermion, gauge boson, and Goldstone boson loops. The additional 2HDM contribution of Fig. 1a and all other diagrams of this group arises from the physical Higgs boson loop contributions to the tadpole, field and mass renormalization constants.
• Counterterm with neutral Higgs bosons: The second group of counterterm diagrams is shown in Fig. 1b. These diagrams are dependent on Yukawa coupling and proportional to δZ ZA , which is the same in both the SM and the 2HDM. The difference arises from the Yukawa coupling and the second neutral CP-even 2HDM Higgs boson. The fermionic loop contribution to δZ ZA are zero, therefore the diagrams of this group do not contribute to the fermion contributions.
The additional 2HDM contribution is obtained when the SM Higgs boson contribution is subtracted from the 2HDM results, a , where the explicit result of Fig. 1b for an arbitrary scalar field S is S can be any of the neutral Higgs bosons in the SM and the 2HDM: The contribution of the CP-odd neutral Higgs A is zero.
The coefficient, C S , is derived from the gauge coupling. It is 1 for the SM Higgs h SM , sin(β − α) for h and cos(β − α) for H. Y S l is derived from the Yukawa coupling constant and listed in Eq. (19). For the SM Higgs h SM , Y h SM l = 1.
• Counterterm diagram with G ± -H ± mixing: The third group consists of the diagram of Fig. 1c which is proportional to the G ± -H ± mixing counterterm vertex. This counterterm diagram does not appear in the SM. The explicit analytic result reads where δt HG = cos(β − α) δT h − sin(β − α)δT H . ζ l -dependency arises from the charged Higgs boson coupling to the muon in the Aligned   [30,39,[45][46][47][48]. The 2-boson diagrams also contain self-energy type diagrams in which the external photon couples to the muon line. The 3-boson diagrams have a more complicated structure and involve three internal bosons which couple to the muon line. Fig. 3 shows all 3boson diagrams which contribute to the difference between the 2HDM and the SM. In addition, diagrams with four bosons coupling to the muon line exist but do not contribute to the difference between the SM and 2HDM at the order of O(m 2 µ ). We especially computed the 3-boson diagrams shown in Fig. 3 for the first time. Figs. 3a and 3b are dependent on the W boson and Figs. 3c and 3d on the Z boson. While computing the diagrams in Fig. 3, we should pay attention to two interactions. One is the muon Yukawa coupling to the neutral scalar bosons, h or H, and the other is the Higgs-gauge interaction of the two neutral Higgs bosons. The gauge interaction to H is suppressed by η. In the SM-limit the contribution from this interaction becomes zero, whereas the gauge interaction to h recovers the SM value. The explicit result of Figs. 3a and 3b reads where S ∈ h, H, and y S ≡ and x S ≡

Analytic results
In this section we present the complete renormalized bosonic 2HDM contribution. The bosonic result is expanded with respect to the parameter η introduced in Sec. 2, and terms up to η 1 are taken. In the SM-limit, η → 0, the interactions of h to the gauge bosons or fermions become just those of the SM. For the discussion of the complete result we do not use the 2-boson and 3-boson separation. Instead, we divide the renormalized bosonic contribution into two parts. One part, a EW add. µ , is defined by the Feynman diagrams containing only gauge/Goldstone/h bosons, i.e. purely SM-like diagrams. The other part is defined by those diagrams which include at least one of the new 2HDM Higgs bosons, H, A, H ± . This part can be again divided into Yukawadependent and Yukawa-independent parts. Considering this classification we can write the bosonic contribution as where a non-Yuk µ denotes Yukawa-independent 2HDM Higgs contributing part and a Yuk µ the Yukawa-dependent part. In the following we explain each of the contributions explicitly.
• a EW add.
The sign of a EW add.
µ is dependent on η ζ l . Even though η must be small, the appearance of ζ l can enhance the contribution of a EW add.    Fig. 3 with H contribute to a Yuk µ , too. Clearly, all diagrams with H or H ± and gauge bosons are suppressed by η but enhanced by ζ l . The diagrams without gauge bosons involve triple Higgs couplings and are of particular interest. A closer look at the triple Higgs coupling constants helps to analyze the t β -dependency. The triple Higgs couplings constants in the 2HDM are The triple Higgs coupling constants show that the t β -dependency comes only in the form of (t β − 1 t β ), which leads to a large t β -enhancement. In the actual Feynman diagrams with triple Higgs couplings, the coupling Eq. (50) appears multiplied with Y h l , and the coupling Eq. (51) is The notation is such that the terms with superscript 0 are independent of η, the terms with superscript 1 are linear in η. The subscript z denotes terms enhanced by ζ l , the subscript 5 denotes terms ∝ Λ 5 . All terms here arise from diagrams with triple Higgs couplings except the a 1 0,z term. The results of the 3-boson diagrams Eqs. (42) and (44)  The plots in Fig. 5 show the complete mass dependence of the coefficients a 0 i,j . a 0 0,0 and a 0 5,0 arise from the Feynman diagrams containing the muon Yukawa interaction to h and the η-independent part of Eq. (50), therefore are dependent only on M H ± and neither enhanced by t β nor by ζ l . In contrast, a 0 0,z and a 0 5,z arise from diagrams involving the triple Higgs coupling Eq. (51) and appear enhanced by large t β and ζ l in Eq. (52).
The plots in Fig. 6 show the change of a 1 0,0 , a 1 0,z , a 1 5,0 and a 1 5,z , the ηsuppressed terms. The coefficient a 1 0,z , which gets contributions from a larger class of diagrams, can be numerically larger than the other coefficients.

Fermionic loop contribution
In this section we present the fermionic loop contribution to a µ . Due to the higher order muon mass suppression (considering terms up to m 2 µ order),    all diagrams contain only one scalar boson, which interacts with the incoming/outgoing muon and the fermion in the inner loop. Thus, the result is always proportional to the product of two Yukawa couplings Y S l Y S f . The fermionic two-loop Feynman diagrams contain either neutral or charged Higgs bosons. Fig. 7a shows the generic diagrams for neutral Higgs bosons while Fig. 7b is related to charged bosons. When the external photon couples with the muon line we obtain self-energy type diagrams, and the sum of these vanishes. The remaining diagrams are Barr-Zee diagrams.
Our result for neutral Higgs bosons is coincident with previous analysis 3 [30,36,37], and the explicit form is where For S = {h, H} we have and for S = A A sum over all types of fermions is implicit. Q f denotes the charge of the respective fermion f , and N f c the color factor. We also define g f v ≡ T 3 2 − Q f s 2 W , and Φ(M S , m f , m f ) is defined in Appendix A. Both γ and Z bosons contribute to the fermionic loop result with neutral Higgs bosons. However, the result from the Z boson is suppressed by factor g f v , which is −1/4 + s 2 W ∼ −0.02 for leptons, compared to the result from the diagrams with photon. Hence the Z contributions are always smaller than those of the photon. Now we turn to the fermionic two-loop contributions with charged Higgs bosons. Figs. 7b and 7c show the corresponding Feynman diagrams. Especially the result of Fig. 7c is divergent, the corresponding counterterm diagram is shown in Fig. 1c. The renormalized two-loop result is obtained by summing up the two-loop and the counterterm diagrams. These diagrams were computed in the context of SUSY models long ago [46,47] in which case a type II structure for the Yukawas needed to be assumed. In the case of a general model (Aligned Model, for instance) the analysis was only recently performed [30]. We also recover the analytic result presented in [30,46,47], explicitly where M f corresponds to pairs of fermions masses as 0), (m µ , 0), (m τ , 0)}, and Eq. (58) contains an implicit sum over pairs. We neglect neutrino masses and generation mixing. We also introduced the definitions below in Eq. (58) Summing the results of Eqs. (53) and (58) and subtracting the corresponding SM-Higgs contribution gives the full renormalized two-loop 2HDM fermionic contribution where (64) is dependent on only one mass parameter, M S , and this enables us to analyze the individual Higgs boson contributions to the fermionic loop contribution in Fig. 8. The first line of Eq. (64) contains terms bilinear in the ζ f , and they are shown in the first three plots of Fig. 8. The terms in the second line are proportional to ζ f η and are illustrated in the fourth plot, Fig. 8d.
In all cases the contribution from the top loop (blue line) is significantly larger, as expected by the factor m 2 t /(M 2 S − M 2 B ) in the analytic formulas (M B is the mass of the internal gauge boson involved). However, as discussed in Sec. 2.2.2, ζ u is constrained to be at most ζ u 1, meaning that the tau loop, enhanced by ζ 2 l , plays the decisive role. Another characteristic shared by Figs. 8a -8c is that they all decrease with the mass of the scalars. Fig. 8d shows the contribution proportional to η which comes from diagrams involving CP-even scalar bosons. As presented in Eq. (64) there is a difference between the h and H results, explaining why the η 1 contribution vanishes as M h approaches M H . For all plots, we have rescaled a µ to the aligned parameters. Finally, in all graphs the contributions can be both positive or negative. The signs depend on the alignment parameters and can be read off from Table 3. Table 3: Relation between signs of the aligned parameters and the functions depicted in Fig. 8.

Numerical Analysis
In this section we present the numerical analysis of our result. Our aim is to study how large the bosonic contribution, fully computed for the first time, can be. We will show that there are regions of the parameter space in which a B µ amounts to (2 · · · 4) × 10 −10 . Although always smaller than the fermionic contribution, it proves to be relevant for a precise determination of the 2HDM contribution to the muon anomalous magnetic moment. We also analyze the impact of deviations from the SM-limit by studying different values for the expansion parameter η.
Since we want to study the impact of the SM-limit deviation to a µ , hereafter we will choose specific values for η and compare how the results differ. We perform a scan over the above region, computing for each point the value of the full a µ as well the contribution only due to two-loop bosonic Feynman diagrams. Our results of the full scan are depicted as blue points in the plots of Fig. 9 (for the three values η = 0, η = 0.1, η = −0.1). We then apply the further experimental/theoretical constraints discussed in Sec. 2.2. The survival sample is depicted as red points in Fig. 9.
As can be readily seen, although the values for the full a µ can be large, the contribution from a B µ can amount to (2 · · · 4) × 10 −10 . One can also notice a difference in behavior between the SM-limit case and the one in which η is negative. In the latter case, one observes that the range of values for a B µ is significantly larger, spreading over the x-axis, while in the former it is constrained inside the region with absolute value 2 × 10 −10 .
In order to obtain a better insight into the bosonic contribution, we choose a sample point for which the muon anomaly can be explained and vary the parameters affecting mainly a B µ . For comparison with the previous analysis [28], we consider as starting value a parameter space point allowed by the type X model.
In the type X model, the explanation to the a µ deviation comes mainly from fermionic contributions containing a tau loop. The reason is that, in this model, only the Yukawa coupling of leptons is enhanced, see Table 1. In the Aligned Model, a type X scenario is recovered if ζ u = ζ d = 1/t β , and ζ l is identified as −t β . In Ref. [28], it was found that the anomaly could be explained for low values of the CP-odd scalar mass (M A < 100 GeV), large t β , and values of the masses of the CP-even and charged scalar of the order of 200 GeV. In that reference, the type X model was considered and only fermionic contributions were included. Therefore, it is particularly simple to translate parameter points to the Aligned 2HDM, by identifying t β as −ζ l . After these considerations, we choose as representative point the one defined by 4 In the Aligned model, the values for t β , λ 1 , and η remain free. The first two are only related to the bosonic contribution via triple Higgs couplings, η affects the bosonic and fermionic contributions. Fig. 13 shows the results from varying these three parameters, and thus particularly the impact of the bosonic contribution to a µ . In all plots we depict a 2HDM,2 µ on the upper graph, and a B µ on the lower one. The upper plots contain as a reference line the value for a µ used in [28], which takes into account only fermionic contributions for η = 0. The η-dependence is depicted by red lines (η = 0), blue lines (η = 0.1), and green lines (η = −0.1). We proceed to explain each of the graphs individually. Figs. 13a -13b show the behavior as a function of t β for λ 1 = 4π. 5 As expected from the scatter plots, the variation of a B µ is in the range (2 · · · 4) × 4 It should be noticed that any other point considered in [28] for which a µ is explained at 1 σ level could be chosen as well. The behavior of all plots as well as all further discussions remain essentially the same. Furthermore the recent references [31,34] also considered τ -decay as a parameter constraint in the 2HDM, which disfavors a significant part of the preferred parameter space. Nevertheless, reference [34] found the general parameter region represented by Eq. (66) to be viable. 5 The analysis is unaltered for other choices to λ 1 , only the absolute value of the bosonic contribution is modified.  There are two regions: small t β and large t β . For small t β , Λ 5 is dominated by the negative term proportional to λ 1 and several bosonic contributions are suppressed by (t β − 1/t β ) which vanishes as t β → 1. This explains the linear behavior in Fig. 13c and the peak in Fig. 13b.
For large t β , Λ 5 2M 2 H /v 2 1.32, and the prefactor (t β − 1/t β ) t β . This explains the linear behavior of the contributions for t β > 20 in Fig. 13a and the independence of λ 1 in Fig. 13d.
Regarding the η-dependence, the dominant terms depending on η are a EW add. µ , Eq. (49), and a 1 0,z , Fig. 6b. For the present parameter region, the coefficients of η ζ l are approximately (2.3 − 1) × 10 −11 . This explains that shifting η by 0.1 decreases a B µ by 10 −10 in all plots. In order to compare our analysis of Fig. 13 with the scatter plots of Fig. 9, we show three representative points in Fig. 9. The first, in green, is the representative point just discussed for the large t β regime (t β = 100, λ 1 arbitrary). The other two, yellow and dark yellow, are related to the low t β regime and have t β = 2 and two different values of λ 1 , λ 1 = 4π (yellow) and λ 1 = 2π (dark yellow). As can be seen, the green triangle for η = −0.1 and η = 0 is close to the border of the constrained sample depicted in red while, for η = 0.1, the green triangle is well inside the allowed area. It is instructive to notice that, for negative values of η, there is a considerable sample of allowed points with similar values for a B µ (2 · · · 3) × 10 −10 . This behavior is explained by observing Figs. 13a -13d which show that for any value of λ 1 there is a large interval for t β , 40 < t β < 100, allowing a B µ > 2 × 10 −10 . This situation should be contrasted with the low t β regime, represented by the yellow triangles. While the η = −0.1 case still has a considerable amount of points with similar values, the (η = 0, 0.1) cases represent rare points in the constrained sample for λ 1 = 4π, and points close to the border of the allowed area for λ 1 = 2π. The explanation can be found in Figs. 13b -13c which show that values for a B µ similar to the ones of the light yellow triangle can only be obtained for a small range of t β , 1 < t β < 5, and large values of λ 1 , λ 1 4π. These observations explain why the scatter plot for negative η, Fig. 9c, has more allowed points with values for |a B µ | of order (2 · · · 4) × 10 −10 if compared with the other cases.
Finally we discuss the plots of Fig. 14. In both cases we study the behavior of a 2HDM,2 µ and a B µ as functions of one of the masses of the scalars (M H ± , and M H respectively) where the region delimited by the dashed lines is allowed by theoretical and EW constraints. The other mass and aligned parameters are kept fixed as in the representative point Eq. (66). Regarding t β , we choose t β = 100, corresponding to a type X parameter point. Since we are in the large t β limit, λ 1 has no significant influence. We adopt λ 1 = 4π. As can be seen in Fig. 14a, there is a slight mass dependence in a 2HDM,2 µ . To illustrate the mass dependence we first remark that, for the parameter region we are considering, only the coefficients enhanced by ζ l and t β are important, namely a 0 0,z and a 0 5,z . Using the definitions in Appendix A and considering the large t β region, one has where x S ≡ M 2 S /M 2 Z , and we used Λ 5 2M 2 H /v 2 . The term containing the functions F 0 m , F ± m is always positive and depends on the inverse of the scalar masses.
Therefore, if M H is kept fixed and M H ± increases, |a B µ | will decrease, explaining the behavior observed in Fig. 14a. In contrast, if M H ± is kept fixed and M H increases, the explicit dependency on M 2 H coming from Λ 5 and the coefficient b(x H , 0) leads to an increase of a B µ with M H in Fig. 14b. Regarding the full a µ , we verified that the fermionic contributions essentially do not depend on M H ± due to the small ζ u , but they depend on M H . As can be noticed analyzing the plots of Fig. 8 and Table 3, the fermionic contributions from H diagrams are negative and decrease in modulus with M H . Therefore, the net result will be an increase in a µ as observed in Fig. 14b.
Finally, it can also be noticed that the plots for non-zero values of η tend to the case η = 0 as M H approaches M h . This behavior can be understood by observing that in this case the two mass-degenerate CP-even scalars together behave exactly SM-like.

Conclusion
We presented the full two-loop 2HDM contributions to a µ , providing the complete analytic result and a numerical analysis. We confirmed the previous results of the fermion-loop and the bosonic Barr-Zee type contributions. We calculated the remaining diagrams including all 3-boson diagrams, which involve three internal boson couplings to the muon line.
The analytic results are expressed in terms of physical parameters. The full bosonic result depends on the three additional Higgs boson masses, t β , sin α, the alignment parameter ζ l and the quartic scalar coupling λ 1 . We always expand in the small parameter η = α − β + π/2, the deviation from the SM-limit. The bosonic contributions are especially dependent on t β and λ 1 , whereas fermionic ones are not. This dependency arises from the triple Higgs couplings in the bosonic Feynman diagrams.
We split the bosonic result into several parts, see Eq. (48) and Eq. (52). Each term has a straightforward dependence on t β and ζ l and depends only on a subset of masses. The compact analytic expression of each term is provided in Appendix A. We documented the parameter dependence in a series of Figures in Sec. 3

.3.
We also confirmed the previous result of the fermionic contribution. Particularly, we presented its analytic form without one-dimensional integral relations in Sec. 3.4 and gave an overview of the numerical behavior. The fermionic result involves all three alignment parameters ζ l,u,d , but the leading contributions are the ζ l dependent terms.
We also investigated the impact of the scenario with a deviation from the SM-limit of the Higgs couplings, η = α − β + π/2 = 0. For this case, we obtain additional contributions from the SM-like Higgs boson, a EW add. µ . This term is proportional to η ζ l and gives the dominant η-dependent bosonic contributions. Its coefficient is dependent only on the SM parameters and can be found in Eq. (49).
In the numerical evaluation we confirmed that the fermionic 2HDM contribution can be of the order of the deviation Eq. (1). A series of plots shows that in parameter regions with large fermionic contributions, the complete bosonic result can yield additional contributions in the range (2 · · · 4)×10 −10 , i.e. at the level of the precision of the planned Fermilab experiment. Allowing the SM-like Higgs couplings to deviate from the SM-limit, i.e. η = 0, and non-zero values of λ 1 , can slightly increase the bosonic contributions.
Here we provide revised versions of the plots in figures 10 and 11 with the modified numerical evaluation using tree-level relations between electroweak parameters: see figures 13 and 14. Here the linear increase in t β is absent and the numerical results in the large t β regime are slightly changed. Again, we refer to Ref. [1] for a more detailed phenomenological evaluation which also takes into account a variety of experimental constraints on the input parameters ζ l , t β , and η.