Prospects of light charged scalars in a three Higgs doublet model with $Z_3$ symmetry

The stringent constraints from the direct searches for exotic scalars at the LHC as well as indirect bounds from flavor physics measurements have imposed severe restrictions on the parameter space of new physics models featuring extended Higgs sectors. In the Type-II 2HDM, this implies a lower bound on the charged Higgs masses of $\cal O$(600 GeV). In this work we analyze the phenomenology of a Z3HDM in the alignment limit focusing on the impact of flavor physics constraints on its parameter space. We show that the couplings of the two charged Higgs bosons in this model feature an additional suppression factor compared to Type-II 2HDM. This gives rise to a significant relaxation of the flavor physics constraints in this model, allowing the charged Higgs masses to be as low as $\cal O$(200 GeV). We also consider the constraints coming from precision electroweak observables and the observed diphoton decay rate of the 125 GeV Higgs boson at the LHC. The bounds coming from the direct searches of nonstandard Higgs bosons at the LHC, particularly those from resonance searches in the ditau channel, prove to be very effective in constraining this scenario further.


Introduction
Two of the major tasks that will be undertaken in the upcoming runs of the LHC and beyond comprise of the precise determination of the properties of the 125 GeV Higgs boson as well as direct searches for additional scalar particles. The remarkable consistency between the predictions from the standard model (SM) and the experimental data from the LHC so far has posed strong challenges for new physics (NP) scenarios beyond the standard model (BSM). A complimentary pathway to explore NP is provided by the low energy precision measurements in flavor physics. Measurements from dedicated flavor physics experiments like BELLE, BABAR and LHCb has so far been largely in agreement with the SM, providing stringent constraints on most of the BSM scenarios.
Introduction of additional Higgs doublets has been one of the most popular choices for new physics extensions beyond the SM. The most minimal choice, the two Higgs doublet model (2HDM) [1,2], has been discussed widely in the literature from both theoretical and phenomenological points of view. In the well-known alignment limit, the lightest CP-even Higgs boson of 2HDM can be SM-like in its tree-level couplings to the fermions and vector bosons and thus can serve as the 125 GeV scalar observed at the LHC [3][4][5][6][7][8]. The additional (pseudo) scalar and charged Higgs bosons can give rise to interesting signatures at the LHC as well as at various flavor physics experiments. Consistency with the strong constraints from the LHC and flavor observables often pushes the charged Higgs boson mass in 2HDM towards the heavier end of the spectrum. It has been observed that a combination of flavor physics measurements can exclude the charged Higgs masses below O (600 GeV) in 2HDM of Type-II [9], where up and down-type quarks obtain their masses from two different Higgs doublets. This bound on the charged scalar masses can be somewhat relaxed in Type-I 2HDM where a single Higgs doublet is responsible for generating masses of the up and down type quarks [10]. This is because in Type-I 2HDM all the fermionic couplings of the charged scalar are proportional to cot β, with tan β being the ratio of the two vacuum expectation values (VEVs), as conventionally defined in 2HDMs. Therefore, the constraints on the nonstandard scalars can be easily evaded by choosing tan β 1. In this work we investigate the possibility of allowing lighter nonstandard scalars without compromising the essential feature of Type-II 2HDM i.e. two different doublets give masses to up and down quarks.
Moving beyond 2HDM, the most natural step ahead is to add one additional Higgs doublet, leading to the three Higgs doublet model (3HDM) [11][12][13][14][15][16]. As in the 2HDM case, it is possible to achieve an alignment limit corresponding to a physical scalar resembling the properties of 125 GeV SM-like Higgs boson [17,18]. In contrast to 2HDM, the scalar spectrum is much broader here, offering a rich phenomenology in both high and low energy experiments. Most importantly, the presence of additional nonstandard Higgs bosons leads to significant modifications in the flavor changing neutral and charged current processes compared to 2HDM. In Ref. [17], it was shown that the conditions for alignment limit in 3HDM can be parametrized by a set of simple equations closely mimicking those of the CP conserving 2HDM. Using the example of Z 3 -symmetric 3HDM (Z3HDM), it was observed that the analytic conditions can be easily implemented in a realistic scenario, making way for efficient numerical analysis. In the present work we explore the phenomenological aspects of the alignment limit in Z3HDM with an emphasis on the effects of flavor physics constraints on its parameter space. We show that the constraints on the parameter space stemming from the interplay of various flavor physics data are notably relaxed compared to those in the Type-II 2HDM. Such a relaxation of constraints transpires from the presence of an additional suppression in the couplings of the charged Higgs bosons in the model compared to Type-II 2HDM. We prescribe a simple analytical set up that automatically guarantees agreement with the ρ-parameter constraints as well as bounds arising from the measurement of Higgs to diphoton decay rate at the LHC. We also study the effect of the bounds coming from direct searches for additional Higgs bosons at the LHC on the parameter space of out interest.
Our paper is organized as follows. In Sec. 2 we describe the scalar sector and Yukawa structure of Z3HDM. The constraints from flavor physics observables are analyzed in Sec. 3. We calculate the diphoton decay rate for this model in Sec. 4. The limits coming from the direct searches at the LHC are discussed in Sec. 5. Finally, we summarize our results in Sec. 6.
fermion Type-I Type-II Type-X Type-Y Democratic with a Z 3 symmetry 1 . By doing so, we are able to find suitable charge assignments for both the fermions and the scalar doublets such that the NFC model ensues.
In our current setup, only scalars and right-handed fermionic fields may transform nontrivially under the Z 3 symmetry. Namely, we require where ω = e 2πi/3 , and n R ( R ) are the right-handed down-type quarks (leptons), as to clearly distinguish between the flavor and mass eigenstates. By taking all other fields to transform trivially under the Z 3 symmetry, it becomes clear that φ 1 couples to the leptons, whereas φ 2 and φ 3 couple to the down-and up-type quarks, respectively. As such, we achieve a Yukawa structure which ensures NFC and the absence of FCNCs at tree-level.

