Flavor constraints on the Two Higgs Doublet Models of $Z_2$ symmetric and aligned types

We give a comprehensive study from flavor observables of pion, kaon, D_(s), and B_(s) mesons for limiting the Two Higgs Doublet Models (2HDMs) with natural flavor conservation, namely, Z_2 symmetric and aligned type of models. With use of the updated studies and analyses of B ->tau nu, D ->mu nu, D_s ->tau nu, D_s ->mu nu, K ->mu nu, Pi ->mu nu, B^0_s ->mu^+ mu^-, B^0_d ->mu^+ mu^-, tau ->K nu, tau ->Pi nu, B ->X_s gamma, K-K bar mixing, B^0_d-B^0_d bar mixing, and B^0_s-B^0_s bar mixing, we obtain constraints on the parameters in the 2HDMs. To calculate the constraints, we pay attention to a determination of CKM matrix elements and re-fit them to experimental data so that new contributions from additional Higgs bosons do not affect the determination. In addition, we discuss excesses of observables in the muon anomalous magnetic moment and the semi-tauonic B meson decays in the context of the 2HDM.


Introduction
Many experiments have investigated the validity of the standard model (SM). It has turned out that the SM can mostly accommodate the present low-energy experimental data. One of the most impressive features of the SM is the flavor structure. It is described by the Cabibbo-Kobayashi-Maskawa (CKM) matrix [1,2] in the Yukawa sector, which in turn requires only one Higgs doublet. Charged currents and CP -violations in the quark sector are controlled by the CKM matrix elements, whereas flavor changing neutral currents (FCNCs) are naturally suppressed due to the Glashow-Iliopoulos-Maiani (GIM) mechanism [3]. Extensions of the SM, especially in the Yukawa sector, easily and/or drastically change such a structure. Thus, flavor observables are good tools to test new physics involving such an extension. Two-Higgs-Doublet Models (2HDMs) are minimal extensions of the SM with a limited number of parameters [4]. In the general 2HDM, there exist tree-level FCNCs in Higgs boson interactions with fermions. A condition of natural flavor conservation [5] is usually imposed on the Yukawa sector to avoid such dangerous FCNCs. In principle, there are two ways to realize natural flavor conservation in the 2HDM. The first one is to introduce a Z 2 symmetry in the Lagrangian so that each fermion doublet couples only to one of the two Higgs doublets. Under this condition, there are four distinct types of Z 2 symmetric model, usually referred to as type I, II, X ("lepton specific"), and Y ("flipped") models (see, e.g., Ref. [6]). A more general method is to impose an alignment condition on the flavor space of the Yukawa matrices [7], which we call as an aligned model in this paper. From the viewpoint of low-energy phenomenology, the differences between these five models appear only in the parametrization of the Yukawa interaction terms that are sources of flavor transitions. The 2HDMs have been under detailed investigation for a long time, in particular the type II model, which corresponds to the Yukawa interaction in the minimal supersymmetric SM. In Ref. [8], a constraint on the 2HDM of type II from flavor observables have been calculated as well as the global fit of the CKM matrix elements, taking the charged Higgs effect into account. In Ref. [9], a comprehensive study for flavor phenomenologies in the aligned model has been done. A study for the Z 2 symmetric models is also given in Ref. [10]. As for studies focusing on individual processes and collider phenomenology in the 2HDMs, there are numerous papers these days (see Ref. [11] and its references and citations for a review).
In recent years, theoretical calculations of higher order corrections on several flavor observables have been developed. The next-to-leading order (NLO) electroweak (EW) and next-to-next-to-leading order (NNLO) QCD corrections on B 0 s,d → + − in the SM have been evaluated in Refs. [12][13][14]. The contribution of the extra Higgs bosons in the 2HDM is also known [15,16]. In Refs. [17,18], the updated NNLO QCD prediction of the branching ratio for B → X s γ in the SM, considering all the available non-perturbative effects, has been reported, and the NNLO contribution in the 2HDM was also obtained in Ref. [19]. The complete one-loop calculation of the EW correction on the B 0 s -B 0 s mixing in the 2HDM is given in Ref. [20]. Moreover, lattice studies on non-perturbative QCD quantities such as the meson decay constants and the bag parameters of neutral meson mixings have also been updated and collected recently [21].
Taking all the recent updates both for the theoretical calculations and the experimental results along with the other significant processes, we perform a comprehensive study of the type I, II, X, Y, and aligned models and obtain the current status of the flavor constraints on the parameters in these models. In our study, we consider: the leptonic meson decays B → τ ν, D → µν, D s → τ ν, D s → µν, K → µν, π → µν, B 0 s → µ + µ − , and B 0 d → µ + µ − ; the hadronic tau lepton decays τ → Kν and τ → πν; the radiative B meson decay B → X s γ; and the neutral meson mixings ∆M s , ∆M d , and | K |. In addition, we discuss several observables in which deviations between SM predictions and experimental results have been reported, such as the semi-tauonic B meson decays R(D ( * ) ) and the muon anomalous magnetic moment ∆a µ , in the context of the 2HDMs. We also summarize formulae for the observables of these processes in the 2HDM and utilize them to calculate constraints on the parameters. To obtain the constraints, we underline all the uncertainties taken into account in our evaluation and, in particular, pay close attention to the determination of the CKM matrix elements since this is one of the dominant uncertainties. In this paper we consider a CP conserved Higgs potential, which means that the CP-odd Higgs boson is not mixed with the CP-even Higgs bosons.
Our paper is organized as follows. In Sec. 2, we review the Yukawa sector and the parametrization of the 2HDM with the hypothesis of natural flavor conservation. The flavor observables used in our analysis are summarized in Sec. 3, and we show useful formulae for them. In Sec. 4 we obtain values of the CKM matrix elements by re-fitting CKM parameters to experimental data, so as to avoid effects from extra Higgs bosons in the 2HDM. In Sec. 5 we show the inputs and experimental data, and then we obtain the constraints on the 2HDMs. We also discuss the anomalies in R(D ( * ) ) and ∆a µ . A summary is given in Sec. 6.

