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

We study the 2HDM contribution to the muon anomalous magnetic moment aμ and present the complete two-loop result, particularly for the bosonic contribution. We focus on the Aligned 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 2HDM contribution to aμ can accommodate the current experimental value, and the complete two-loop bosonic contribution can amount to (2⋯4) × 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.
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 JHEP01(2017)007 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 [36][37][38][39][40][41][42] 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. [43] and for a µ in the 2HDM in refs. [44][45][46][47]; the most complete calculation is presented in ref. [38]. 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 [48][49][50]. These promise to reduce the experimental uncertainty significantly, in particular the Fermilab measurement plans to obtain ∆a Fermilab µ = 1.6 × 10 −10 . (

1.2)
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,[51][52][53][54].
In other models, such as the MSSM, several classes of two-loop contributions have been evaluated [54][55][56][57][58][59][60][61]. 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 [62,63]. This paper is divided as follows: in section 2 we review the 2HDM and introduce the phenomenological constraints adopted in our analysis. In section 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 section 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 section 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 , j = 1, 2. (2.1) 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 [64]. In this work, we consider the CP-conserving case in which all parameters are real. We also impose an approximate Z 2 symmetry which demands that two parameters, commonly called λ 6 and λ 7 , vanish, while the soft Z 2 symmetry breaking term m 2 12 is allowed. This 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 ) [64,65] V (φ 1 , φ 2 ) = m 2 11 Through a rotation with angle tan β ≡ t β ≡ v 2 /v 1 , we can choose new scalar doublets In the new basis only the doublet Φ v contains the VEV and the Goldstone bosons, and the components are explicitly .

(2.4)
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. [65]. The LHC data allow a small deviation η [66], which we define as In this work we present the results away from the SM-limit, where η = 0.

JHEP01(2017)007
Seven of the eight parameters, m 2 11 , m 2 22 , m 2 12 , λ 1 , · · · , λ 5 , introduced in the 2HDM potential eq. (2.2) 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 [64,65]. The tree-level relations for the λ i can be written as (2.7) 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. [69,70]. 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 proportional to each other with proportionality constant, ζ f [68]. The aligned Yukawa Lagrangian reads Sf y S f P R f + h.c., (2.14) Table 1. Relation between the Yukawa aligned parameters ζ f and the usual type I, II, X, and Y models.

JHEP01(2017)007
where P R,L = 1 2 (1 ± γ 5 ), and V CKM is the Cabibbo-Kobayashi-Maskawa matrix. The Yukawa coupling matrices are defined as Since we focus on small deviations from the SM-limit, i.e. small η, it is useful to expand the coefficients of eq. (2.16) 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 section 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. [36], 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 retain. Regarding EW constraints, we assure that the allowed range for masses of the new scalars does not violate the experimental measured values of EW precision observables such as M 2 W or sin θ W .

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

JHEP01(2017)007
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. (2.2), all these requirements are translated into relations between the different λ i introduced on the potential as below [64,71,72]: • Stability • Perturbativity As refs. [36,71], we adopt λ max = 4π. In the phenomenological analysis we employ eqs. (2.7)-(2.11) to translate the constraints of eqs. (2.18)-(2.20) into 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. [71].

Electroweak and experimental constraints
Regarding electroweak precision data, we will include the constraints on the Peskin-Takeuchi parameters S, T and U [73,74]  To implement them in our phenomenological analysis, we use 2HDMC-1.7.0 [75,76] 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 [74] M H ± ≥ 80 GeV.

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 JHEP01(2017)007 and a ∆r-shift µ are discussed in section 3.1, the counterterm parts in section 3.2. The bosonic and fermionic results are presented in section 3.3 and section 3.4, respectively.
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. [54,55] using TwoCalc [77] 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 [78][79][80] a 2HDM,1 where S ∈ {h, H, A, H ± }. G F is the muon decay constant. The Y S l are given in eq. (2.17). For each Higgs boson the loop-function F S is defined as The right column shows the approximations in the small x limit [36].
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 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 (3.8)