Scalar sector
While there are more than one ways to achieve a democratic Yukawa structure, the different choices will lead to different scalar potentials. The most general scalar potential for a Z3HDM obeying the symmetry in Eq. (1a) is given by [17,[24][25][26] where we allow the presence of the soft-breaking terms m 2 ij , i = j, since they will be of some importance for regulating the charged Higgs contribution to the diphoton decay amplitude [27]. For simplicity, we take all the parameters of the scalar potential to be real, so that the neutral scalars can be easily classified as CP-even and CP-odd bosons.
After the electroweak symmetry breaking (EWSB), the scalar doublets can be decomposed in terms of the component fields as where v k denotes the VEV of the field φ k ( φ k = v k / √ 2). For notational convenience, the VEVs are expressed as where v = v 2 1 + v 2 2 + v 2 3 is the usual electroweak (EW) VEV. The inclusion of three scalar doublets will give rise to four charged scalar particles, H ± 1,2 , three CP-even neutral ones h, H 1,2 , as well as two CP-odd neutral particles A 1,2 , where the remaining fields are the usual Goldstone bosons w ± , ζ. These physical particles can be obtained by rotating the fields onto the mass basis. For the charged and pseudoscalar sectors, we can obtain the physical scalars by performing the following 3 × 3 rotations, where, the rotation matrices are defined as and For the CP-even sector, we can obtain the physical mass basis through with For more details on the analysis of the scalar potential, we refer the reader to Appendix A.
The existence of nonstandard neutral CP-even scalars in the Z3HDM, in general, leads to a deviation of the couplings of the physical scalar h from the respective SM predictions. However, the data obtained from the LHC runs shows a good agreement of the experimental data to the SM prediction for the Higgs signal strengths [28,29]. This motivates us to work in the alignment limit which is a set of conditions such that the lightest CP-even scalar mimics the SM-Higgs in its tree-level couplings, automatically respecting the agreement between the experimental data and the corresponding SM predictions for the Higgs signal strengths. For our Z3HDM, the conditions for alignment are given by [17] As more data accumulate in the future runs of the High-Luminosity LHC (HL-LHC), the possibility of deviating from the alignment limit will become increasingly constrained, if no BSM signals are detected.

Yukawa sector and charged Higgs couplings
The quark Yukawa Lagrangian of the Z3HDM can be written as where Q L ≡ (p L , n L ) T denotes the SU (2) L left-handed quark doublet field, p R the up-type right-handed quarks, and Y d,u are the respective 3 × 3 Yukawa matrices in flavor space. After EWSB, the mass matrices of the down and up-type quarks are given by As usual, we can redefine the quark fields to rotate into the mass basis through which, in turn, will diagonalize the mass matrices through the bi-unitary transformation Similar to the SM, the CKM matrix is defined as V = U L D † L . As intended, our model does not have any FCNC at the tree-level, and the Higgs signal strengths will also be compatible with the corresponding SM expectations in the alignment limit. However, the presence of charged scalars brings forth new channels for loop contributions to several flavor observables such as neutral meson oscillations and b → sγ. In fact, these processes are quite restricted from experiments, and thus are usually used to place lower bounds on the nonstandard scalar masses, as their contributions must be kept in check. Thus, it becomes important to study the charged scalar couplings to the fermions, as these will govern the vertices responsible for these processes at the one-loop level.
Given its importance, we focus on the original quark Yukawa Lagrangian containing the charged Higgs couplings, where, in the last step, we have rotated into the quark mass basis. Our goal is to arrive at couplings among the physical fields, as so we further rewrite the Lagrangian in the scalar mass basis. Using which describes the vertices between the physical charged scalars to the physical quarks. The same process can be repeated to obtain the leptonic couplings: where, ≡ (e, µ, τ ) T , ν ≡ (ν e , ν µ , ν τ ) T and D = diag(m e , m µ , m τ ). In the following, we will focus mostly on the consequences of quark flavor observables. Hence, to better grasp the model's implications, it is helpful to substitute the X ij elements explicitly following Eqs. (6) and (7), recasting the charged Higgs couplings to quarks as One noteworthy observation is the similarity between the Z3HDM and the type-II 2HDM. In fact, both are NFC models, where the difference lies in the fact that, in the Z3HDM, the lepton Yukawa couplings have a dedicated doublet, whereas in the type-II 2HDM the leptons share the doublet responsible for the down type quark masses. The resemblance can be made more explicit by noting that due to the Z 3 charge assignments of the scalar doublets, φ 1 is responsible for the lepton masses, which are generally much lower than the quark masses. Combining this with the relation between each individual VEV and the EWSB seen in Eq. (4), it seems reasonable to assume v 1 v 2 , v 3 , which is achieved by taking large values of tan β 1 , while still remaining in a perturbative regime for the τ -Yukawa coupling. In this regime, where cot β 1 1, the cot β 1 dependency of the charged Higgs couplings of Eqs. (18a) and (18b) can be neglected, and the couplings become similar to those of the type-II 2HDM, relaxed by either cos γ 2 or sin γ 2 , which are always less than one. Indeed, by comparing with the corresponding couplings in the type-II 2HDM [1], we can identify tan β of 2HDM-II with tan β 2 of Z3HDM, since both control the ratio v u /v d , where v u (d) are the VEVs of the scalars that couple to the up (down) quarks, respectively. If we further consider a scenario where either H ± 1 or H ± 2 is relatively heavy ( 5 TeV), while keeping the other relatively light ( 1 TeV), then the heavy particle decouples and its contribution will be negligible, and our effective theory becomes similar to a type-II 2HDM scenario. The striking difference is that while one of the scalars is decoupled, the effective theory still retains some consequences of the full theory. In order to exemplify, we consider a scenario where H ± 2 is decoupled and tan β 1 1. In this case the H ± 1 couplings of Eq. (18a) can be approximated as: Comparing with Eq. (19), we notice the remarkable similarity with the type-II 2HDM except for the fact that the couplings are reduced in strength by a factor of sin γ 2 . This will play an important role in diluting the constraints from flavor data compared to those in the type-II 2HDM, which we will discuss in the next section.

Constraints from flavor data
Since compliance with flavor data is continuously pushing the lower bound on the mass of the charged Higgs of the type-II 2HDM upwards, the relaxation due to γ 2 in this effective 2HDM can easily justify lower masses for new charged particles, while still remaining within the experimental limits for the new physics (NP) contributions to the flavor processes.
In order to make the discussion concrete, we analyze the resulting bounds coming from flavor data. We restrict ourselves to the analysis of the NP contributions to the radiative decay b → sγ, as well as the bounds coming from the B meson oscillations, ∆M B s,d 2 . We make use of the FlavorKit [30] functionalities within SPheno [31,32], compiled by SARAH [33], explicitly retaining contributions up to one-loop only. In order to gain some qualitative insights into the processes and phenomenologies at hand, we refer the reader to appendix B, where we provide analytic expressions for the relevant processes. It is, however, easy to note that in models with no tree-level FCNCs, the only one-loop NP contributions to both b → sγ as well as ∆M B s,d will come from the charged Higgs couplings. Therefore, these observables will be governed by a set of five parameters, namely, As we mentioned earlier, the Z3HDM where one charged Higgs is decoupled from EW scale dynamics becomes a relaxed type-II 2HDM effective scenario. Namely, a remnant of the full theory survives as a damping of the usual type-II 2HDM charged scalar couplings, which will in turn result in a relaxation of the bounds that are found for the type-II 2HDM. As such, we initially focus on this case where one of the charged scalars is decoupled, featuring the relaxation of the bounds.
Our point is clearly exemplified in Fig. 1 where we note that the type-II 2HDM bounds coincide with the more restrictive case of this Z3HDM limit (γ 2 = π/2 for the bounds on H ± 1 , and γ 2 = 0 for H ± 2 ). As we can see, for our benchmark of tan β 1 = 10, the constraints on the charged scalar masses are, at worst, comparable to the corresponding bounds in type-II 2HDM for appropriate values of γ 2 . But the important point is that by changing the values of γ 2 the bounds can be considerably diluted. Even while keeping away from the extremal cases, the bounds can be easily relaxed by a factor of 2, by taking γ 2 = π/4, as clearly seen in the plots. From Fig. 1 we also note that there is an asymmetry in the bounds on H ± 1 and H ± 2 when we are away from the type-II 2HDM limit. This feature can be attributed to the tan β 1 dependency of the charged-Higgs couplings. Now that we have established that relatively light charged scalars can successfully pass through the stringent constraints imposed by the flavor data, it will be interesting if we can say something about the masses of the neutral nonstandard scalars in relation to those of the charged scalars. This is where the constraints from the electroweak ρ-parameter become useful. The neutral scalars are expected to have masses such that the impact of NP on the ρ-parameter is minimized. Using the general expressions in Refs. [34,35] we have calculated the NP contribution to the ρ-parameter in the alignment limit of our model. The relevant expression is particularly clean and intuitive in the limit γ 1 = γ 2 = −α 3 = α (say) as we display below: where, One easy way to circumvent the ρ-parameter constraint will be to impose m H + 1 ≈ m H1 ≈ m A1 = M 1 (say) and m H + 2 ≈ m H2 ≈ m A2 = M 2 (say) as ∆ρ becomes zero in this limit. Under this assumption, the scalar spectrum conveniently breaks down into two degenerate tiers of nonstandard masses. This spectrum of masses and mixings can be easily achieved with a simplified scalar potential of the following form, which has an enhanced symmetry in its quartic part [36]: In the above potential there are seven parameters which can be traded in favor of the seven physical parameters, (v, β 1 , β 2 , m h , M 1 , M 2 , α). The relevant reparametrizations are given below: where c x and s x are shorthands for cos x and sin x respectively.
At this point, we wish to remark that the potential of Eq. (23) contains only one quartic parameter, λ. Thus, both unitarity and stability of the scalar potential can be ensured by requiring 0 < λ < 4π 3 . Looking at Eq. (24d), we can easily see that the potential of Eq. (23) is manifestly compatible with the unitarity and vacuum stability constraints. Next, we extract the top, bottom, and τ Yukawa couplings as which follow from our convention that φ 3 , φ 2 , and φ 1 couple to up-type quarks, down-type quarks, and charged leptons respectively. For the perturbativity of Yukawa couplings, we should have |y t |, |y b |, |y τ | < √ 4π. The resulting constraint from perturbativity has been displayed in Fig. 3. Throughout our paper, we have used values of tan β 1,2 which are consistent with this perturbative region.

Implications for diphoton decay rate
At this point one might naturally wonder whether such light charged scalars would leave observable imprints in loop induced Higgs decays such as h → γγ. After the 13 TeV run of the LHC, updated constraints on the Higgs to diphoton signal strength has been reported by both the ATLAS [38] and CMS [39] collaborations at 139 fb −1 luminosity. It is thus important that we check whether such light charged scalars can negotiate the bound arising from the measurement of the Higgs to diphoton signal strength. To do that, we need to calculate the hH + i H − i couplings which, for the potential of Eq. (23) are given below: Using this, we can easily write down the expression for the diphoton signal strength as follows: where, and the loop functions are given by [40], with, f (x) = sin −1 1/x 2 for x > 1. It is interesting to note that in the limiting potential of Eq. (23), the charged Higgs couplings to the SM-like Higgs in Eq. (26) are completely independent of any mixing angles and fixed to a constant value. Therefore, the charged Higgs contribution to the loop-induced Higgs to diphoton channel will always be suppressed by the charged Higgs masses when the charged Higgses are much larger than the SM-like Higgs. We display our results in Fig. 2  ) mass plane, where we see that the current Higgs data mainly discards the parameter space where both or any one of the charged Higgs masses are below O (200 GeV). In Fig. 2, the region below the red dashed line is excluded by the current data at 95% C.L [38].

