$R(D^{(*)})$ in a general two Higgs doublet model

Motivated by an anomaly in $R(D^{(*)})={\rm BR}(\bar{B}\rightarrow D^{(*)} \tau^-\bar{\nu})/{\rm BR}(\bar{B}\rightarrow D^{(*)} l^-\bar{\nu})$ reported by BaBar, Belle and LHCb, we study $R(D^{(*)})$ in a general two Higgs doublet model (2HDM). Although it has been suggested that it is difficult for the 2HDM to explain the current world average for $R(D^{(*)})$, it would be important to clarify how large deviations from the standard model predictions for $R(D^{(*)})$ are possible in the 2HDM. We investigate possible corrections to $R(D^{(*)})$ in the 2HDM by taking into account various flavor physics constraints (such as $B_c^-\rightarrow \tau^- \bar{\nu}$, $b\rightarrow s\gamma$, $b\rightarrow s l^+l^-$, $\Delta m_{B_{d,s}}$, $B_s\rightarrow \mu^+\mu^-$ and $\tau^+\tau^-$, and $B^-\rightarrow \tau^- \bar{\nu}$), and find that it would be possible (impossible) to accommodate the 1$\sigma$ region suggested by the Belle's result when we adopt a constraint ${\rm BR}(B_c^-\rightarrow \tau^- \bar{\nu})\le30~\%$ (${\rm BR}(B_c^-\rightarrow \tau^- \bar{\nu})\le10~\%$). We also study productions and decays of heavy neutral and charged Higgs bosons at the Large Hadron Collider (LHC) experiment and discuss the constraints and implications at the LHC. We show that in addition to well-studied production modes $bg\rightarrow tH^-$ and $gg\rightarrow H/A$, exotic productions of heavy Higgs bosons such as $cg\rightarrow bH^+,t+H/A$ and $c\bar{b}\rightarrow H^+$ would be significantly large, and the search for their exotic decay modes such as $H/A\rightarrow t\bar{c}+c\bar{t},~\mu^\pm\tau^\mp$ and $H^+\rightarrow c\bar{b}$ as well as $H/A\rightarrow \tau^+\tau^-$ and $H^+\rightarrow \tau^+\nu$ would be important to probe the interesting parameter regions for $R(D^{(*)})$.


I. INTRODUCTION
The standard model of elementary particles has been very successful to explain phenomena in nature. Recently, however, there are several observables in experiments, which may be suggesting the existence of new physics. For example, measurements of b → c transition processesB → D ( * ) l −ν at BaBar [1, 2], Belle [3][4][5] and LHCb [6,7] indicate a discrepancy between the experimental and theoretical values of where l = e, or µ. Measurements of b → s transition process B → K * µ + µ − at Belle [8], ATLAS [9], CMS [10],BaBar [11] and LHCb [12] also suggest an anomaly in the angular observable and moreover the LHCb has reported deviations from the SM predictions in R ( * ) [13]. A discrepancy between experimental and theoretical values of the muon anomalous magnetic moment has been also well-recognized [14]. These anomalies may indicate some hints of new flavor structure beyond the standard model.
One of the simplest extensions of the standard model is a two Higgs doublet model (2HDM). When both Higgs doublets couple to all fermions, it is notorious that the Higgs bosons have flavor violating interactions even at tree level, and it tends to induce the large flavor violating phenomena. Therefore, without any experimental supports, it would be difficult to justify this kind of extension. As we mentioned above, however, currently there are several indications which may support the extra flavor violation beyond the standard model, and hence we seriously consider a general type of two Higgs doublet model as a possibility to explain some of the anomalies mentioned above. In Refs. [15][16][17], it has been pointed out that µ−τ flavor violations in the general 2HDM may be able to explain the muon g −2 anomaly. 1 In this paper, we concentrate on R(D ( * ) ) in the general 2HDM.
1 Although recently the CMS collaboration reported a stronger constraint on the lepton flavor violating Higgs boson decay h → µτ [18], the µ − τ flavor violations in the general 2HDM can still explain the muon g − 2 anomaly, which is consistent with the recent CMS constraint, as seen in Refs. [16,17].
In order to study the anomaly, a model-independent effective operator approach [26][27][28][30][31][32] is very useful to understand what type of interaction is relevant to the anomaly. On the other hand, to identify the new physics model and to know constraints from other processes and possible correlation among other phenomena, it is necessary to study the anomaly in a model-dependent way. Furthermore, the anomaly would be explained by the new particles whose masses are of O(1) TeV scale. The current Large Hadron Collider (LHC) experiment is powerful to investigate such O(1) TeV scale physics as stressed in Refs. [33,34]. Therefore the interplay between the flavor physics and LHC physics would be also very helpful to probe a source of the anomaly.

II. GENERAL TWO HIGGS DOUBLET MODEL
In a two Higgs doublet model, when the Higgs potential is minimized in the SMlike vacuum, both neutral components of Higgs doublets get vacuum expectation values (vevs) in general. Taking a certain linear combination, we can always consider a basis (so called Higgs basis or Georgi basis [41,42], and see also, for example, [43][44][45][46][47]) where only one of the Higgs doublets has the vev as follows: where G + and G are Nambu-Goldstone bosons, and H + and A are a charged Higgs Here θ βα is the mixing angle. Note that when cos θ βα → 0 (sin θ βα → 1), the interactions of φ 1 approach to those of the SM Higgs boson. In this paper, we adopt Higgs basis. In Appendix A, we summarize a relation between the Higgs basis and the general basis.
If any discrete symmetries are not imposed, both Higgs doublets can couple to all fermions. 2 In the mass eigenbasis for the fermions, the Yukawa interactions are expressed by Here where where c βα ≡ cos θ βα and s βα ≡ sin θ βα . Note that when c βα is small, the Yukawa Here we also stress that the interactions of the charged Higgs boson are simply parameterized by ρ f Yukawa couplings in the Higgs basis. In order to analyze the effects on R(D ( * ) ), we adopt the Higgs basis in our analysis because it is convenient to effectively understand how large deviation from the SM prediction on the R(D ( * ) ) is possible in the 2HDM, which is the aim of this paper. To understand the effects on R(D ( * ) ), we can consider some simple flavor violation in the Higgs basis. On the other hand, if we consider the simple flavor violation in the original basis, it may correspond to the complex flavor violation in the Higgs basis (as shown in Appendix A), and hence it induces the effects not only on R(D ( * ) ) but also on other processes, which may generate strong constraints in the model. In that sense, we expect that our approach is conservative to see the possibly large effects on R(D ( * ) ), but consistent with other constraints. Therefore, in our analysis, using the Higgs basis, we try to clarify how large deviations are possible within the framework of the 2HDM in general.