JHEP01(2017)007
Mass renormalization constants:  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,[51][52][53][54]. 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 [81,82].
In accordance with ref. [82] 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 [83,84]. 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 CP-conserving 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 (3.11) 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. Figure 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.
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. [51,52]. 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 figure 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. (3.11) and (3.14) we find that the tadpole renormalization constants with gauge/Goldstone bosons and fermion loops drop out from the result, eq. (3.12). Consequently no new 2HDM contributions are obtained from the fermion, gauge boson, and Goldstone boson loops. The additional 2HDM contribution of figure 1a and all other diagrams of this group arises from the physical Higgs boson loop contributions to the tadpole, field and mass renormalization constants.

JHEP01(2017)007
• Counterterm with neutral Higgs bosons. The second group of counterterm diagrams is shown in figure 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 figure 1b for an arbitrary scalar field S is • Counterterm diagram with G ± -H ± mixing.
The third group consists of the diagram of figure 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 2HDM. Figure 1c is the only counterterm diagram which contributes to the fermionic two-loop result. Figure 2. Generic 2-boson Feynman diagrams. The gray loops denote any bosonic loop.  Figure 3 shows all 3-boson 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 figure 3 for the first time. Figures 3a and 3b are dependent on the W boson and figures 3c and 3d on the Z boson. While computing the diagrams in figure 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 figures 3a and 3b reads

JHEP01(2017)007
where S ∈ h, H, and y S ≡ For the SM Higgs boson, Y h SM l = 1, and C h SM = 1. The divergent part of eq. (3.17) drops out in the final result of the difference of the SM and 2HDM. Note that the result of figure 3b alone is finite. In the off-SM scenario, η = 0, the result of figure 3 for S = h results in additional EW contributions.
The additional 2HDM contribution from the diagrams of figures 3a and 3b is obtained when the SM Higgs boson result of eq. (3.17) is subtracted from the sum of the h and H contributions, After employing the known SM parameters we obtain the numerical result . C S and Y S l for h and H are the same as in eq. (3.17). Like eq. (3.18) the additional 2HDM contribution from the diagrams of figures 3c and 3d is obtained as

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 section 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 Yukawa-dependent 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.   quartic scalar boson couplings. The 3-boson diagrams of figure 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. (3.25) appears multiplied with Y h l , and the coupling eq. (3.26) is multiplied with Y H l and Y A l . This allows to read off which combinations of the parameters ζ l , η, Λ 5 appear in these diagrams. With these considerations, we can rewrite a Yuk µ as a Yuk 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. (3.17) and (3.19) for H are included in a 1 0,z . The parameter dependence of each coefficient a k i,j is rather simple. a 0 0,0 and a 0 5,0 are dependent only on M H ± and the rest dependent only on M H and M H ± . In appendix A we present the explicit expression of the coefficients a k i,j as well as a non-Yuk µ , and in appendix B we show that there is no dependence on M A .
The plots in figure 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. (3.25), 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. (3.26) and appear enhanced by large t β and ζ l in eq. (3.27).
The plots in figure 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. Figure 7a shows the generic diagrams for neutral Higgs bosons while figure 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.

JHEP01(2017)007
Our result for neutral Higgs bosons is coincident with previous analysis 3 [38,44,45], and the explicit form is and for S = A

(3.32)
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. Figures 7b and 7c show the corresponding Feynman diagrams. Especially the result of figure 7c is divergent, the corresponding counterterm diagram is shown in figure 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 [56,57] 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 [38].

(3.39)
After applying the Aligned 2HDM Yukawa coupling constants in eq. (2.17) we can rewrite eq. (3.39) with ζ f , and the result reads  Table 3. Relation between signs of the aligned parameters and the functions depicted in figure 8.
The terms in the second line are proportional to ζ f η and are illustrated in the fourth plot, figure 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 section 4, ζ 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 figures 8a-8c is that they all decrease with the mass of the scalars. Figure 8d shows the contribution proportional to η which comes from diagrams involving CP-even scalar bosons. As presented in eq. (3.40) 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.