Yukawa sector in 2HDM
When we consider two Higgs doublet fields Φ 1 and Φ 2 in a model with the SM fermion field contents, the most general Yukawa Lagrangian is given by where Y a f are flavor mixing complex matrices andΦ a ≡ iσ 2 Φ * a . The vacuum expectation values (VEVs) are defined as Φ a = (0 v a ) T . In general, this term immediately induces FCNCs via neutral Higgs bosons even at the tree level. To see it clearly, we can change the basis of the Higgs fields Φ i into Ψ i so that where tan β = v 2 /v 1 and H ± (A) is a charged (CP odd) Higgs boson. The neutral Higgs fields are indicated as h 1 and h 2 , which are not yet in the mass eigen basis. In this basis, the SM Higgs VEV (v = v 2 1 + v 2 2 ) and the NG bosons (G ± , G 0 ) are contained in Ψ 1 . Thus, we can rewrite the Yukawa Lagrangian in (2.1) as whereŶ ij f are the SM Yukawa matrices which derive fermion mass matrices, and ρ ij f are new couplings which are in general not diagonalized in the mass eigen basis. Note that ρ ij f do not contribute to the fermion masses. The forms ofŶ f and ρ f are described in terms of the original matrices Y a f aŝ We can see that the off-diagonal elements of ρ f cause FCNCs in the neutral Higgs sector at the tree level. So, it is required to impose a natural condition so that such FCNCs are suppressed, namely natural flavor conservation [5], or directly constrain parameters inducing FCNCs by experiments. It is known that there are two ways to archive the flavor conservation in the neutral current.

Z 2 symmetry
It is well-known that the FCNC can be prohibited by imposing a Z 2 symmetry on the fields in the Yukawa sector [4]. This is realized so that two Higgs doublets Φ 1 and Φ 2 are assigned to be Z 2 -even and -odd respectively such as Φ 1 → +Φ 1 and Φ 2 → −Φ 2 . Due to this assignment, each field f (= u, d, ) cannot have one of two original Yukawa matrices Y a f , which immediately leads to the relation ρ ij f ∝Ŷ ij f and thus the FCNC term does not appear in the Lagrangian. The protection against FCNCs is valid at any scale in the Z 2 symmetric models [15], as can be seen later in the formula for B 0 q → + − . This procedure leads to Table 2. The relation between the Z 2 models and the aligned model.
four distinct types of Yukawa interactions. In the mass eigen basis, it is summarized as where V ud is the CKM matrix element, P R L = (1 ± γ 5 )/2, h and H are the neutral CP even Higgs bosons in their mass eigen basis obtained by (h 1 h 2 ) T = R(α − β)(H h) T , and ξ φ f are scaling factors of the Yukawa couplings, dependent on the type of the model. The explicit expression for ξ φ f is listed in Table 1. We see that the Yukawa sector is controlled by tan β and sin(α − β) in this model.

Alignment
Another way to naturally forbid the tree level FCNC is worked out by taking the two Yukawa matrices to be aligned such as where ζ f for f = u, d, and are family universal scaling factors 1 . The parameters ζ f can be complex values, in principle. If ζ u,d are complex values, determinations of the CP phase in the CKM matrix are affected. In this paper, we take ζ f to be real. In this case, the scaling factors of the Yukawa couplings in (2.6) are written as (2.10) We can see from Table 2 that the Z 2 symmetric types can be considered as the particular cases of the aligned model. In this model, the alignment condition as in (2.7) is only 1 We can also define the alignment condition such as Y 2 f =ζ f Y 1 f . Under this condition, ρ f is also proportional toŶ f . The relation betweenζ f and ζ f is written by guaranteed at the scale where the model is defined. Namely, non-zero contributions from running effects can be generated at a low energy scale, see e.g., Ref. [15]. We see such effects in B 0 q → + − later. We have practically sin(α − β) and ζ f as the model parameters in the Yukawa sector and tan β is irrelevant to this sector. A typical difference of this model from the Z 2 symmetric models is that ζ f is nothing to do with the fermion mass and the VEVs. Hence it can be arbitrary.

Type III
Indeed, we can consider the case that the FCNC interactions appear at the tree level. This class of model is called as a type III model and directly obtained from (2.4) without taking any condition. In this model, there are a lot of couplings which can induce the FCNC transitions and thus, they are severely constrained by experiments. For an overview of this model and flavor constraints, see Ref. [22]. Recent studies for the top quark FCNC in the type III are given in Refs. [23][24][25][26][27][28][29][30]. As for the lepton flavor violating decay of the Higgs, there are several investigations in Refs. [31][32][33][34]. The decay can be related to a muon anomalous magnetic moment in this class of model [35]. Some flavor structures can be derived if a global symmetry is imposed on the fields. The BGL model has been proposed by considering such a symmetry [36]. This model prohibits FCNCs in the up-type quarks. On the other hand, FCNCs in the down-type quarks are controlled by the CKM matrix elements. This class of models has been recently analyzed in Refs. [37,38].

Flavor observables
In the 2HDM, the charged Higgs boson H ± can contribute to flavor observables via the charged current. By considering the Z 2 symmetry or the alignment condition, interesting features appear in the Yukawa interaction term as follows. The tree-level interaction term of H ± with quarks has the same CKM structure with that of W ± . Higher oder corrections due to loop diagrams through extra Higgs bosons provide additional contributions to FCNC processes in the quark sector. On the other hand, the lepton flavor violation is quite suppressed due to a tiny neutrino mass contribution.
A sensitivity of an observable to a new physics model, from the viewpoint of limiting a new physics model, depends on the situation such that how much precisely the observable is measured and how much the parameters involved in the new physics model can enhance the effect on the observable. To clarify which observables are important to constrain the 2HDM, we classify flavor observables as follows. Obviously, observables in the category (a) are significant in order to limit the 2HDM. The category (b) plays an important role when we fit the CKM matrix elements in the 2HDM. In the SM, they can be determined by a fit to experimental data for each element, or we can also perform a global fit to a specific parametrization. In the 2HDM, however, some of observables are affected by the extra Higgs bosons and thus the CKM matrix elements must be re-fitted by taking such an effect into account. In particular, the CKM matrix elements should be determined by the observables classified as (b-1) and not (b-2). The category (c) means that a discrepancy between a SM prediction and an experimental result exists in an observable. Such an observable is also important for an indirect evidence of the additional scalar bosons.
In the following subsection, we summarize theoretical formulae of flavor observables classified as (a) and (c), which are useful to constrain the 2HDM. After that in the next section, we discuss the way to obtain the CKM matrix elements by the global fit with use of observables in (b).