A. Higgs mass spectrum
A scalar potential in the general 2HDM is given by In the basis shown in Eq. (2), Higgs boson masses are related as follows: Especially, when c βα is close to zero (or λ 6 ∼ 0), we approximately get the following expressions for the Higgs boson masses: Note that fixing the couplings λ i , the heavy Higgs boson masses are expressed by the CP-odd Higgs boson mass m A . We also note that a dangerous contribution to Peskin-Takeuchi's T-parameter are suppressed by the degeneracy between m A and m H + or m A and m H as well as the small Higgs mixing parameter c βα .
In the 2HDM, the charged Higgs boson generates a new contribution to B → Dl −ν . The relevant hadronic matrix elements are parameterized by the form factors f 0 (q 2 ) and f + (q 2 ) as follows [59,60]: where p B and p D (m B and m D ) are momenta (masses) of B and D mesons, respectively, and q is a momentum transfer , and m b and m c are b and c quark masses, respectively. In Appendix B, we summarize an information on form factors and in Appendix D, we list numerical values of various parameters we use in our numerical analysis.
Then the differential decay rate in the 2HDM is given by where The charged Higgs boson induces the corrections which are proportional to the form factor f 2 0 (q 2 ), and the contributions δ Dl H + are expressed by where m H + is a charged Higgs boson mass, and ρ f (f = e, u, d) are Yukawa couplings introduced in Eq. (4). The first term in Eq. (13) is an interference between the SM W-boson and charged Higgs boson contributions suppressed by m 2 H + , and the second term comes from the charged Higgs boson contribution suppressed by Notice that the charged Higgs boson contribution to Γ(B → Dl −ν ) is proportional to ρ ll e or (ρ † e ρ e ) ll . We also note that if (ρ † d V † − V † ρ u ) bc is small, the charged Higgs boson contribution is suppressed because the scalar couplingbcH + is small.
On the other hand, if (ρ † d V † − V † ρ u ) bc is sizable, the charged Higgs contribution can be significant especially in the large q 2 region.
The relevant hadronic matrix elements for this process are parameterized by the form factors V (q 2 ) and A i (q 2 ) (i = 1 − 3) [59,60]: where Here p D * and m D * are momentum and mass of D * , respectively, and q is a momentum In Appendix B, we summarize an information on the form factors.
The differential decay rate in the 2HDM is given by where Note that compared to the processB → Dl −ν , various form factors contribute toB → D * l −ν . However the charged Higgs boson only generates the corrections proportional to A 2 0 (q 2 ), and the corrections δ D * l H + (q 2 ) are given by Notice that similar toB → Dl −ν process, the charged Higgs boson contributions to the decay rate is proportional to ρ ll e or (ρ † e ρ e ) ll . We also note that contrary toB → Dl −ν process, the charged Higgs contributions are proportional to (ρ † d V † + V † ρ u ) bc since the pseudo-scalar couplingcγ 5 bH + contributes toB → D * l −ν process.
R(D ( * ) ) are defined to measure a lepton flavor universality between τ and l (l = e, µ) modes inB → D ( * ) l −ν : where l = e or µ. So far, an apparent flavor non-universality between µ and e modes has not been reported [61]. Therefore, we expect that the deviation from the SM prediction of R(D ( * ) ) mainly comes fromB → D ( * ) τ −ν mode in the 2HDM.
As we discussed in the previous section, the charged Higgs boson contributions toB → D ( * ) l −ν are proportional to ρ ll e or (ρ † e ρ e ) ll (l = e, µ, τ ). Therefore, ρ τ τ e , ρ eτ e and ρ µτ e induce the corrections toB → D ( * ) τ −ν , on the other hand, ρ iµ e and ρ ie e (i = e, µ, τ ) break the lepton flavor universality between µ and e modes, and hence here we assume ρ iµ e and ρ ie e are negligibly small. The charged Higgs boson contributions inB → D ( * ) τ −ν are proportional to Because of the hierarchical structure of the CKM matrix as effects of Yukawa couplings ρ bb d , ρ cc u and ρ uc u are small even if these are of the order of one. Therefore, we consider the effects of ρ tc u , ρ sb d and ρ db d . In order to see how large deviations of R(D ( * ) ) from the SM prediction are possible in the 2HDM, we consider the following three typical scenarios: 1). non-zero ρ tc u 2). non-zero ρ tc u and ρ sb d (and ρ db d ) 3). non-zero ρ tc u and ρ µτ e (and ρ eτ e ) in addition to the non-zero ρ τ τ e . We will also show the predictions in the Type II 2HDM as a comparison.
Before we study scenarios listed above, first we discuss the constraints from measurements of Higgs boson couplings, a top decay t → hc and lepton flavor violating Higgs boson decays on the relevant Yukawa couplings in these scenarios.
The result of the constraint on ρ τ τ e c βα at the 95% confidence level (C.L.) is given by

Constraint from the top decay t → hc
Non-zero ρ tc u as well as ρ ct u generate an exotic top quark decay t → hc. The decay branching ratio is obtained by Here we have adopted Γ t = 1.41 GeV for the total decay rate of the top quark and we assumed that ρ ct u is negligible in the numerical estimate of Eq. (28). The current experimental limit is set as at the 95% C.L. [64].