Direct search constraints
The presence of two charged and four additional neutral Higgs bosons places this model under the scrutiny of direct searches for nonstandard Higgs bosons at the LHC. In the parameter region favored by the flavor physics constraints (m H + 1,2 > m top ), the dominant production mode of a charged Higgs boson at the LHC is in association with tb-quark pairs. Both ATLAS and CMS collaborations have provided model-independent upper bounds on the production cross-section times branching ratio for this mode with the charged Higgs boson decaying to tb [41,42] and τ ν [43,44] final states 4 . On the other hand, the search for heavy scalar and pseudoscalar resonances yields the most stringent constraints in the τ + τ − final state. In this case model-independent bounds are available for production via gluon-gluon-fusion and in association with a b-quark pair [45,46]. In the following, we discuss the impact of the various bounds mentioned above on the parameter space of Z3HDM.   To elucidate the relevance of different direct search constraints on the parameter region of our interest, we plot in Fig. 4 the branching ratios of H ± 1 , A 1 and H 1 as a function of tan β 2 . Keeping in mind the precision constraints from electroweak ρ-parameter, we choose to work in the limit m H + i = m Hi = m Ai , i = 1, 2 and γ 1 = γ 2 = −α 3 . Furthermore, we consider the case of one decoupled charged Higgs boson for simplicity, which in this case we take to be H ± 2 5 . The branching ratios are calculated for a fixed tan β 1 = 10 and two different values of sin γ 2 shown as solid (γ 2 = π 6 ) and dashed (γ 2 = π 4 ) lines. For tan β 2 1 the leading decay mode of H ± 1 is H ± 1 → tb, because of the dominance of the first term in Eq. (18a), which is proportional to cot β 2 . As tan β 2 increases, the terms proportional to tan β 2 in Eqs. (17) and (18a) promptly takes over. In the tan β 2 > 1 region favored by the flavor physics data, H ± 1 → τ ν becomes the dominant decay mode. A somewhat similar pattern is observed in the branching ratio of A 1 and H 1 , with A 1 /H 1 → τ + τ − being the dominant decay mode in the tan β 2 > 1 region.
The implications of various direct search constraints on the parameter space of Z3HDM becomes even clearer by looking at Fig. 5 where we show the production cross-section times branching ratio of the Higgs bosons H ± 1 , A 1 and H 1 at the 13 TeV LHC as a function of their masses. For this calculation, we implemented our model in 4 The vanishing H ± 1 W ∓ h coupling in the alignment limit leads to the absence of W ± h final state in H ± 1 decay. Furthermore, the decay to W ± A 1 final state is kinematically disfavored as a result of the assumed degnearcy between H ± 1 and A 1 . 5 Similar bounds may be imposed on H ± 2 in the case of a decoupled H ± 1 .
FeynRules [48,49] to generate files in the UFO format [50]. These files are then used by MadGraph5 aMC@NLO [51] to compute the signal cross-section at the LHC. The gray-shaded region denotes the parameter space excluded by the corresponding bound from ATLAS. We consider a relatively small value of tan β 2 , namely, tan β 2 = 2, to comply with the stringent bounds from the LHC. The value of tan β 1 is kept fixed at tan β 1 = 10. As can be seen from Fig. 5a, the value of σ H ± 1 →τ ν tbH ± 1 remains comfortably within the upper limit set by the direct search for charged Higgs boson mass within the region of interest. Thus, even the relatively low m H + 1 region allowed by the flavor physics constraints remains effectively safe from the direct LHC constraints on charged Higgs bosons. However, the bounds from scalar and pseudoscalar resonance searches in the ditau channel can be more constraining in this case, especially when taken together with the indirect bounds from the precision measurements of the electroweak ρ-parameter and diphoton decay rate of the 125 GeV Higgs boson as discussed in Sec. 4. It can be clearly seen from Fig 5b that the associated production of the pseudoscalar A 1 with b-quark pairs does not impose any significant restrictions on the relevant parameter space. However, the production cross-section via gluon-gluon-fusion (ggF) process, as shown in Fig. 5c, can be significantly larger in this case, allowing only m A1 450 GeV. The similar LHC bounds on H 1 comes out to be weaker than those on A 1 in most of the parameter space for this scenario.  Before we conclude, it should be noted that this seemingly strong constraint from the direct searches should be interpreted with some caution. The strong constraints on the neutral bosons are essentially due to relatively high BRs into the τ τ channel, which, in turn, may be attributed to our choice of tan β 1 1. But one should also remember that we have been working in a simplified limit of Z3HDM where the nonstandard scalars come with two tiers of degenerate masses, which we motivated from the ρ-parameter constraints. However, we can lift the degeneracy and allow one of the neutral scalars to have a different mass while still keeping the NP contributions to the ρ-parameter under control. In this way, it will be possible to open up channels like A 1 → H 1 Z which will reduce BR(A 1 → τ τ ). Moreover, there can be additional decay modes in the scalar sector too. As an example, if there is a dark singlet coupling to the other scalars, then decay modes such as H 1 → SS where S is the dark singlet, can open up. Keeping all these possibilities in mind, we can say that the bounds from the direct searches in our simplified analysis can be considerably relaxed.

Summary
To summarize, we have analyzed a 3HDM with Z 3 -symmetry featuring NFC, where the Z 3 -symmetry ensures a democratic Yukawa structure requiring each type of SM fermion to be coupled to a particular Higgs scalar doublet, thus eliminating FCNCs at tree level. We have discussed the characteristics of the scalar and Yukawa sectors in detail, focusing on the alignment limit where the lightest CP-even Higgs boson of the model possesses SM-like tree-level couplings and hence can serve as a candidate for the 125 GeV scalar observed at the LHC. This alignment limit can be characterized by a set of simple analytic conditions closely resembling that of a 2HDM.
The alignment limit is also phenomenologically well-motivated in view of the increasingly precise measurements of the signal strengths of the 125 GeV SM-like Higgs boson at the LHC.
The presence of two pairs of charged (H ± 1,2 ) and additional two neutral CP-odd (A 1,2 ) and two neutral CP-even (H 1,2 ) Higgs bosons in the model gives rise to distinctive signatures in various experiments looking for direct or indirect signals of BSM physics. From the phenomenological point of view, we have put an emphasis on analyzing the effect of the flavor physics constraints on the parameter space of Z3HDM. The leading BSM contribution to flavor observables like BR(b → sγ) and ∆M B s,d comes from the loops containing the charged Higgs bosons H ± 1 and H ± 2 . The Yukawa coupling structure of the charged Higgs bosons in this model bear close likeness to those of the Type-II 2HDM. However, the key difference from Type-II 2HDM is that the fermionic couplings of the charged Higgses feature an additional suppression effect which is essentially non-decoupling in nature. Thus, even in the limit of an effective Type-II 2HDM, with one of the charged Higgs bosons taken to be decoupled from the spectrum, the couplings of the other charged Higgs retains the suppression factor. This produces a significant relaxation of the bounds coming from flavor observables in this model compared to Type-II 2HDM. It is observed that charged Higgs masses as low as 200 GeV are allowed by the flavor data in Z3HDM, whereas in the case of Type-II 2HDM the lower bound on charged Higgs mass from the same flavor physics constraints stands at O (600 GeV). We have discussed the combined contribution of H ± 1 and H ± 2 to flavor observables when both of them are taken to be light. We have also taken into account the precision constraints from EW ρ-parameter which can be easily satisfied in a simple set up with two sets of degenerate masses for the nonstandard scalars, m H + i = m Hi = m Ai , i = 1, 2. We show that in this limiting case, the couplings of the charged Higgs bosons to the 125 GeV SM-like Higgs assumes a constant value. Therefore, the contribution to h → γγ decay from charged Higgs loop, being suppressed by a factor of m 2 h /m 2 H + i , i = 1, 2, does not produce any additional constraint on the relevant parameter space.
Finally, we have analyzed the constraints coming from the direct searches of the nonstandard Higgs bosons at the LHC. We show that the bounds from direct charged Higgs boson searches can be satisfied relatively easily in most of the parameter region satisfied by the flavor data. However, the constraints on neutral CP-even and CP-odd Higgs boson masses coming from ditau resonance searches put somewhat stringent bounds on our parameter space. We also note that the strong constraints from the ditau resonance searches are a consequence of the particular parameter choice we make to satisfy the bounds from EW ρ-parameter. Some alternative choice of parameters to satisfy the ρ-parameter constraint may lead to a significant reduction in the branching ratios of the neutral CP-even and CP-odd scalars to τ τ final states. This, in turn, may lead to a dilution of the LHC bounds by a considerable amount.

A Appendix
In this Appendix we give a detailed description of the scalar sector. It is to be noted that the scalar potential of Eq. (2) contains 18 parameters including the three soft-symmetry breaking terms m 2 12 , m 2 13 and m 2 23 . Among all, the bilinear parameters m 2 11 , m 2 22 and m 2 33 can be traded for the three VEVs, v 1 , v 2 and v 3 or equivalently v, tan β 1 and tan β 2 .The remaining twelve quartic couplings will correspond to the seven physical masses (three CP-even scalars, two CP-odd scalars and two pairs of charged scalars) and five mixing angles (three in the CP-even sector, one in the CP-odd sector and one in the charged scalar sector). Below, we demonstrate this relations by examining the potential of Eq. (2) in more detail.
The minimization conditions used to replace the bilinear parameters in terms of the VEVs are given below: Now let us demonstrate the diagonalization of the mass matrices in different sectors following the same prescription as in [17] but in the presence of the soft terms.

A.1 CP-odd scalar sector
The mass term for the pseudoscalar sector can be extracted from the scalar potential as, where M 2 P is the 3 × 3 mass matrix that can be block diagonalized as follows: The elements of B 2 P are given by, The matrix B 2 P can be fully diagonalized by an orthogonal transformation where, O γ1 is given in Eq. (6), which entails the following relations Using Eq. (31), Eq. (33) can be inverted to solve for λ 10 , λ 11 and λ 12 as