M
In the 2HDM, decay processes M ± → ± ν and τ ± → M ± ν occur at the tree level and their branching ratios are given by where M + is the meson which consists of (ud), f M is its decay constant, τ M is its lifetime, and V ud is the relevant CKM matrix element in the process. In the following part, we omit the notation of charge assignment unless otherwise stated. The contribution from the charged Higgs boson is encoded in C H , written by where ξ A f is defined in (2.6) and applied to each specific model. The main sources of non-negligible uncertainties for theoretical evaluations are f M and V ud as shown later.
The branching ratios of B → τ ν, D → µν, D s → τ ν, and D s → µν are notable observables in the 2HDM. The processes B → µν and D → τ ν are less sensitive to new physics since only upper limits of the branching ratios are given by experiments for now. For the K, π, and τ decays, the ratio can be significant thanks to smaller theoretical uncertainties and precise experimental data. In this case, the long distance electromagnetic corrections enter as (1+δ K/π EM ) and (1+δ K/π,τ EM ) in (3.4), respectively. These corrections are very small, but include uncertainties taken into account. Input values are shown in Sec. 5.

M 0 → + −
Pure leptonic decays of neutral meson M 0 → + − can probe the effect of the 2HDM via loop contributions. For instance, the branching ratio of B 0 q → + − can be written by where the Wilson coefficients C 10 , C P , and C S show the contributions from the effective respectively. Since the branching ratio is quite small, the effect of neutral meson mixing is not negligible. To involve such an effect, the averaged time-integrated branching ratio is defined [39][40][41] as where τ H(L) B 0 q is the life time corresponding to the heavier (lighter) eigenstate of B 0 q , and ∆Γ q is the decay width difference. This is the actual observable which can be compared with experimental data. A brief review of this observable is shown in Appendix B.1.
The coefficient C SM 10 includes the dominant SM contribution from O 10 . At the leading order, it is evaluated as C SM-LO where we define x q as From recent studies in Refs. [12][13][14] evaluating C SM 10 up to the NNLO QCD and NLO EW corrections, we obtain the following fit formula: where α s (m Z ) is the QCD running coupling at the m Z scale and M t shows the pole mass of the top quark. In our analysis, we use (3.11) as the SM prediction. The scalar and pseudo-scalar type effects in C S and C P are suppressed in the SM. The explicit one-loop order calculation for B 0 q → + − in the aligned 2HDM has been done in Ref. [15]. The similar calculation in the Z 2 symmetric 2HDMs is also given in Ref. [16]. In the 2HDM, the charged Higgs boson contributes to C 10 via Z-penguin diagrams. It is described as (3.12) The neutral Higgs bosons give contributions in C S and C P . The contributions are divided by two parts such as The first terms C c S and C c P come from box diagrams in the Unitary gauge (box, Z-penguin, and Goldstone-penguin diagrams in the Feynman gauge) and are given by where C c, SM S,P show the SM contributions that are given in Appendix B.1. The second terms C n S and C n P indicate the contribution from scalar boson exchange diagrams. They can be expressed as where c θ = cos θ, s θ = sin θ, and The analytic expressions for G i and F i are given in Appendix B.1. Note that the effects from C P and C S on the branching ratio is suppressed by the factor m 2 004 as seen in (3.7). In the Z 2 symmetric model, however, C P,S can be enhanced for large tan β, relative to C 10 . In the aligned model, we can obtain dominant constraints on ζ d and ζ from this observable. Further details for this analytic formula are found in Appendix B.1. The dominant uncertainties for this process are the decay constant f B 0 q and the product of CKM matrix elements V tb V * tq as well as the M → ν case. The formulae for D 0 → + − and K L → + − are easily given by replacing masses and CKM components as appropriate. In these two cases, however, there exist non-negligible long distance contributions and thus we must concern that effect in addition to the short distance contributions given above. Then they are less significant to constrain the 2HDMs.

Neutral meson mixing
The 2HDM gives additional contributions to neutral meson mixings and affects measurements of the CKM matrix elements from mixing observables. For the B 0 q -B 0 q mixing, the mass difference between two CP eigenstates of B 0 q andB 0 q defined as is used to determine the parameters in the CKM matrix in the case of the SM, since this is less sensitive to a long distance effect. In the 2HDM, the exchanges of the charged Higgs boson in the box diagrams contribute to ∆M q [42,43]. A complete one-loop calculation without neglecting the term proportional to x b is given in Ref. [20]. We express the analytic formula as   [20] by taking nonzero external momenta into account and thus they are different from those given in Ref. [43]. We independently obtained the same result for (B.39) and (B.40). In addition, we found the x b terms in A W H and A HH as shown in (B.37) and (B.38). We stress that they are new contributions which are not described in Ref. [20].
As for the K 0 -K 0 mixing, the mass difference ∆M K is not a good observable due to a non-negligible long distance effect. Instead, the parameter which is defined as is the measurement of the CP violation in the K 0 -K 0 system. In the 2HDM, the analytic formula can be written [43] as where A V V are the same forms with those for ∆M q , B V V are additional functions expressed in Appendix B, and η qq indicates the QCD corrections for each pair of the internal quarks. In this formula, x s = 0 is assumed and thus there is no contribution from the SRR and T RR operators 2 . Note that we use the relation ∆Γ K 2∆M K in (3.23) and the experimental data of ∆M K as the input value in our numerical analysis. For K and ∆M q , the uncertainty from the bag parameterB is taken into account as well as the decay constant and the CKM matrix element. Later we show and summarize 2 It is often stated that the matrix element for the scalar type operator can give a larger contribution than that for the SM vector operator due to a chiral enhancement. But, the loop diagram for the coefficient of the matrix element provides the suppression by xs in the 2HDM and then in total the contribution from the scalar sector is suppressed by ∼ m 2 the detail of input parameters. As for the QCD corrections in the V LL operator, we use the following values [44], The QCD correction of η Bq from the extra Higgs bosons has been obtained in Ref. [45]. In our study, we simply neglect that effect. Practical input values of the bag parameter and QCD correctionB ST Bq η ST Bq in the SRR and T RR operators are shown in Sec. 5. In the D 0 -D 0 mixing, there is no feasible observable in which a long distance effect is sufficiently suppressed. In principle, we can give a constraint in the D 0 -D 0 mixing by taking such an unknown effect as one of uncertainties. This strategy can be applied to other observables such as K , ∆M K and so on, but in this paper we do not consider such a case.

3.4B → X q γ
Inclusive radiative decays of B meson,B → X q γ, are one of the most interesting FCNC processes, and hence they have been precisely evaluated to bound on several new physics models. Since the hadron transition occurs in this process, perturbative QCD corrections are much important and non-perturbative effects must be concerned as well. According to the recent summary that collects all the available and relevant contributions presented in Refs. [17,18], the SM predicts for E 0 = 1.6 GeV, where E 0 indicates the photon cutoff energy and B(b → qγ) is the CP-and isospin-averaged branching ratio ofB → X q γ. For the evaluation of (3.25), perturbative QCD corrections up to the NNLO [46][47][48][49][50] and calculable long-distance effects, (see Ref. [18] and its references for more detail), are taken into account. The uncertainty in (3.25) comes from non-perturbative part (±5%), input parameters (±2%), and others (∼ ±4%). The product of the CKM matrix elements, defined as is contained in the input parameters. In Refs. [17,18], the latest result of the CKM fit in Ref. [51] are applied. To separate the CKM product from the observable, we employ the following expression for the SM prediction: where δr V denotes the uncertainty derived from r V . We note that δr V obtained from Ref. [51] is negligible. Differently from the other observables in the present paper, we utilize the uncertainty already evaluated by the other well-sophisticated study as shown above.
The contribution from the charged Higgs boson is provided through the effective oper- The LO correction is gained by one loop EW and QCD penguin diagrams and the NLO correction is calculated as in Refs. [52][53][54][55][56]. To incorporate these effects, we obtain the useful numerical formula in terms of the Wilson coefficients at the scale µ t = 160 GeV as follows, show the new physics contributions from O i for i = 7, 8 at the LO (NLO) part. The SM contributions are separated in advance. This numerical formula is estimated based on Refs. [57][58][59] and we confirmed that the C LO 7 term, which is the most dominant, is consistent with Ref. [17]. The explicit form of the coefficients C LO i and C NLO i is given by for i = 7, 8 and µ t = 160 GeV, with respect to the mass ratio, The loop functions G i a , C i a , and D i a are described as in Appendix B.3. The NNLO QCD correction to this process has been studied in Ref. [19]. Within the present status on the uncertainties in the experimental result and the theoretical prediction, ignoring such a correction does not much change our result. As for B(b → dγ), the uncertainty in the prediction is still large and thus not of importance in limiting the 2HDM yet.

Anomalies
In this subsection we focus on formulae of observables in which discrepancies between the SM prediction and the experimental result have been reported, categorized as (c). The current status on the discrepancies are summarized in Sec. 5.

3.5.1B → D ( * ) τν
Semi-tauonic B meson decaysB → D ( * ) τν are sensitive to the effect of the charged Higgs boson since its contribution is proportional to [60,61]. The results from the Belle, BaBar, and LHCb collaborations are nowadays available despite that it is difficult to identify the tau in these processes. Useful observables of these decays are given by , (3.34) in which we can reduce uncertainties coming from input parameters. The SM predicts precise values of R(D ( * ) ) with the help of the heavy quark effective theory to evaluate form factors [62,63]. The effect on R(D ( * ) ) in the context of the 2HDMs has been calculated as seen in Ref. [64][65][66][67][68]. Based on our previous study in Ref. [69], we give the numerical formulae of the branching ratios for the 2HDMs as follows. For the branching ratios of the light-leptonic modes, where ρ 2 1 , ρ 2 A 1 , R 1 , and R 2 are the form factor parameters fitted by energy distributions; V 1 (1) and A 1 (1) are overall normalizations of the form factors; and Γ D ( * ) i are the coefficients of polynomial expansion with respect to ρ 2 1 (ρ 2 A 1 ). The explicit forms are given in Appendix B.4. In this formula, we neglect the charged Higgs contribution to B(B → D ν) and B(B → D * ν) since it is suppressed by the factor, m q m /m 2 H + . In the muonic mode, the contribution can be potentially a few % and it can affect the determination of |V cb |. We will discuss it in the next section. In the 2HDM, the numerical formulae for B(B → D ( * ) τν) are written as where the charged Higgs contributions included in C S 1 and C S 2 are written as As seen in the formulae, the overall normalizations V 1 (1), A 1 (1) and the CKM matrix element |V cb | are irrelevant to R(D ( * ) ). Later we show the fitted result for the input parameters.

Muon anomalous magnetic moment
The muon anomalous magnetic moment a µ provides a sensitive test of quantum loop effects in the electroweak sector. The SM contributions are evaluated as in Refs. [70][71][72] including several higher oder corrections [73][74][75][76]. Recent studies for higher order corrections are also obtained in Refs. [77][78][79][80]. A discrepancy between the experimental result reported in Ref. [81] and the SM prediction, a exp.
µ − a SM µ , can be compared with a new physics contribution.
In the 2HDM, the complete one-loop diagrams and the Barr-Zee type two-loop diagrams can be significant. The analytic formula for the one-loop diagrams is given [82][83][84] by where y f φ is defined in (3.33) and the loop functions F φ are calculated as The result for the Barr-Zee type two-loop diagrams is given [73,[85][86][87][88] where the index f represents the fermion in the loop, Q f and N c f are the electric charge and color degrees of freedom of f . The functions G φ are obtained by In the SM prediction, the contribution from the SM Higgs boson is already taken and thus we must care about this part when considering the 2HDM. Substituting the corresponding contribution, the 2HDM contribution which can be compared to a exp.
µ − a SM µ is represented as

Determination of CKM in the 2HDM
As stressed in the previous section, it is necessary to concern the effect of the extra Higgs bosons when we determine the CKM matrix elements by fitting to experimental data, in the 2HDM. This is expected to be more crucial for the future flavor experiments at the SuperKEKB/Belle II [90]. The global fit of the CKM matrix elements, together with the parameters of the 2HDM, to all the relevant experimental data is one of the approaches for the analysis [8]. In this paper, we employ a more visualized approach as follows. For the re-fit of the CKM matrix elements, we use the Wolfenstein parametrization which is defined as where we neglect O(λ 4 ) ∼ O(0.001). Then, we obtain fitted values of λ, A, ρ, and η by using observables in which contributions from the extra Higgs bosons are negligible. As for r V needed in the evaluation forB → X s γ, we take r V 1 − λ 2 (1 − 2ρ).

λ and A
The most precise value of the Cabibbo mixing parameter λ is provided from the determination of |V ud | by the super allowed (0 + → 0 + ) nuclear beta decays. The experimental result |V ud | = 0.97425 (22) [91] implies λ = 0.2269 ± 0.0010. In the SM, λ is also determined from leptonic K decays such as K → (π) ν and τ → Kν for = e, µ. Among them, K → eν is usable to determine λ in the 2HDM, since the effect of the extra Higgs bosons is safely negligible and its experimental data is available separately from the muonic mode. The experimental result B(K → eν) = (1.581 ± 0.008) × 10 −5 [92] is translated into λ = 0.2221 ± 0.0014, where the decay constant of K we used is listed in Table 3. Therefore the combined result is given as and we use this value for the following analysis in this paper.
The parameter A is included in V ub , V cb , V td , and V ts , and usually obtained from the determination of |V cb |. It is, however, known that the values of |V cb | obtained from inclusive (B → X c ν) and exclusive (B → D ( * ) ν) decay modes are not in good agreement [92,93]. In the 2HDM, although the charged Higgs boson affects the muonic modes ( = µ), it is hard to compensate this discrepancy. In the present paper, we simply obtain a combined value of |V cb | considering the charged Higgs effect. For the determination fromB → X c ν, a combined fit to moments of several variables (a hadronic-mass, a lepton-energy spectrum, and a photon-energy spectrum) are required. Calculating the charged Higgs effect on its distribution is beyond the scope of this paper. Instead, we roughly estimate such an effect by using the expression defined as where δ H indicates the contribution of the charged Higgs boson to the muonic decay mode and |V cb | obs. is the experimental result of the fit by assuming the SM. The coefficient C stems from the difference of the effective operator between the SM and the charged Higgs contributions. Considering the quark level process, which is the leading order contribution involved inB → X c ν, we obtain C 0.05 for = µ. Then we find that the correction from C δ H is less than 1% for m H + > 150 GeV and ξ A d = ξ A = tan β < 100 in the type II and aligned models as shown in Fig 1. As for the other types, it is completely negligible. The similar estimation can be done forB → D ( * ) ν. In these exclusive processes, we obtain C 0.15 (0.017) with use of the formula forB → D ( * ) τν replacing m τ with m in Ref. [69]. In Fig 1, the correction to the measurement of |V cb | is shown in the type II model. We can see that the charged Higgs effects inB → Dµν are more important than that in the inclusive process but, in any case, for m H + 300 GeV they are not sizable. As for the combined experimental value of |V cb |, we refer to the latest determination by the CKMfitter [51] (not the result from the global fit in this reference). To conclude, we take A = 0.808 ± 0.017 , (4.5) in the case of the type I, X, and Y models. We add an additional uncertainty to (4.5) in accordance with δ H in the case of the type II and aligned models.

ρ and η
The CP phase in the CKM matrix is given by ρ and η. For the actual observables,ρ + iη = −(V ud V * ub )/(V cd V * cb ) is defined and measured by experiments, where it is related as In the SM,ρ andη are fitted by several variables such as K , ∆M d , ∆M s , |V ub |, and the angles of unitarity triangle, defined as In the 2HDM, we note that measuring these angles is not affected by the extra Higgs bosons as long as ξ A f is real, whereas the others are potentially harmed. Thus we use only the unitary triangle to determineρ andη. The latest world averages of the angles are provided in Ref. [51] and then related as and we show the (ρ,η) plot for the fit in Fig. 2.
We note that the SM global fit reported in Ref. [51] Figure 2. Allowed regions on the (ρ,η) plane obtained from the measurements of α = φ 2 (yellow), β = φ 1 (blue), and γ = φ 3 (green). The red region indicates the combined result from these three measurements. The right panel is zoomed version of the left panel.

Constraint
Here we give current constraints on the Z 2 symmetric and aligned types of the 2HDM with GeV , ∆M s , ∆M d , and | K |. The way to evaluate uncertainties and exclusion confidence levels (CLs) for the above observables is shown in Appendix A. To begin with, we exhibit input required to evaluate the observables and the experimental results. After that, we obtain the constraints and comment on them. We also discuss the anomalies of R(D ( * ) ) and a µ in the context of the 2HDMs with the natural flavor conservation.

Input and experimental data
We apply our fit result for the Wolfenstein parametrization to the CKM matrix elements. Obtained from the previous section, we can express the result as The lattice studies for the meson decay constants and the bag parameters are summarized in Ref. [21], and the recent updates for f D , f Ds , and f K /f π are available in Ref. [94]. The EM corrections of B(K → µν)/B(π → µν) and B(τ → Kν)/B(τ → πν) are given in Refs. [95,96]. The values are listed in Table 3. As for the parameterB ST Bq η ST Bq in (3.19), the scale dependent expression defined aŝ  Table 3. Lattice results of the meson decay constants, the bag parameters and the electromagnetic correction evaluated in Refs. [21,[94][95][96], and input values of the initial conditions for the evaluation of the running quark masses [92].
are only evaluated. The bag parameters at the µ b = m b scale are given as [97][98][99] f Bs B (s) assuming ∆Γ d /Γ d 0 [102], where uncertainties less than 1% are neglected for these parameters. The other numerical input for our numerical analysis are shown in the Appendix A. In addition, we summarize the experimental data for the relevant observables in Table 4, along with the SM contributions which we evaluated with use of the input values shown above.

Setup of model parameters
Here, we summarize setup for the parameters of the 2HDMs in our numerical analysis. We assume the same masses for the extra Higgs bosons, m H = m A = m H + . This is favored by the truth that this relation satisfies a theoretical bound from perturbativity [88,104], and it is also allowed by the EW precision tests [4]. In this case, constraints on m A given by the ATLAS and CMS collaborations [105,106] Table 4. Experimental results of the observables combined by the PDG and/or HFAG collaborations in Refs. [92,93]. As for B(B 0 q → µ + µ − ), the combined results from the LHCb and CMS collaborations are shown as in Ref. [103].
For the mixing angle of h and H, we take the SM-like limit sin(β − α) = 1 in which the heavier CP-even Higgs boson H can not decay into W + W − and ZZ. This is justified by current Higgs boson searches at the Large Hadron Collider (LHC). The current combined fit of sin(β − α) to the LHC results has been studied in Refs. [104,107,108].
The case that sin(β − α) is close to, but not exactly, one is interesting for collider searches. From the viewpoint of flavor physics, B(B 0 q → µ + µ − ) can be affected as varying sin(β − α). But, the small difference of sin(β − α) from one changes only a few % of B(B 0 q → µ + µ − ), much smaller than the current experimental and theoretical uncertainties. For example, one finds 1.5% reduction of B(B 0 s → µ + µ − ) in the type II model for sin(β − α) = 0.9, tan β = 30, and m H = m A = m H + = 500 GeV from the sin(β − α) = 1 case. Changing sin(β − α) = 1 and m H = m A = m H + also affect ∆a 2HDM µ . Later, we will loosen these assumption and see the effect.

Constraint on the Z 2 symmetric models
In Fig. 3, we show constraints on the plane of (m H + , tan β) at 95% CL from the individual observables in the type I, II, X, and Y models. These constraints are given by evaluating χ 2 for each observable. Comments on the results for each model are as follows.

Type I:
The region tan β 1 is strongly constrained by B(B 0 s → µ + µ − ) and ∆M s in the type I model. The other observables also provide the constraints for small tan β in this type. We can see that the extra Higgs boson mass is not constrained by the flavor observables in the large tan β range.
in the lower panels. The black line contour in the type II and Y is the boundary of 95% CL exclusion fromB → X s γ. The dashed horizontal lines are ones for tan β = 1 and 0.057, corresponding to the top Yukawa coupling to be 1 and 4π, respectively. The gray region is the minimal exclusion from LEP searches [109]. The exclusion from τ → µνν is also shown in the type X [110].

Type II:
In the type II model, the dominant constraint comes from B(B → τ ν) and B(B 0 q → µ + µ − ) for large tan β. The branching ratio B(b → sγ) Eγ >1.6 GeV gives the lower limit on the mass. Our result shows that m H + < 493 GeV is ruled out at 95% CL and close to what was reported in Ref. [17]. Moreover, m H + < 408 GeV is excluded at 99% CL. The loop induced processes such as the neutral meson mixings exclude the region for small tan β as well as the type I model.

Type X:
As for the type X model, the processes M → ν and τ → M ν provide no constraint on the (m H + , tan β) plane from the current data, whereas the loop induced processes exclude the range for small tan β as well as for the type I case. Indeed, as the tan β enhancement can be seen only in the lepton sector, we can put a constraint for large tan β region from the measurement of the Fermi constant G F from τ → µνν [110]. In the figure, we show the result from τ → µνν, where we obtain the theoretical formula (at the one loop level) and the experimental data based on Ref. [110]. A similar bound is obtained in the type II model but we have checked that it is smaller than the one from B(B → τ ν) and B(D s → τ ν). To conclude, the type X model does not have a significant exclusion for the mass.

Type Y:
The type Y model is constrained by B(B 0 s → µ + µ − ) and ∆M s for small tan β as is the same with the other models. The lower mass limit is obtained by B(b → sγ) Eγ >1.6 GeV , similarly to the type II model, as m H + < 493 GeV (408 GeV) at 95% CL (99% CL), because of the same couplings, ξ A u = 1/ tan β and ξ A d = tan β.

Constraint on the aligned model
Next, we see constraints on the aligned model. The tree processes M → ν are insensitive to the aligned model unless the products of ζ u ζ and/or ζ d ζ are very large. However, the large value of ζ u is constrained by the neutral meson mixings as shown in Fig. 4. The results do not much depend on ζ d since the term including ζ d is proportional to x b in ∆M s and ∆M d and there is no dependence in | K |. We can see that the small value of ζ u is only allowed, e.g., |ζ u | < 1.5 for m H + = 1000 GeV and |ζ u | < 3.5 for m H + = 4000 GeV, as has been pointed out in Ref. [9]. On the other hand, ζ d can be limited by B(B 0 q → µ + µ − ) and B(b → sγ) Eγ >1.6 GeV . In Fig. 5, we show the constraints on (m H + , ζ d ) from these observables as varying ζ u and ζ . The upper and lower figures are the results for ζ = ζ d and ζ = 0, respectively. Constraints in the case of the negative value of ζ u are obtained by replacing ζ d to −ζ d in the vertical axis of these plots. The parameter ζ is irrelevant for b → sγ. For ζ u = 0 and ζ = 0, the constraint from B(B → τ ν) becomes dominant but is insensitive for large mass. There is (trivially) no significant constraint for ζ u = ζ = 0. We can see that, for ζ u = 0, the combination of the constraints from B(B 0 q → µ + µ − ) and B(b → sγ) Eγ >1.6 GeV provides the lower mass limit. For example, we obtain the exclusions such as m H + < 1500 GeV for |ζ u | = 2 and m H + < 3700 GeV for |ζ u | = 4 at 95% CL. This is the updated result of Refs. [9,15,111].

Analysis of anomalies
In this subsection, we study the 2HDM effect on a µ and R(D ( * ) ), in which the deviations between the SM predictions and the experimental results have been reported. The muon anomalous magnetic moment has been measured by the Muon G-2 collaboration as in Ref. [81]. Since this is a high precision test of the EW corrections, higher oder contributions in the SM are important and have been evaluated. Discrepancies between the experimental result reported in Ref. [81] and the SM predictions in Refs. [70][71][72] are presented as (282 ± 91) × 10 −11 from Ref. [70] (287 ± 85) × 10 −11 from Ref. [71] (261 ± 80) × 10 −11 from Ref. [72] .

(5.7)
Even though there has been only one experimental measurement up to now, the deviation between the SM prediction and the experimental value is around 3σ as shown in (5.7). As a reference value, we take a exp.
Excess of the observables R(D ( * ) ) in the semi-tauonic B meson decays has been reported by the BaBar, Belle, and LHCb collaborations in Refs. [112][113][114][115]. The latest combined result suggests that the deviations from the SM predictions are described as  Figure 6. Allowed region plots on the (m H + , tan β) plane from a µ , R(D), and R(D * ) in the type II model. The darker (lighter) blue, orange, and green regions are the results for a µ , R(D), and R(D * ) within 1σ (2σ). The excluded regions from B 0 s → µ + µ − and B → τ ν at 95% CL are shaded with red and yellow dashed boundaries, respectively.
where the discrepancy reaches around 4σ taking experimental correlations into account. We note that asB → D ( * ) τν occur at the tree level in the SM, these deviations have an impact on limiting new physics.
In the 2HDM these three observables are affected by the Yukawa interactions of the extra Higgs bosons in (2.6). The formulae for a µ and R(D ( * ) ) are shown in Sec. 3. As for the input parameters in R(D ( * ) ), we use [93] ρ 2 1 = 1.186 ± 0.054, ρ 2 A 1 = 1.207 ± 0.026, R 1 = 1.403 ± 0.033, and R 2 = 0.854 ± 0.020. One can easily see that the type I and Y models cannot explain the present anomalies of a µ and R(D ( * ) ) at all, since no large contributions to these processes are available. The type II model can explain these anomalies individually, however, it is inconsistent with each other and also contradictory to the other constraints obtained in Sec. 5.3. In Fig. 6, we exhibit the allowed regions from a µ and R(D ( * ) ) along with the excluded regions from B(B 0 s → µ + µ − ) and B(b → sγ) Eγ >1.6 GeV as indicated in the legend. We can see that the small value of m A is required to explain the anomaly in a µ . But, in any case for the relation among the extra Higgs boson masses, the three anomalies cannot be explained at the same time, and are not consistent with the present constraints.
The type X is often discussed as one of the solutions for the anomaly in a µ . The recent studies for the type X model aiming at this anomaly are given in Refs. [88,110,[116][117][118][119].
In the upper panels of Fig. 7, we review the allowed region plot for a µ in the type X model. With small m A , it can explain the a µ anomaly. According to the study in Ref. [110], however, the constraint from τ → µνν has turned out to be severe. In the figure, we also show the 95% CL exclusion dashed line from τ → µνν. As can be seen, the explanation of the a µ anomaly at the 1σ level is not possible, but that at the 2σ level is accessible in the type X model. The result for the case of sin(β − α) = 0.9 is also shown with the black lines and one find it does not change the above conclusion. Remind that the constraints from the meson observables are negligible. Then, in any case, this model cannot accommodate In the aligned model, the parameters ζ u , ζ d , and ζ are independent and thus there is a larger parameter space than that in the Z 2 symmetric models. Nevertheless, ζ u is severely limited by ∆M s and then the constraints on ζ d and ζ come from B(B 0 s → µ + µ − ) and B(b → sγ) Eγ >1.6 GeV , which are correlated to the value of ζ u . They are less bounded if ζ u = 0 as we have seen in Figs. 4 and 5. To see how the aligned model affect a µ and R(D ( * ) ), we simply take ζ u = 0 for clarity. The possible significant constraint comes from B(B → τ ν) in this model as we can see in the leftmost panels of Fig. 5. We surveyed several parameter set and found that the anomalies in R(D) and a µ can be explained simultaneously in a specific region, such as ζ u = 0, ζ d = 5, ζ ∼ −70, m H = m H + ∼ 200 GeV, and m A = 20 GeV. The results are plotted in the lower panels of Fig. 7. We can also see that it is hard to accommodate the excess in R(D * ) without any contradiction to the other constraints. This specific point is not excluded by the other constraints yet, but close to the excluded region from B(B → τ ν). We note that the sizable value of ζ u ζ are required to explain both the excesses in R(D) and R(D * ) simultaneously as studied in Ref. [69].

Summary
We have given the comprehensive study from the observables of π, K, D (s) , and B (s) for limiting the 2HDMs with the hypothesis of natural flavor conservation, namely, the Z 2 symmetric and aligned models. Then we have obtained the significant constraints on the masses and couplings of the extra Higgs bosons, and shown the possibilities to accommodate the anomalies in the muon anomalous magnetic moment and the tauonic B meson decays.
We have considered the following flavor observables; GeV , ∆M s , ∆M d , and | K |, and collected the formulae of them, in some of which we have taken the updated calculations into account. In addition, we have also obtained the updated formula of ∆M q .
We have re-fitted the CKM matrix elements to the experimental data to which the extra Higgs bosons do not give large contributions, and evaluated the effect on the determination of |V cb |. As a result we have seen no significant difference between the re-fitted values and the global SM fit values. With the use of the re-fitted CKM matrix elements and the latest combined experimental results summarized by the PDG and HFAG collaborations, we have evaluated the excluded regions on the model parameters of the 2HDMs, with carefully considering the uncertainties from the input parameters. To obtain the results, we have assumed the same masses among the extra Higgs bosons and the SM-like limit sin(β − α) = 1, favored by the theoretical bounds, the EW precision tests, and the collider searches.
As a consequence of our work, in the Z 2 symmetric models, it has been found that B 0 s → µ + µ − plays a significant role to constraint tan β as well as the B 0 s -B 0 s mixing. The charged Higgs boson mass is constrained by the processB → X s γ in the type II and Y models. The updated theoretical evaluation and the experimental result of B(b → sγ) Eγ >1.6 GeV suggest that m H + < 493 (408) GeV is excluded at 95% (99%) CL in these two models. There is no severe constraint on the mass in the type I and X models.
In the aligned model, there are three free parameters in the Yukawa interaction term of the charged Higgs, ζ f for f = u, d, . The neutral meson mixings constrain ζ u and we have obtained severe bound as |ζ u | < 1.5 for m H + = 1000 GeV and |ζ u | < 3.5 for m H + = 4000 GeV, which are mostly independent on the other couplings, ζ d and ζ . With a non-zero value of ζ u , the parameters ζ d and ζ are limited by B(b → sγ) Eγ >1.6 GeV and B(B 0 s → µ + µ − ). For example, |ζ d | 5 is excluded for m H + = 1000 GeV and ζ u = 1. We have also shown that the combination of these two observables gives the lower mass limit as m H + < 1500 GeV for |ζ u | = 2 and m H + < 3700 GeV for |ζ u | = 4.
In addition, we have summarized the current status of the anomalies in the muon anomalous magnetic moment a µ and the tauonic B meson decays R(D ( * ) ), in the context of the 2HDMs. We have shown that the type II model can explain each anomaly, however, the allowed regions are not only inconsistent with each other but also contradictory to the other constraints. The type X model is often considered as one of the good candidates to accommodate the excess of a µ . We have reconfirmed that this model can solve the excess of a µ individually, but it is not consistent with the constraint from τ → µνν at the 1σ level. Note that this model cannot resolve the excesses in R(D ( * ) ), in any case. We also surveyed the possibility to explain the anomalies in the aligned model. We have pointed out that these three anomalies cannot be explained simultaneously, whereas the excesses of a µ and R(D) can be explained for ζ u = 0, ζ d = 5, ζ ∼ −70, m H = m H + ∼ 200 GeV, and m A = 20 GeV. This parameter set is allowed by the other flavor constraints yet, but close to the excluded region from B(B → τ ν).
We have not considered semi-leptonic meson decays such asB → πτν, B → (K ( * ) , φ)µ + µ − , and others to obtain the constraints. Although form factors in B → π, K ( * ) transitions still include large uncertainties in fit parameters, these decays can provide constraints on new physics and will become more significant at the future experiments. The recent studies for B → πτν are given in Ref. [120] and for B → (K ( * ) , φ)µ + µ − in Refs. [121,122]. Exclusive radiative B meson decays are also important, see e.g., Ref. [123].
The bounds obtained in this work are expected to be the last updated status before starting the future flavor experiments such as the SuperKEKB/Belle II [90] and the LHCb run II [124]. Future searches at the Belle II and the LHCb run II will improve sensitivities to the 2HDMs and may reveal the source of the excesses in the semi-tauonic B meson decays. Future muon g − 2 searches at the J-PARC [125] and the Fermilab [126] will also change the present situation on the anomaly. The requirement for explaining the excess of a µ implies the mass of the CP odd Higgs boson should be small. Therefore, collider signatures from h → AA → 4τ, 4b can be important.

A Evaluation of the uncertainties and input parameters
Here, we explain the way to evaluate the uncertainty in the observable coming from the one in the input parameters. Suppose the observable is expressed as F (y; {x i }), where x i is an input parameter measured (or calculated) as x i = x 0 i ± δx i , {x i } shows a set of parameters for i = 1, 2, · · · , k, and y indicates model parameters. We define the uncertainty of F (y; {x i }) as To obtain excluded and allowed regions of a parameter space y, we evaluate the χ 2 function. In our analysis it is defined as where the experimental result is shown as F exp. ± δF exp. . The parameters taken as {x i } in our analysis are listed in Table 3, and (5.3)-(5.5). The other input values used in our numerical analysis are shown in Table 5.

B Analytic formulae for flavor observables
Here we give the functions of analytic formulae for the flavor observables, which are used to obtain the constraints in this paper.
The analytic formula for the averaged time-integrated branching ratio is given in terms of the Wilson coefficients C 10 , C P , and C S . The SM contributions in (3.14) and (3.15) are written as The functions G i and F i written in (3.16) and (3.17) are described as (B.10) (B.14) (B.17) The SM Higgs contributions for C n S,P are included in (3.16) and (3.17), which can be extracted as a SM-like limit. Taking cos(β − α) → 0, sin(β − α) → 1, ξ A f → 0, and m φ → ∞ for φ = H + , H, A, the SM-like limit is obtained as

B.1.2 Assumptions
In the formulae of C n S,P , we have ignored FCNC contributions induced by a running effect of the Yukawa interaction term at the low energy scale. In Ref. [15] such contributions in C n S,P are estimated as in the SM-like limit, where C R (µ t ) shows the renormalized coupling of FCNC term in the Yukawa Lagrangian ((2.15) in Ref. [15]). We can see from (B.20) and (B.21) that R S,P = 0 in the Z 2 symmetric models. This is because that the Z 2 symmetry can protect the alignment condition at any scale. On the other hand, in the aligned model the condition is guaranteed only at the scale where the model is set and thus the non-zero contribution can appear at the low energy scale. In this paper, we simply ignore this effect in all types of 2HDM. We also neglect contributions proportional to light quark mass m q and Higgs self couplings λ 3,7 [15]. As for the Higgs self couplings, we have confirmed that the effect is negligible.

B.1.3 Averaged time-integrated branching ratio
The averaged time-integrated branching ratio B(B 0 q → + − ) can be understood as follows. The "untagged" decay rate for P → f is defined and described as with Γ P = (Γ L P + Γ H P )/2 = 1/τ P , ∆Γ P = Γ L P − Γ H P , and A f = (A H − A L )/(A H + A L ), where "H" and "L" denote two mass eigenstates with difference lifetimes, 1/Γ L P and 1/Γ H P . In experiment, a branching ratio is usually extracted from the total event yield. It means that the lifetime of neutral mesons is nothing to do with the measurement of branching ratio. Thus, the experimentally measurable branching ratio can be defined as On the other hand, the theoretical branching ratio is considered as Therefore, when we define (the normalization is adjusted as appropriate), where it can be derived with use of (B.27), (B.28), and (B.30). In the 2HDMs of Z 2 symmetric types and of aligned type with real ζ f , we trivially see φ P = φ S = 0. In this case, finally we can derive (3.8)