Constraints from lepton flavor violating Higgs boson decays
A search for the lepton flavor violating Higgs boson decays h → eτ, µτ constrains on the lepton flavor violating Yukawa couplings ρ eτ e and ρ µτ e . The latest data from the CMS collaboration at the LHC with an integrated luminosity of 35.9 fb −1 at √ s = 13 TeV is given by at 95% C.L. [18]. The branching ratio BR(h → µτ ) in the general 2HDM is given by where Γ h is the total decay rate of the Higgs boson and we have used Γ h = 4.1 MeV, and we have neglected ρ τ µ e for this numerical estimate. Similarly we get Note that if the Higgs mixing parameter c βα is small (|c βα | ≤ 10 −2 − 10 −3 ), these ρ f Yukawa couplings relevant to R(D ( * ) ) can be still of O(0.1) − O(1).
B. Scenario (1): non-zero ρ tc u and ρ τ τ e This scenario has been suggested by Refs. [26,35] as an interesting solution to explain the R(D ( * ) ) anomalies in the effective theory language. In this scenario, the charged Higgs boson contributions δ D ( * ) τ H + in Eqs. (13,18) are expressed by where double sign corresponds to a case for D and D * , respectively. For the parameterization of the Cabbibo-Kobayashi-Maskawa (CKM) matrix, we use the Particle Data Group (PDG) convention. We note that in the PDG convention, V cb and V tb are real and positive. When ρ τ τ e is real and negative, one can increase R(D * ) by increasing ρ tc u positively, but R(D) is decreased. However, as one increases ρ tc u further, R(D) also starts to increase so that the current world average values of R(D ( * ) ) can be explained. As noticed in Ref. [35], this scenario requires relatively large Yukawa coupling ρ tc u in order to explain R(D ( * ) ). However, the constraints from the other processes in this scenario would be important. Therefore, we will clarify the allowed parameter space consistent with the experimental results.
In general, the effective operators forB → D * τ −ν which are generated by the charged Higgs boson also contribute to B − c → τ −ν process. Therefore, large extra contributions toB → D * τ −ν process would be strongly constrained by this process [40,65]. 4 The branching ratio for B − c → τ −ν is obtained by where τ Bc is the lifetime of B − c meson and the charged Higgs boson contribution δ D * τ H + is given in Eq. (18). For non-zero ρ τ τ e and ρ tc u , the numerical estimate shows As studied in Ref. [40], the life time of the B − c meson provides a constraint BR(B − c → τ −ν ) ≤ 30%. Furthermore, the recent study suggested that the LEP data taken at the Z peak requires BR(B − c → τ −ν ) ≤ 10% [65]. We simply adopt both constraints as references. 5 In the Scenario (1), the constraint BR(B − c → τ −ν ) ≤ 30% (10%) corresponds to, for example, −0.8 ≤ ρ tc u ≤ 0.5 (−0.5 ≤ ρ tc u ≤ 0.2) for m H + = 500 GeV and ρ τ τ e = −0.5 if ρ tc u is real. In order to explain the current world average for R(D * ) [R(D * ) = 0.304], large Yukawa coupling ρ tc u is required (ρ tc u ∼ +1 for m H + = 500 GeV and ρ τ τ e = −0.5) [35]. Therefore, the constraint severely restricts the possibility to have a large deviation from the SM prediction for R(D * ).
In order to study the effects of the complex phases in the Yukawa couplings, we parameterize the Yukawa coupling ρ τ τ e as where ρ τ τ is a real parameter and the phase δ τ τ is assumed to be 0 ≤ δ τ τ ≤ π 2 . The phase of Yukawa coupling ρ tc u is effectively absorbed into the phase of δ τ τ in R(D ( * ) ), and therefore we assume that ρ tc u is a real (positive or negative) parameter. 6 The leading effect of δ τ τ comes from the interference between the SM and charged Higgs contributions, and it is proportional to "ρ tc u ρ τ τ cos δ τ τ " in R(D ( * ) ). Therefore, the 5 The constraint obtained in Ref. [65] relies on the theoretical value of BR(B − c → J/ψlν), which may be subject to debate. Therefore, we use both constraints BR(B − c → τ −ν ) ≤ 30% and 10% as references. 6 As we show in Appendix C, the phase of ρ tc u affects an imaginary part of B d,s −B d,s mixing. Therefore, the measurements of the imaginary part of B d,s −B d,s mixing have a potential to distinguish the origin of the phase. range of δ τ τ (0 ≤ δ τ τ ≤ π 2 ) is sufficient to see the possible predictions for the R(D ( * ) ) when ρ tc u ρ τ τ is allowed to be positive or negative. In Figure 1, the correlations between R(D * ) and BR(B − c → τ −ν ) are shown for δ τ τ = 0 (dashed blue line) and π 2 (solid orange line). When the phase δ τ τ is zero, the predicted value of R(D * ) can not be larger than 0.279 (0.263) because of constraints [40,65]. We found that when the phase δ τ τ is nonzero, the constraint on B − c → τ −ν becomes much stronger. For example, for δ τ τ = π 2 , the upper limit on . Note that the correlations do not depend on the detail of values of model parameters other than the phase parameter δ τ τ in the current scenario.
Because of the constraints on BR(B − c → τ −ν ), the predicted R(D * ) are strongly restricted. When ρ tc u is positive and the phase parameter δ τ τ is non-zero, the imaginary value of ρ τ τ e ρ tc u can increase R(D). As one increases δ τ τ and ρ tc u , R(D) can get larger before the value of R(D * ) reaches the constraint from B − c → τ −ν as one can see in Figure 2. The predicted R(D) and R(D * ) can not explain the current world average values. They can (can not) be within the 1σ region of the combined Belle Here we have fixed a charged Higgs mass to be 500 GeV (m H + = 500 GeV) and Yukawa coupling Therefore, even when we change m H + and ρ τ τ , the predicted values of R(D ( * ) ) (and Figure 2. For the LHC physics study discussed later, we choose two reference points shown by a black circle in Figure 2. As a reference point 1, we select a point with ρ tc u ρ τ τ = −0.2, δ τ τ = 5π 12 and m H + = 500 GeV, which satisfies the constraint BR( The reference point 1 predicts (R(D), R(D * )) = (0.361, 0.262). As a reference point 2, we choose a point with ρ tc u ρ τ τ = −0.156, δ τ τ = π 2 and m H + = 500 GeV, which satisfies the constraint As stressed in Refs. [26,31], the charged Higgs boson contribution significantly affects the q 2 distribution in dΓ(B → Dτ −ν )/dq 2 . In Figure  the difference from the SM distribution increases in large q 2 regions as the phase of ρ τ τ e (δ τ τ ) becomes larger. The BaBar [2] and Belle [3] have provided results of the q 2 distribution. However, the general 2HDM can not explain the BaBar's results of R(D ( * ) ), and the Belle's result for the q 2 distribution still has large uncertainty.
Furthermore, their analyses for the q 2 distribution measurement rely on the theoretical models. Therefore, here we do not explicitly impose a constraint from the However we would like to stress that the precise measurement of the q 2 distributions in dΓ(B → Dτ −ν )/dq 2 would have a significant impact on this scenario.
As one sees, the charged Higgs boson contributions with non-zero ρ tc u prefer the negative values. The measured values ∆m exp B d,s [19] are and the SM predictions ∆m SM B d,s (95 % C.L. region) [67] are given by Since the uncertainty of the measured values is very small in contrast to the SM predictions, we take into account the uncertainty of the SM predictions, and we apply the 95% C.L. allowed region to the 2HDM contribution as follows: Compared with Eqs. (42,43), interesting parameter regions (ρ tc u ∼ 0.2 − 0.5 for m H + = 500 GeV) in Scenario (1) are consistent with the current limit on ∆m B d,s .
The effective operators relevant to b → sγ and b → sl + l − are given by where l is a charged lepton l = e, µ or τ . We summarize the detail expressions induced by the charged Higgs boson via ρ u Yukawa couplings for C 7,8,9,10 in Appendix C. The charged Higgs boson contributions with non-zero ρ tc u are expressed by where x c = m 2 c /m 2 H + 1 for m H + = 500 GeV and various functions are defined in Appendix C and the approximate expressions are obtained as for x c 1. Here Q u and Q H + are charges of up-type quarks (+2/3) and a charged Higgs boson (+1). We note that a term proportional to G γ (x c ) in C 2HDM 9(l) originates from γ penguin contribution, on the other hand, terms proportional to G Z (x c ) in C 2HDM 9,10(l) come from Z penguin contribution. Therefore, C 2HDM 9,10(l) are universal to all lepton flavor l. We also note that in G γ (x) there is a log-enhancement when x is small.
are estimated as The low-energy Wilson coefficient C 7 (µ b ) is evaluated by using the renormalization group equation and the charged Higgs boson contribution C 2HDM where η = α s (µ)/α s (µ b ) and we have taken µ = m W , µ b = 5 GeV and α s (m Z ) = 0.118 in our numerical estimate. As discussed in Ref. [68], for example, the allowed regions of the new physics contributions C NP 7 (µ b ) from a global fit analysis of b → s transition observables are constrained as at 1σ level. Therefore, the interesting regions in Scenario (1) (ρ tc u ∼ 0.2 − 0.5 for m H + = 500 GeV) are consistent with this constraint.
For b → sl + l − process, the contributions to C 9,10 play an important role. The numerical estimates for the charged Higgs boson contributions with non-zero ρ tc u are given by where the Z penguin contribution is suppressed, and hence only C 2HDM 9(l) receives the non-negligible correction from the γ penguin contribution, which has a "log x c " enhancement.
As we discussed in the introduction, the current experimental situation on b → sl + l − processes is involved. In the angular observable of B → K * µ + µ − , the discrepancies between the SM prediction and the measured value have been reported. (See Refs. [8][9][10][11][12].) The LHCb results have also indicated deviations from the SM predictions in lepton flavor universality measurements R ( * ) . In Ref. [69], for example, a global fit analysis of angular observables in B → K * µ + µ − as well as branching ratios of B → K ( * ) µ + µ − and B s → φµ + µ − suggests that the best fit values of new physics contributions to C 9,10 are preferred to be away from zero, (C NP 9(µ) , C NP 10(µ) ) = (−1.15, +0.26). As shown in Eqs. (60,61), the interesting regions in Scenario (1) (ρ tc u ∼ 0.2−0.5 for m H + = 500 GeV) can not explain these discrepancies and the fit may be only slightly improved, compared to the SM predictions. 7 Since the anomalies in b → sl + l − are still subject to discuss, we do not impose the constraints from this process in our analysis.
In addition to these b → s(d) transition processes, in the presence of ρ tc u and ρ τ τ e , B s → τ + τ − process is induced. However, we have checked that the constraint from this process is still very weak, and hence we only list the expression of B s → τ + τ − process in Appendix C.
In addition to the non-zero ρ τ τ e and ρ tc u , we introduce other non-zero Yukawa couplings ρ sb d (and/or ρ db d ) in this scenario. This scenario has been suggested by, for example, Refs. [2,28,31] in the effective theory description.
To study the typical parameter space which is consistent with the other experimental data, we parameterize the Yukawa couplings as follows: Note that Therefore for non-zero r s (and/or However, as we will see below, the strong constraints come from not only as well as non-zero ρ τ τ e and ρ tc u , the branching ratio for B − c → τ −ν is given by the same expression shown in Eq. (36): From the constraint on BR(B − c → τ −ν ), the upper limit on ρ tc u is modified by about 1/|1 + a * |, compared with the Scenario (1). In Figure 4, the constraint from shown by horizontal dashed black lines (30% and 10%). As discussed in Scenario (1), we should note that R(D * ) can not be larger than 0.279 (0.263) The detail expression is shown in Appendix C, and numerical values of the parameters are listed in Appendix D, and here we have assumed that To satisfy the constraint shown in Eq. (46), the degeneracy between H and A (that is, small ∆ 2 ) has to be realized 8 . For example, ∆ 2 /m 2 H has to be much smaller than 10 −3 for m H 500 GeV and ρ sb d ∼ 0.1. Therefore, we simply assume ∆ 2 = 0 in our studies below. To satisfy the constraint, the coupling ρ sb d c βα has to be small. For example, for c βα = 0.01, ρ sb d ≤ 0.08. In Figure 4, the constraints from ∆m Bs are shown by long dashed pink line. The regions above the long dashed pink line are excluded.
For B d mass difference, one can get the expression from one for B s by replacing its flavor index s with d. To satisfy the constraint on ∆m B d shown in Eq. (46), the coupling ρ db d c βα also has to be small. For example, for c βα = 0.01, ρ db d ≤ 0.02. Since the constraint is very severe, ρ db d can not generate the significant contributions to R(D ( * ) ). Therefore, we do not include the effects of ρ db d in our studies.

Constraint from
The measured value [70,71] of branching ratio for The SM prediction [72] is given by which agrees with the current measured value, and hence the new physics contributions would be strongly constrained.
The branching ratio for B s → µ + µ − in the current scenario is given by Here note that the small Higgs mixing parameter c βα can suppress the 2HDM contribution. In Figure 4, the constraints from B s → µ + µ − are shown by dotted cyan lines in each r s case. The regions above the line are excluded.
In addition to B s → µ + µ − process, B s → τ + τ − process could also put a strong constraint on this scenario. In the presence of non-zero ρ sb d and ρ τ τ e , the 2HDM contributions to Wilson coefficients are given by where we have assumed that ρ sb d and ρ τ τ e are real. The current experimental bound on this process is given by the LHCb collaboration [73]: at 95% C.L. Since the SM prediction is (7.73 ± 0.49) × 10 −7 [74], we neglect it in our analysis below. The branching ratio is given by where Note that ∆ τ is small when the Higgs mixing parameter c βα is small (c βα ∼ 0.01) and ρ τ τ e ∼ O(0.1). The numerical estimate shows where we have assumed m A = m H . In Figure 4, the constraints from B s → τ + τ − are shown by dashed-dotted purple lines in each r s case. The regions above the line are excluded.
The charged Higgs boson can also contribute to B − → τ −ν process. 9 The current experimental value [3,78] is given by The decay branching ratio for B − → τ −ν in the 2HDM can be obtained from one for B − c → τ −ν by replacing its flavor index c with u. Especially, if non-zero ρ sb,db d are assumed (and ρ iu u = 0 for i = u, c, t), the branching ratio in the 2HDM is given by where the SM ratio and Y ub are expressed by Here, to calculate the SM prediction, we have used the best fitted values of the CKM matrix elements (|V ub | = 0.409 × 10 −2 ). Since |V us /V ub | 6 × 10 and |V ud /V ub | For example, for m H + =500 GeV and ρ τ τ e = −0.2, which corresponds to the same parameter set as in Figure 4, the constraint from B − → τ −ν is given by at 95% C.L. In Figure 4, the constraint from B − → τ −ν is shown by the solid green line. The regions above the line are excluded.
Since the experimental constraints are severe, the effects on R(D ( * ) ) from the Yukawa couplings ρ sb,db d in the Scenario (2) are very limited. Therefore, we will not study the Scenario (2) furthermore.
D. Scenario (3): non-zero ρ tc u , ρ τ τ e and ρ µτ e (and/or ρ eτ e ) The lepton flavor violating Yukawa couplings ρ µτ e and ρ eτ e can also affectB → D ( * ) τ −ν . Therefore we study their effects. The charged Higgs contributions are given by where double sign corresponds to R(D) and R(D * ), respectively. As discussed in the Scenario (1), negative ρ τ τ e and positive ρ tc u increase R(D * ), but decrease R(D). However, since the flavor violating contributions do not interfere with the SM contributions, the effects behave like the imaginary part of ρ τ τ e discussed in the Scenario (1), and always increase both R(D) and R(D * ). Here, to study the effects of ρ µτ e , we take m H + = 500 GeV and ρ τ τ e = −0.15 as a reference set of parameters. We parameterize the flavor violating coupling ρ µτ e as In Figure 5, the predicted values of R(D) and R(D * ) are shown in the Scenario A new interesting process is τ → µγ (τ → eγ) when both ρ τ τ e and ρ µτ e (ρ eτ e ) Yukawa couplings are non-zero. We will discuss the τ → µγ process in the next section.

τ → µγ
In the presence of non-zero ρ τ τ e and ρ µτ e , τ → µγ process is generated. As discussed in Ref. [16], not only 1-loop contribution but also Barr-Zee type 2-loop contribution are important. The 1-loop contribution tends to be suppressed when H and A are degenerate. On the other hand, the Barr-Zee type 2-loop contribution is affected by ρ tt u , which does not have an impact on R(D ( * ) ). In Figure 6, the predicted branching ratio BR(τ → µγ) is shown as a function of ρ tt u and r 1 (= ρ µτ e −0.27 = ρ τ τ e −0.1 ). Here we have set c βα = 0.001, m H/A = m H + = 500 GeV and r τ = |ρ µτ e /ρ τ τ e | = 2.7. We note that this parameter set is consistent with the constraint from h → µτ . As seen in Since the Yukawa coupling ρ τ τ e plays a crucial role to induce the important contribution to R(D ( * ) ), the heavy neutral Higgs boson search in H/A → τ + τ − mode at the LHC experiment would be important [33]. When the Higgs mixing parameter c βα is small, the gluon-gluon fusion production for heavy neutral Higgs bosons H and A relies on the interaction via ρ tt u , which does not affect R(D ( * ) ). Here we use HIGLU [79] to calculate the production cross section for H and A via the gluon-gluon fusion process at the NNLO and they are given by at √ s = 13 TeV for m H/A = 500 GeV.
For the LHC physics study in Scenario (1), we consider reference point 1 and 2, shown by black circle points in Figure 2. As discussed in Scenario (1), the reference point 1 is defined by a point with ρ tc u ρ τ τ = −0.2, δ τ τ = 5π 12 and m H + = 500 GeV, and the reference point 2 is ρ tc u ρ τ τ = −0.156, δ τ τ = π 2 and m H + = 500 GeV. The decay branching ratios of H and A are shown as a function of ρ tc u in Figure 8. 11 For solid lines (dashed lines), we have taken the reference point 1 (reference point 2).
Here we assume that m H/A = m H + and c βα = 0.001 as a reference parameter set for Scenario (1). We also fix ρ tt u = 0.1. 12 When we fix m H + (m H + = 500 GeV) and a product of the Yukawa couplings ρ tc u ρ τ τ e , the predicted values of R(D ( * ) ) are also fixed. In the reference parameter 11 When H, A and H + are not degenerate, two body decays such as H → AZ and H → H ± W ∓ may be allowed. However, we note that for the heavy Higgs boson like m H 500 GeV, nondegeneracy for such two body decays requires large Higgs quartic couplings λ 4,5 . 12 We note that the phase δ τ τ does not affect the decay H/A → τ + τ − . set, on the other hand, the decay branching ratios for H/A depend on ρ tc u (fixing a value of ρ tc u ρ τ τ e ) as shown in Figure 8. When ρ tc u is small (large), the dominant decay mode is H/A → τ + τ − (H/A → tc + ct) because ρ τ τ (ρ tc u ) is large. For ρ tt u = 0.1, a decay branching ratio for H/A → tt is always small. For a fixed value of ρ tc u , an absolute value of ρ τ τ e in the reference point 1 (|ρ τ τ e | = 0.2/ρ tc u ) is larger than one in the reference point 2 (|ρ τ τ e | = 0.156/ρ tc u ), and hence the branching ratio of H/A → τ + τ − (H/A → tc + ct) in the reference point 1 is larger (smaller) than one in the reference point 2, as seen in Figure 8.
In Figure 9 yellow shaded region (lower)), respectively. As can be seen from Figure 9, the Run 2 results already put a stronger constraint on ρ tt u . In Figure 10 [82,83].
In addition to the gluon-gluon fusion production for the neutral Higgs bosons (H and A), in the existence of the sizable ρ tc u Yukawa coupling, the production via ρ tc u (gc → t + H/A)) shown in Figure 11 would be important. We calculate σ(pp → t + H/A) via gc → t + H/A process by using the calchep [84] with CTEQ leading order PDF. The result is given by at √ s = 13 TeV for m H/A = 500 GeV. In Figure 10 A main difference between Scenario (1) and (3) is the existence of non-zero ρ µτ e coupling. Because of non-zero ρ µτ e , the size of ρ τ τ e coupling can be smaller than one for the Scenario (1) in the interesting regions for R(D ( * ) ). For the LHC study in Scenario (3), we consider two reference points, which are denoted by a black circle in Figure 5 for Scenario (3). As a reference point 1, parameters with ρ tc u ρ τ τ e = −0.075, r τ (= |ρ µτ e /ρ τ τ e |) = 2.7 and m H + = 500 GeV are set, and as a reference point 2, those with ρ tc u ρ τ τ e = −0.045, r τ (= |ρ µτ e /ρ τ τ e |) = 2.7 are taken. As a reference parameter set for Scenario (3), we also assume that m H/A = m H + = 500 GeV and c βα = 0.001. In Figure 14, the decay branching ratios of H/A are shown as a function of ρ tc u . Here we have taken the reference parameter set with the reference point 1 (solid lines) and the reference point 2 (dashed lines) as mentioned above. A value of ρ tt u is fixed to be −0.1. For this ρ tt u , BR(H/A → tt) is small (not shown in the figure). As we can see from the figure, BR(H/A → τ + τ − ) can be much smaller than one in the Scenario (1), on the other hand, BR(H/A → µ ± τ ∓ ) can be of order one. Since absolute values of ρ τ τ e and ρ µτ e in the reference point 1 (|ρ τ τ e | = 0.075/ρ tc u ) are larger than those in the reference point 2 (|ρ τ τ e | = 0.045/ρ tc u ) with the fixed value of ρ tc u , BR(H/A → τ + τ − ) and BR(H/A → µ ± τ ∓ ) (BR(H/A → tc + ct)) in the reference point 1 are larger (smaller) than those in the reference point 2.
In Figure 15

B. Productions for charged Higgs boson
In the presence of non-zero ρ tc u and ρ τ τ e , the charged Higgs boson production via cg → bH + shown in Figure 16, decaying into τ + ν (H + → τ + ν) would be very important. As pointed out in Ref. [34], this process is directly related toB → D ( * ) τ −ν process, and hence the search for this charged Higgs boson production and decay mode would be a crucial probe for the anomaly of R(D ( * ) ). The production cross section via cg → bH + is calculated by the calchep [84] with CTEQ leading order PDF, Therefore the production cross section would be significantly large. In Figure 17, decay branching ratios for charged Higgs boson H + are shown as a function of ρ tc u . Here we have assumed a reference parameter set with the reference point 1 (solid lines) and the reference point 2 (dashed lines) for Scenario (1) as mentioned in Figure 8. Here the coupling ρ tt u is also fixed to be 0.1. We note that since the ρ τ τ e coupling in the reference point 1 for Scenario (1) is larger than one for the reference point 2, BR(H + → τ +ν ) (BR(H + → cb, tb)) for the reference point 1 is larger (smaller) than one for the reference point 2.
We also note that in Scenario (3), in the presence of ρ µτ e coupling, the charged Higgs boson can decay via the ρ µτ e coupling. The decay product, however, is τ  Figure 18. In Figure 18, the size of cross sections σ × BR with BR(H ± ) = BR(H ± → τ ± ν), BR(H ± → cb, bc) and BR(H ± → tb, bt) are shown in solid green lines, dashed purple lines and dotted orange lines, respectively.
Here the same parameter set is taken as in Figure 10. For the reference point 2 in Scenario (1), as seen in Figure 17, the cross sections σ(pp → bH + +bH − ) × BR(H ± ) with BR(H ± ) = BR(H ± → τ ± ν) are slightly smaller, on the other hand, σ(pp → bH + +bH − ) × BR(H ± ) with BR(H ± ) = BR(H ± → cb, bc) and BR(H ± → tb, bt) are slightly larger than those for the reference point 1 with fixed values of ρ tc,tt u . As one can see from the figure, the cross sections are significantly large in the interesting regions for R(D ( * ) ).
As suggested in Ref. [34], for pp → bH + → bτ + ν process, the main SM background comes from pp → jW + → jτ + ν where j stands for a light quark or a gluon jet misidentified as a b-quark jet, and another background process is pp → bW + → bτ + ν. Although the cross sections of these background processes are also large, the p T distribution of τ for the signal process would be significantly different from those for the background events as discussed in Ref. [34]. In Figure 19, we show p T distributions of τ for the signal process pp → bH + → bτ + ν (signal) and a background process pp → jW + → jτ + ν (BG1) and another background process pp → bW + → bτ + ν (BG2). For BG1, we have shown numbers of events multiplied by 0.015 because it is motivated by the fact that the light-parton misidentification probability is 1.5% as studied in Ref. [34]. We have imposed p T > 20 GeV for all j, b and τ as motivated by experimental studies. We have taken m H + = 500 GeV, ρ tc u = 0.3 and ρ τ τ = −0.67 (the reference point 1 for Scenario (1)). As we can see from the figure, the p T distribution for the signal events would be significantly different from those for the background events. Therefore, we expect that even current data of the LHC has a potential to discriminate the signal events from the backgrounds.
In addition to the cg → bH + production process, the production processes bg → cH − and bg → tH − would be also important when we have non-zero ρ tc u and ρ tt u . The production cross sections via these processes are computed by the calchep [84], In Figure for m H + = 500 (600) GeV. Especially when H ± → cb, bc, this production mode receives constraints from searches for dijet resonances [89]. The constraints are only available for a resonance mass larger than 600 GeV, and the CMS limit at the √ s = 13 TeV using 27 fb −1 data is obtained by for the resonance mass to be 600 GeV. Here A is the acceptance for narrow resonances with the kinematic requirements |∆η ij | < 1.3 and |η| < 2.5 and we take A = 0.6 as suggested by the CMS collaboration [89]. From Eq. (93), the dijet production cross section in our 2HDM is given by for m H + = 600 GeV. Therefore, the regions with |ρ tc u | ≤ 1 are still consistent with this limit. 13 The improvement of the dijet resonance search limit in future would also have important impact on this scenario.
In conclusion, the searches for these exotic production and decay processes would be crucial to probe the scenarios for the anomaly of R(D ( * ) ). Here we only present interesting theoretical predictions for the production cross sections of heavy neutral and charged Higgs bosons at √ s = 13 TeV at the LHC. Apparently, more realistic detailed analyses for the LHC physics would be important. 13 Other charged Higgs boson production modes such as σ(pp → bH ± ) × BR(H ± → cb, bc) and σ(pp → cH ± ) × BR(H ± → cb, bc) also receive a constraint from the dijet resonance searches. However, the limit is weaker than one in Eq. (95), because of the smaller cross section.

VI. SUMMARY
The current experimental results have indicated a discrepancy between the measured values and the SM predictions of R(D ( * ) ). Since the current knowledge of the SM can not explain the discrepancy, it would be worth studying the possibilities of the extension of the SM. In this paper, we have studied the R(D ( * ) ) in a general 2HDM to clarify how large deviations from the SM predictions of R(D ( * ) ) are possible by taking into account various flavor physics constraints. We found that it is possible (impossible) to accommodate the 1σ region of the Belle's results if we . To obtain the large deviations, the masses of the heavy Higgs bosons in the 2HDM should be less than O(1) TeV if the size of new Yukawa couplings is less than O(1). Therefore, a study for a direct production of these heavy Higgs bosons at the LHC experiment is very important.
We have studied the productions and decays of the heavy Higgs bosons at the LHC, and discussed the constraints from the current LHC results and implications at the current and future searches at the LHC. Especially we found that the exotic productions such as cg → t + H/A and cg → bH + would be significantly large, and the searches for the exotic decay modes such as H/A → tc + ct, µ ± τ ∓ and H + → cb as well as H/A → τ + τ − and H + → τ + ν would be quite important to probe the interesting parameter regions which generate the sizable deviations from the SM predictions of R(D ( * ) ) in the general 2HDM. We have also shown that the p T distribution of τ lepton in signal events of pp → bH + → bτ + ν would be significantly different from those in the background events. Therefore, we expect that even the current data of the LHC would have a sensitivity to probe the interesting regions for R(D ( * ) ). Therefore, the interplay between the flavor physics and the LHC physics would play a crucial role to reveal the origin of the anomaly of R(D ( * ) ).

Note added
After we finished this work, we were aware of new (preliminary) LHCb result on R J/Ψ = BR(B c → J/Ψτν)/BR(B c → J/Ψµν) [90], which is another interesting measurement of lepton flavor universality in b → clν process. The result is R LHCb J/Ψ = 0.71 ± 0.17 ± 0.18 and suggests about 2σ deviation from the SM prediction R SM J/Ψ = 0.283 ± 0.048 [91]. Our study on R J/Ψ and R ηc = BR(B c → η c τν)/BR(B c → η c µν) (which is possibly another measurement of lepton flavor universality in b → clν process, and the SM prediction is R SM ηc = 0.31 +0.12 −0.07 [92]) in the general 2HDM shows that (R where i, j denote flavor indices and the summations are taken into account. Herẽ Φ a = iσ 2 Φ * a (a = 1, 2) where σ 2 is a Pauli matrix. We parameterize Higgs boson doublet fields as follows: where both Higgs boson doublets get vevs. We can always perform the following transformation: where we define a mixing angle β to satisfy so that only one of Higgs boson doublets (H 1 ) receives the vev (v) shown in Eq. (2). This is called Higgs (Georgi) basis. Here φ + 1,2 , G 0 1,2 , φ 0 1,2 in Φ 1,2 are related to G + , H + , G, A, φ 1,2 in H 1,2 (shown in Eq. (2)) via the transformation in Eq. (A3).
Under the Higgs basis, we rewrite Yukawa interactions in Eq. (A1) as follows: where Yukawa couplings Y f 1,2 are written in terms of original Yukawa couplings y f 1,2 Since only H 1 receives vev, masses of quarks and leptons are induced from the Yukawa couplings Y f 1 (f = d, u, e). Therefore, we diagonalize the Yukawa couplings for f = d, u, e in order to obtain the fermion's mass eigenbasis. The mass eigenbasis for the fermions (f L,R ) is defined by Note that the CKM (V ) and MNS (V MNS ) matrices are defined by V = U u L U † d L and V MNS = U e L U † ν L (where we assume that U ν L diagonalizes the 3 × 3 left-handed Majorana neutrino mass matrix in the low-energy effective theory), respectively.
Comparing these with the Yukawa interactions in Higgs basis (Eq. (2)), the relation between Yukawa couplings in both bases are obtained by where y f and ρ f (f = d, u, e) are defined in Higgs basis (Eq. (2)) and y f is diagonalized and it is related to the fermion mass y i f = √ 2m i f /v, on the other hand, y f 1,2 are defined in general basis shown in Eq. (A1).
We stress that any bases can be transferred to the Higgs basis. In general, if there is no symmetry to distinguish Higgs boson doublets Φ 1 and Φ 2 , it is difficult to define the most natural basis without considering the theory beyond the 2HDM.
However, it would be useful to discuss the relations between the original basis and Higgs basis in some specific cases. For well-known type I, type II, type X (leptonspecific) and type Y (flipped) 2HDM, for example, an extra Z 2 symmetry restricts the Yukawa interactions in the model, • Type I (y d 1 = y u 1 = y e 1 = 0), • Type II (y d 2 = y u 1 = y e 2 = 0), • Type X (lepton-specific) (y d 1 = y u 1 = y e 2 = 0), • Type Y (flipped) (y d 2 = y u 1 = y e 1 = 0).
As a consequence, the flavor violating Yukawa couplings are not allowed. However, in many cases, such a Z 2 symmetry may not be exact, and hence the symmetry breaking may induce the flavor violating Yukawa couplings. For example, in the Type II 2HDM, the exact Z 2 symmetry forbids the Yukawa couplings y d 2 , y u 1 and y e 2 , as shown above, and hence these Yukawa couplings may be induced as corrections to the Z 2 symmetry breaking. When we write these Z 2 breaking Yukawa couplings as ∆y A (A = d 2 , u 1 , e 2 ), the relations between the Yukawa couplings in the original basis and Higgs basis (Eq. (A9)) are written by for f = d, e and for f = u. Therefore, we can express the flavor violating Yukawa coupling ρ f in the Higgs basis in terms of the symmetry breaking Yukawa coupling ∆y A as follows: where y f is diagonalized in the Higgs basis and it is related to the fermion mass In this Appendix, we summarize hadronic matrix elements forB → D ( * ) lν which we use in our numerical analysis. The formula we use is taken from Refs. [59,60]. 14 The hadronic matrix elements relevant toB → Dl −ν in the 2HDM are written as where Here p D * and m D * are momentum and mass of D * , respectively, and q is a momentum Here the form factors f 0 , f + , V and A i (i = 1 − 3) are given by The q 2 dependence of these form factors comes through w D ( * ) (q 2 ) = w D ( * ) (q 2 ) + 1 + √ 2 .
where i = 1 and 2 for B d −B d and B s −B s , respectively. The contributions induced by the neutral Higgs boson mediations at the tree level are given by where y f φij is shown in Eqs. (6). The charged Higgs boson generates the contribution to C V LL via ρ u Yukawa couplings at the one loop level as follows: (C10) The numerical values of the decay parameters are listed in Appendix D.
2. b → sγ and b → sl + l − (l = e and µ) The effective operators relevant to b → sγ and b → sl + l − are given by e 2 16π 2 C 9(l) (sγ µ P L b)(lγ µ l) + e 2 16π 2 C 10(l) (sγ µ P L b)(lγ µ γ 5 l) , where l is a charged lepton l = e or µ. The charged Higgs boson contributions with non-zero ρ u are expressed by where Q u,H + ,l are electric charges of up-type quark, charged Higgs boson and charged lepton (Q u,H + ,l = 2/3, +1, −1), respectively and T 3l = −1/2, and various functions are given by We note that a term proportional to G γ1, γ2 in C 2HDM

9(l)
originates from γ penguin contribution, on the other hand, terms proportional to G Z in C 2HDM 9,10(l) come from Z penguin contribution. Therefore, C 2HDM 9,10(l) are universal to all lepton flavor l. We also note that in G γ1 (x) there is a log-enhancement when x is small. In a general 2HDM, the Yukawa couplings ρ sb(bs) d induce B s → l + l − process at the tree level, and furthermore non-zero ρ u Yukawa couplings also generate B s → l + l − process at the one loop level.
The decay rate is given by (C26) At the one loop level, the charged Higgs boson via non-zero ρ u Yukawa couplings can induce the contribution to C bs A(l) which is the same as the effective operator proportional to C 2HDM 10(l) discussed in b → sl + l − process [Eq. (C15)]: which is induced by Z penguin contribution. We note that the expression in Eq. (C15) is correct for l = τ because it is lepton flavor universal. In addition, in the presence of ρ τ τ,µτ e Yukawa couplings, the box diagram generates the following contribution to C bs A(τ ) : where x i = m 2 u i /m 2 H + and the function B(x) is defined by