A.2 Charged scalar sector
Similar to the pseudoscalar case, the 3 × 3 charged sector mass matrix M 2 C can also be block diagonalized as: where, Further, the charged scalar mass matrix can be completely diagonalized with the use of the rotation matrix O γ2 (given in Eq. (6)) as Thus, we will have the following relations: cos γ 2 sin γ 2 (m 2 These equations in conjunction with Eq. (35) will enable us to solve for λ 7 , λ 8 , and λ 9 as given below: where, the other three couplings (λ 10 , λ 11 & λ 12 ) can be replaced using Eq. (34).

A.3 CP-even scalar sector
The mass terms in the neutral scalar sector can be extracted from the potential as, where, M 2 S is the 3 × 3 symmetric mass matrix whose elements are given by, This mass matrix should be diagonalized via the following orthogonal transformation where, O α has already been defined in Eq. (9). Inverting the above Eq. (40), we get, which enables us to solve for the remaining six lambdas as follows: NFC 2HDM [52,53] for a scenario with two different H + . Following Ref. [54], the branching ratio for b → sγ is controlled by the C eff 7L and C eff 7R Wilson coefficients: where the normalization by Br (b → ceν) helps canceling some of the hadronic uncertainties. The effective Wilson coefficients read where, as in the usual analysis of 2HDMs [54], the leading log QCD corrections in the SM are described by and η = α s (M Z )/α s (µ), where µ is the QCD renormalization scale, µ ≈ 221 MeV. Taking into account the absence of tree-level FCNCs, the coefficients in Eqs. (44a) and (44b) can be recast as where the A + terms correspond to our NP (charged-Higgs) contributions. These contributions can be further broken down into with y q = m 2 q /M 2 H + and in which X and Y are the charged-Higgs coupling to left-and right-handed quarks, respectively, and the loop functions are given by where p F q (a, b; c; d) is the Hypergeometric Function. Finally, the SM amplitude is given by (keeping only the top contribution) where x t = m 2 t /M 2 W . So far, we have presented the analysis of the b → sγ processes in a 2HDM where FCNCs are absent. To extend these results to our model, we redefine Eqs. (46) and (47) to account for both charged-Higgs contributions: where now y i q = m 2 q /M 2 H + i , and where we can see the X i and Y i couplings now carry an index, denoting the H + 1 and H + 2 chiral (P L and P R ) couplings to quarks. In the present model, these couplings can be extracted from Eqs. (18a) and (18b): Y 1 = − tan β 2 cot β 1 cos γ 2 sin β 2 + sin γ 2 , (52b) We now have all the relevant information needed to compute the b → sγ branching ratio in our model. As advertised, the only dependencies on the BSM degrees of freedom is through tan β 1 , tan β 2 , γ 2 , which control the couplings, the charged-Higgs masses, m H + 1 and m H + 2 , which will affect the loop functions. Finally, the SM prediction for the b → sγ branching ratio can be found in ref [55], and the experimental values in [19].

B.2 Neutral Meson Mixing: ∆M Bq
A very restrictive aspect of BSM models comes from neutral meson oscillations. These processes, for models without tree-level FCNCs, are forbidden at tree-level, but may have sizable one-loop contributions. The left panel in Fig. 7 represents the SM contribution for such processes, whereas the other two diagrams represent the additional contributions in our Z3HDM scenario. To obtain some qualitative intuitions we write the effective a,b=u,c,t i,j=H ± 1 ,H ± 2 λ a λ b ω a ω b S(y a , y b ) 4 + X ia X ib I 1 (y a , y b , y i ) + X ja X jb I 2 (y a , y b , y i , y j ) O F .
The SM contribution is encoded in S(y a , y b ), normalized by a factor of 4 to account for the summation on the charged Higgs. The I 1 (y a , y b , y i ) contributions are due to the mixed W ± − H ± i boxes, and I 2 (y a , y b , y i , y j ) are the H ± i − H ± j boxes in Fig. 7. The above expression is valid in the zero external momenta approximation, where the down-type quark masses are taken to be zero. We use X ia to denote the coupling between the charged-Higgs H ± i and the up-quark a, which, as seen in Eq. (52a), are flavor universal, i.e., X 1a = X 1 = − cot β 2 sin γ 2 and X 2a = X 2 = cot β 2 cos γ 2 for H ± 1 and H ± 2 , respectively. The quantities y a and y i stand for the ratios m 2 a /M 2 W and m 2 H + i /M 2 W respectively. The specificities of the neutral meson under consideration are contained in the CKM elements λ a , and the dimension-6 operators O F . For a generic meson P = (q 1 , q 2 ), these are defined as Finally, the loop functions are given by g(x, y, z) = x(x − 4) ln x (x − 1)(x − y)(x − z) , I 1 (y a , y b , y i ) = g(y a , y b , y i ) + g(y b , y i , y a ) + g(y i , y a , y b ), (55b) h(x, y, w, z) = I 2 (y a , y b , y i , y j ) = h(y a , y b , y i , y j ) + h(y b , y a , y i , y j ) + h(y i , y a , y b , y j ) + h(y j , y a , y b , y i ).
The limiting cases where, for instance, the same Higgs runs in the I 2 box diagram should be carefully dealt with, as the loop functions are only apparently divergent for x i = x j , but indeed have a well-defined limit.
Finally, we can obtain ∆M P from the effective Lagrangian, where M P is the meson mass, f P its decay constant, and B P is its bag parameter. The 2HDM limit (with no tree-level FCNCs) of Eq. (53) can be easily extracted, taking some care on the symmetry factors. As in the b → sγ computations, it would be possible to parametrize these results to match numerical results with higher-order corrections. The experimental values which will determine the experimentally allowed region are taken from [19], whereas the relevant hadronic parameters can be found in [56].