JHEP01(2017)007 4 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 η.
For the analysis we choose physical free input parameters, the masses of the different scalars (M H , M A , M H ± ), the alignment parameters (ζ l,u,d ), t β , the expansion parameter η, and Λ 5 . As presented in section 2, the last parameter can be expressed in terms of λ 1 , which is directly constrained by stability and perturbativity. Therefore, for the numerical analysis, it will be useful to replace the parameter Λ 5 with λ 1 . For ζ f we adopt the parameter range in ref. [42] 0 < |ζ u | < 1.2, 0 < |ζ d | < 50, 0 < |ζ l | < 120, (4.1) and the rest of the parameters satisfy 125 < M H < 500 GeV, M A < 500 GeV, 80 < M H ± < 500 GeV, 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 figures 9a-9c (for the three values η = 0, η = 0.1, η = −0.1). We then apply the further experimental/theoretical constraints discussed in section 2.2. The surviving sample points are depicted as red points in figures 9a-9c.
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 . Figure 9d shows the influence of the upper limit on the parameters λ i . The blue points in this plot are defined like the blue points of figure 9c, except that now 0 < λ 1 < 1. The black points of figure 9d satisfy in addition |λ i | < 1 for all i = 1 · · · 5. Like in the previous plots, the red points are the allowed points after all other constraints are applied. We find that the possible range of a B µ is only slightly smaller. Hence very large values of the quartic couplings λ i are not essential to obtain significant a B µ . We show only the plot for η = −0.1, but the discussion is valid for other values of η as well.
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 JHEP01(2017)007 Figure 9. Scatter plots showing possible values for a µ and a B µ evaluated at different η values. For plots (a)-(c), blue dots represent points in the general allowed parameter space eq. (4.2) while red dots represent the remaining ones after the constraints are applied. Green and yellow triangles are representative points discussed in the text. In plot (d), we show the influence of the maximal value allowed for the quartic couplings. Blue dots represent points in a modified allowed parameter space (the range of λ 1 in eq. (4.2) is replaced to 0 < λ 1 < 1) while black dots represent the remaining ones after imposing that all quartic couplings, λ i (i = 1 · · · 5), are less than 1. Red dots are the surviving sample points after the further electroweak/theoretical constraints are applied. mainly a B µ . For comparison with the previous analysis [36], 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. [36], 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 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. Figure 10 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 [36], 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. Figures 10a-10b 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) × 10 −10 , and the contribution can either be negative or positive. The behavior can be understood by analyzing the formula for Λ 5 , eq. (2.13), the general formula for a B µ , eq. (3.23), and the values of the different coefficients, figures 4-6.
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 figure 10c and the peak in figure 10b.
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 figure 10a and the independence of λ 1 in figure 10d.
Regarding the η-dependence, the dominant terms depending on η are a EW add. In order to compare our analysis of figure 10 with the scatter plots of figure 9, we show three representative points in figure 9. The first, in green, is the representative point just discussed for the large t β region (t β = 100, λ 1 arbitrary). The other two, yellow and dark yellow, are related to the small t β region 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 4 It should be noticed that any other point considered in [36] 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 [39,42] also considered τ -decay as a parameter constraint in the 2HDM, which disfavors a significant part of the preferred parameter space. Nevertheless, reference [42] found the general parameter region represented by eq. (4.3) to be viable. 5 The analysis is unaltered for other choices to λ1, only the absolute value of the bosonic contribution is modified. values for a B µ (2 · · · 3) × 10 −10 . This behavior is explained by observing figures 10a-10d 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 small t β region, 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 figures 10b-10c 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 η, figure 9c, has more allowed points with values for |a B µ | of order (2 · · · 4) × 10 −10 if compared with the other cases. 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 figure 11a, 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 and decrease in modulus with M H . Therefore, the net result will be an increase in a µ as observed in figure 11b.
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. (3.23) and eq. (3.27). 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 section 3.3.
We also confirmed the previous result of the fermionic contribution. Particularly, we presented its analytic form without one-dimensional integral relations in section 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. (3.24).
In the numerical evaluation we confirmed that the fermionic 2HDM contribution can be of the order of the deviation eq. (1.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.

A Analytic results
Here we provide the full analytic result of the complete renormalized bosonic two-loop contributions a B µ , in the decomposition of eq. (3.23). We begin with required loop function (defined first in ref. [86]): The coefficient a non-Yuk µ of the contribution without Yukawa couplings is given by The abbreviations appearing in a non-Yuk JHEP01(2017)007 , 14) The following coefficients depend only on known SM parameters; we provide therefore also approximate numerical values: