Current bounds on the Type-Z $Z_3$ three Higgs doublet model

Type-Z models, where charged leptons, up type quarks, and down type quarks each couple to a different scalar, are only possible when there are three or more Higgs doublets. We consider the Type-Z three Higgs doublet model imposed by a softly broken $Z_3$ symmetry. We take into account all theoretical and experimental constraints, highlighting the importance of perturbative unitarity and bounded from below conditions that we develop here, and which have a strong impact on $B \rightarrow X_s \gamma$. Since there can be cancellations between the two charged Higgs in $B \rightarrow X_s \gamma$ (and in $h \rightarrow \gamma \gamma$), the lower bounds obtained on the charged Higgs masses are alleviated. We also discuss in detail the important physical differences between exact alignment and approximate alignment, and present some useful benchmark points.

After the observation in 2012 by ATLAS and CMS [1, 2] of a new scalar particle closely resembling the Standard Model (SM) Higgs boson, the search for physics beyond the Standard Model (BSM) is now the main goal of the LHC experiments. Popular extensions where only Higgs doublets are added to the SM have been extensively studied and allow for both the agreement with the experimental results and the possibility of new features; for reviews see [3][4][5].
The simplest extension, the two-Higgs-doublet model (2HDM), can provide new sources of CP-violation necessary to fulfill the Sakharov criteria for baryogenesis [6]. However, the most general Higgs-fermion Yukawa couplings generically yield Higgs-mediated flavorchanging neutral "currents" (FCNCs) at tree level, in conflict with experimental observations. A common method to have FCNCs sufficiently suppressed is to impose symmetries on the Lagrangian: tree-level FCNC effects can be completely removed by establishing how the fermion and scalar fields have to transform under the chosen symmetry. In the two Higgs doublet model (2HDM) this can be achieved by imposing a Z 2 symmetry [7,8]. Reference [9] showed that in general N Higgs doublet models (NHDM) the Yukawa coupling matrices to fermions of a given electric charge remain proportional (thus removing FCNCs) under the renormalization group running if and only if there is a basis for the Higgs doublets in which all the fermions of a given electric charge couple to only one Higgs doublet. The models are then classified based on these choices. The four (five) distinct types of Yukawa couplings in models with two (more than two) doublets that fit this requirement were introduced in [9] and denoted in [10] by Types I, II, X (also known as lepton-specific), Y (flipped), and Z, according to Type-Y: φ u = φ e = φ d , where φ u,d,e are the single scalar fields that couple exclusively to the up type quarks, down type quarks, and charged leptons, respectively. In this work, we set our attention on the Type-Z that can only appear for NHDM with N > 2. It is interesting to see what differences there are in this new type of model, since it decouples completely the up quark, down quark and charged lepton sectors from one-another. There are have been implementations of Type-Z in three-Higgs-doublet models (3HDM) using a Z 2 × Z 2 symmetry [11,12] or Z 3 [13,14]. For this work, we choose to use a Z 3 symmetric potential. This symmetry is realizable through the following representation, the theoretical constraints coming from perturbative unitarity, discussed explicitly for the Z 3 3HDM model in Ref. [16] and BFB conditions, which we develop here. Compatibility with the bounds coming from Higgs searches is also checked with the newest version of the HiggsBounds-5. 9.1 (HB5) code [17]. In particular, we show that recent LHC measurements exclude all points in Fig. 2 of Ref. [15], for the same parameter choices. We then show that by scanning for a larger range of parameters (away from exact alignment, but still consistent with all experimental data) we can obtain viable points corresponding to smaller masses for the additional particles.
In Sec. I we describe succinctly the scalar and Yukawa sectors of the Z 3 3HDM model, discussed also in [13][14][15]. The theoretical and experimental constraints are described in Sec. II. In Sec. III we describe the impact of current LHC measurements on the 125GeV scalar decays, both excluding and including the impact of HB5 bounds. In particular, we discuss the fact that the couplings of the two charged scalars may have different signs, thus allowing for canceling contributions to h → γγ. A similar effect is possible in B → X s γ, thus alleviating the lower bounds on charged scalar masses. This is discussed in Sec. IV and Sec. V, where we explore the regions of parameters allowed by the different constraints imposed, starting from the experimental limits on the BR(B → X s γ) and progressively varying the ranges on our parameter scans. Our work highlights the importance of going beyond strict alignment, when procuring the full range of possibilities existent within the Z 3 3HDM. We present illustrative benchmark points in Sec. VII and discuss our conclusions in Sec. VIII, leaving the appendix for the full expression of some couplings required in our calculations.

A. Scalar sector
Taking the potential defined by [13], the terms invariant under the chosen transformation, with the quartic part and the quadratic part also including terms, m 2 12 , m 2 13 and m 2 23 , that break the symmetry softly. After spontaneous symmetry breaking (SSB), the three doublets can be parameterized in terms of its component fields as: where v i / √ 2 corresponds to the vacuum expectation value (vev) for the neutral component of φ i . It is assumed that the scalar sector of the model explicitly and spontaneously conserves CP. 2 That is, all the parameters in the scalar potential are real and the vevs v 1 , v 2 , v 3 , are also real. With this assumption, the scalar potential of Eq. (3) contains eighteen parameters. The vevs can be parameterized as follows: leading to the Higgs basis [19][20][21] to be obtained by the following rotation, The scalar kinetic Lagrangian is written as and contains the terms relevant to the propagators and trilinear couplings of the scalars and gauge bosons.
We can now define orthogonal matrices which diagonalize the squared-mass matrices present in the CP-even scalar, CP-odd scalar and charged scalar sectors. These are the transformations that take us to the physical basis, with states possessing well-defined masses. Following Ref. [13,14], the twelve quartic couplings can be exchanged for seven physical masses (three CP-even scalars, two CP-odd scalars and two pairs of charged scalars) and five mixing angles. The mass terms in the neutral scalar sector can be extracted through the following rotation, 1 Notice that we use x i in place of Ref. [13]'s h i , because for us h i are the physical neutral scalar mass eigenstates. 2 Strictly speaking, it is not advisable to assume a real scalar sector while allowing the Yukawa couplings to carry the phase necessary for the CKM matrix. This is also a problem with the so-called real 2HDM [18]. One can take the view that the complex terms and their counterterms in the scalar sector exist, with the former set to zero.
where we take h 1 ≡ h 125 to the be the 125GeV Higgs particle found at LHC. The form where For the CP-odd scalar sector, the physical basis is chosen as G 0 A 1 A 2 T and the transformation to be where O γ 1 is defined in order to diagonalize the 2x2 submatrix that remains in the Higgs basis, with the form For later use, we define the matrix P as the combination For the charged scalar sector, the physical basis is G + H + 1 H + 2 T and the transforma- where We write the masses of H + 1 and H + 2 as m H ± 1 and m H ± 2 , respectively. The matrix Q is then defined as the combination Considering that the states in the physical basis have well-defined masses, we can obtain relations between the set and the parameters of the potential in Eq. (3), as shown in Ref. [13,14]. We performed an extensive scan of the parameter space in Eq. (19). Our fixed inputs are v = 246 GeV and m h1 = 125 GeV. We then took random values in the ranges: These parameter ranges will be used in all scans and figures presented below, except where noted otherwise. The lower limits chosen for the masses satisfy the constraints listed in Ref. [22]. 3

B. Higgs-Fermion Yukawa interactions
One can now impose the Type-Z model on the Yukawa Lagrangian, by establishing how the fields behave under the Z 3 transformation. For this, there are multiple possibilities that differ on which of the scalars gives mass to each type of fermion. We follow the choice made by Das and Saha [13]. The scalar doublets φ 1 and φ 2 transform nontrivially as: where ω = e 2π i/3 . For the fermionic fields, we consider that under Z 3 while the rest of the fields remain unaffected. It follows that the Yukawa coupling matrices are now restricted: φ 1 only has interaction terms with the charged leptons, giving them mass; φ 3 and φ 2 are responsible for masses of the up and down type quarks, respectively. When taking into account the restrictions imposed by the symmetry, the Yukawa couplings to fermions can be written in a compact form. For the couplings of neutral Higgs to fermions, where we group the physical Higgs fields in a vector, as The coefficients are given in Eq. (25), where we introducev i = v i /v, with the vevs in Eq. (7). Note how the coupling of each type of fermion depends on entries of the diagonalization matrices in Eqs. (11) and (15). The couplings of the charged Higgs, H † 1 and H † 2 , to fermions can be expressed as while for leptons, V ij = δ ij since we are considering massless neutrinos. The couplings are for leptons and quarks, respectively.

II. CONSTRAINTS ON THE PARAMETER SPACE
In this section we study the constraints that must be applied to the model parameters in order to ensure consistency.

A. Theoretical Constraints 1
We impose perturbativity unitarity, sufficient bounded from below (BFB) conditions, and the oblique parameters S, T , and U .

BFB conditions on the 3HDM
As basic requirements for any physical theory, the Higgs potential must satisfy conditions that ensure it possesses a stable minimum, around which one can perform perturbative calculations. That is, it must be bounded from below, meaning that there is no direction in field space along which the value of the potential tends to minus infinity. This need of a non-trivial minimum is then translated to conditions on the parameters of the potential.
Focusing on the study of the 3HDM constrained by a Z 3 symmetry, the quartic terms in Eq. (4) can be written as where V 0 has the terms in λ 1→9 and V 1 the terms λ 10→12 . If the potential were just V 0 in Eq. (28), then the BFB necessary and sufficient conditions would be simply those given by Klimenko in Ref. [23]. The problem, not yet solved for the 3HDM with a Z 3 symmetry is the V 1 part. We will introduce sufficient conditions for BFB by bounding the potential by a lower potential. To do that we follow [23,24], checking for neutral minima. Neutral directions in the Higgs space correspond to situations when all φ i are proportional to each other 4 . Along these directions, we can then define It then follows that for V 0 , and for V 1 , where δ i are some combination of the phases θ i . Considering that x, y, z > 0 by definition, we can start our strategy of bounding the potential by a lower one with Notice that for non-negative x, y, z one has Therefore, 4 Other directions, along which the strict proportionality of all three doublets does not hold, are called charge-breaking (CB) directions. In recent works [25,26], it has been proven that these directions can lead to pathological situations for other symmetries in the 3HDM. It is then required to consider these directions when doing a complete work of looking for necessary and sufficient BFB conditions. Our contribution to the analysis of the Z 3 symmetry is to specify sufficient conditions along the neutral direction. and combining Eq. (34) with Eq. (30), it follows that where with the definitions, Now, for the potential V BFB the necessary and sufficient conditions are obtained from Ref. [23]: where As V 0 + V 1 > V BFB , these conditions are sufficient conditions for the original potential. They are not necessary, and therefore might be throwing away part of the parameter space. However, it still gives us a very good sense of the possibilities within the Type-Z 3HDM.

Unitarity
In order to determine the tree-level unitarity constraints, we use the algorithm presented in [16]. As described there, we have to impose that the eigenvalues of the scattering Smatrix of two scalars into two scalars have an upper bound (the unitarity limit). As these arise exclusively from the quartic part of the potential, the eigenvalues obtained for a Z 3 symmetric potential in Section 4.4 of [16] can also be used for the potential with quadratic soft-breaking terms, Eq. (3). The conversion between the notation of the algorithm and the potential chosen, Eq. (4), is as follows: Denoting by Λ i the eigenvalues of the relevant scattering matrices, we have 21 Λ's to calculate for each set of physical parameters randomly generated, and the condition to impose is that

Oblique parameters ST U
In order to discuss the effect of the S, T, U parameters, we use the results in [27]. To apply the relevant expressions, we write the matrices U and V used in [27] with the notation choices that we made when obtaining the mass eigenstates in section I A. We start with the and find, by comparison with Eqs. (10) and (13), that V is The 3 × 3 matrix U defined as gives us the correspondence U = Q T from Eq. (16).
Having applied the expressions for S, T, U , the constraints implemented on S and T follow Fig. 4 of Ref. [28], at 95% confidence level. For U , we fix the allowed interval to be

B. Theoretical Constraints 2
As we want to explore the range of low tan β 1 and tan β 2 we should avoid that the Yukawa couplings become non-perturbative. We have, in our model We require We see from Fig. 1 of Ref. [15] that the constraints coming from ∆M b,s tend to exclude very low values on tan β. Thus, we take

D. LHC Constraints
For comparison with experiment, we consider only the contributions of the lowest nonvanishing order in perturbation theory. The decays that require one-loop calculations are those of neutral scalars into two photons (h j → γγ), one Z and one photon (h j → Zγ), and two gluons (h j → gg). The final formulas for the first two widths are given in Ref. [29], only having to adapt the particles and their couplings to our case. The formula for the width h j → γγ reads, where, noticing that for scalars the Y terms in [29] vanish, We used where m f (m H ± k ) is the mass of the relevant particle in the loop, while m h j is the mass of the decaying Higgs boson. The function f (τ ) is defined in the Higgs Hunter's Guide [3], and the couplings C j and λ h j H + k H − k for this model are written in the appendix. They were derived with the help of the software FeynMaster [30,31], that uses QGRAF [32], FeynRules [33,34] and FeynCalc [35,36] in an integrated way.
The decay into gluons can be obtained from the expression for the γγ decay, where and the sum runs only over quarks q. For the 125GeV scalar, the coupling modifiers, are calculated directly from the random angles generated and constrained to be within 2σ of the most recent ATLAS fit results, [37, Table 10]. Having chosen a specific production and decay channel, the collider event rates can be conveniently described by the cross section ratios µ h if , Starting from the collision of two protons, the relevant production mechanisms include: gluon fusion (ggH), vector boson fusion (VBF), associated production with a vector boson (VH, V = W or Z), and associated production with a pair of top quarks (ttH). The SM cross section for the gluon fusion process is calculated using HIGLU [38], and for the other production mechanisms we use the results of Ref. [39]. Each of the 3HDM processes is obtained by rescaling the SM cross sections by the relevant relative couplings. As for the decay channels, we calculated the branching rations for final states f = W W, Z Z, b b, γ γ and τ + τ − . Finally, we require that the µ h if for each individual initial state × final state combination is consistent, within twice the total uncertainty, with the best-fit results presented in the most recent study of data collected at √ s = 13 TeV with the ATLAS experiment [37, Figure 5].
For the heavier neutral and charged scalars, we use HiggsBounds-5.9.1 in Ref. [17], where a list of all the relevant experimental analyses can be found. For the decays allowing for off-shell bosons, we use the method explained in [40]. We also consider the constraints coming from b → sγ, as we explain in sections IV and V.

III. DECAYS OF h 125 IN THE Z 3 3HDM
In this section, we use the scan ranges defined in Eq. (21), pass them through all theoretical and experimental constraints, and we study the impact on the decays of the 125GeV Higgs h 1 = h 125 found at LHC.
The contribution from the two charged scalars to the h 125 → γγ decay process is shown in Fig. 1. There are two interesting regimes. To the left (right) of the vertical line at coordinate zero, the two charged Higgs conspire to decrease (increase) the branching ratio into γγ. Most of the points are on the left and correspond to a significant reduction of the decay width. However, there are indeed points on the right, which allow for an increase which could be up by 20%. We have also confirmed the existence of allowed results where the destructive interference between the two charged Higgs leads to a null X H , occurring when the signs of the couplings λ h j H +  The points of Fig. 1 where |X H | 2 is large, for which the charged Higgs provide a considerable contribution to the overall h 125 → γγ decay rate (the latter, still within current bounds) is only obtained for very fine tuned points in parameter space with some charged Higgs mass below 200GeV. As we will see in Figs. 12-13 below, this is a very constrained (fine tuned) region.
The set of points that are consistent with all the bounds is now plotted in the sin (α 2 − β 2 ) − sin (α 1 − β 1 ) plane as shown in Fig. 2. Comparing with the plot in the same plane shown in [13, Fig.1], it can be seen that the use of more recent experimental data for the simulated results leads to us being closer to the alignment limit, defined by α 1 = β 1 and α 2 = β 2 .
However, as we will illustrate below, points in parameter space close to the alignment limit exhibit physical properties which differ significantly from the exact alignment limit. Figure 2: Results of the simulation in the sin (α 2 − β 2 ) − sin (α 1 − β 1 ) plane. The green points passed all constraints including HB5, while the red points did not pass HB5.
To study the allowed regions for the cross section ratios µ h if , we follow [29,41] and calculate each µ h if using all production channels. Our set of points is then shown in Figs. 3 -6. Similar to the complex 2HDM analyzed by Fontes, Romão and Silva in [29], there is a strong correlation between µ Zγ and µ γγ in our Type-Z model, as shown in Fig. 6.    Such a correlation is also visible between µ ZZ and µ γγ in Fig. 3. It is less apparent in correlations with τ + τ − and bb, as shown in Figs. 4 and 5.

A. Introduction
It is well known that the experimental bounds on B → X s γ place stringent restrictions on the parameter space of models with charged scalars [11,[42][43][44][45]. Most notably, there is a bound on the mass of the only charged Higgs boson present in the Type-II 2HDM which, at 95% CL (2σ), is according to [44] m H + > 580 GeV .
The exact value for this bound depends on both the theoretical approximations [46] and the experimental errors. The experimental average gives [47] BR while the NNLO calculation within the SM yields [11,48] with an error of about 5%. As explained below, we will take an error of 2.5% around the central value of the calculation and, following [11], we consider 99% CL (3σ) for the experimental error: 2.87 × 10 −4 < BR(B → X s γ) < 3.77 × 10 −4 . (66)

B. The calculation
We follow closely the calculation by Borzumati and Greub in Ref. [42]. There, the new contributions from the charged Higgs bosons are encoded in the Wilson coefficients, where we are using the notation in Ref. [42] which should be consulted for the definitions and also for the procedure used in evolving the coefficients to the scale µ b = m b . The dependence on the charged Higgs mass appears because the functions C 0,eff i,YY , C 0,eff i,XY , C 1,eff i,YY , and C 1,eff i,XY depend on y = m 2 t /m 2 H + , while the SM coefficients depend on x = m 2 t /M 2 W . For models with multiple charged Higgs there is one contribution (and one parameter y k ) for each particle. A model with two charged Higgs is discussed in [11,12], with interesting earlier work highlighting the possible cancellation between the two charged Higgs contributions appearing in Refs. [49,50]. We obtain, for example, where we wrote explicitly the dependence on the charged Higgs masses, and used We took the input parameters from Ref.
We find that much of the parameter space considered in Ref. [15] is forbidden. This is most apparent by considering their Fig. 2, which we turn to next.

B. The effect of other constraints
From the previous plots, the conclusion that we can have one of the charged Higgs relatively light if the other is sufficiently heavy seems correct. However we now show that for this choice of parameters this is not the case. With the choice of Eqs. (72)-(73, the bounds from the decays of the 125 GeV Higgs are simply satisfied. However the same is not true for current bounds on heavier scalars. Indeed, every single point in Fig. 7 is excluded by HB5; not a single point remains. This will be explained in detail in the following section.

C. Enlarging the region of good points
We discovered that the situation described in the previous section is a consequence of the small range chosen for γ 2 . To illustrate this, we kept the other conditions in Eqs. (72)-(73), but allowed for and (for Fig. 8) also varied tan β 1 . The points which survive HiggsBounds-5.9.1 are shown in dark green on the left panel of Fig. 8. The allowed points for tan β 1 = 10 are concentrated around γ 2 = 0, ±π/2, excluding γ 2 = π/6, π/4, π/3. Taking the interval in Eq. (74) one can indeed find regions of good points. 6 It is interesting to compare with what happens with the previous version of HiggsBounds-5.7.1, shown on the right panel of Fig. 8. For that case there are many points allowed for all values of γ 2 , even for tan β 1 = 10. We have found that this is due to the recent bounds on h 2,3 → τ + τ − decay in Ref. [52], included in HiggsBounds-5.9.1 but not in HiggsBounds-5.7.1, which used the previous bounds [53,54]. 7 To better illuminate this point, we show σ(pp → h 2 ) × BR(h 2 → τ τ ) versus m h 2 in Fig. 9. In this figure, the parameters are as in Eqs. (72)-(73, except that γ 2 ∈ [−π/2, π/2]. Points in cyan are points that pass all constraints before HiggsBounds. In light green are the points in the restricted interval γ 2 ∈ [π/6, π/3]. In the left panel points in dark green are those who survided after HiggsBounds-5.7.1. In the right panel we have the same situation but now we used HiggsBounds-5.9.1. We see that there were good points in the restricted interval γ 2 ∈ [π/6, π/3] in the left panel, but they disappeared with the newer version HiggsBounds-5.9.1. We have confirmed that similar plots can be obtained for h 3 . This is a good point to stress again the role that the LHC is having in constraining models with new scalar physics. One sees the strong impact that the updated LHC results have in constraining the Z 3 3HDM. This highlights the importance that the new LHC run will have in constraining the parameter space of extended scalar sectors. To better understand the behaviour of σ(pp → h i ) × BR(h i → τ τ ) (i = 2, 3) we can make the simplified assumption 8 that this product is proportional to 7 In Ref. [15] the strong constraints from neutral scalar decays into τ τ still seemed to allow points with the choices in Eqs. (72)-(73). 8 We are neglecting the dependence of the cross section on the mass.
where we are assuming that the production occurs mainly via gluon fusion with the top quark in the loop. Now, using the assumptions of Eq. (73) in Eq. (25), we have where, for Fig. 9, β 1 , β 2 are fixed and γ 2 ∈ [−π/2, π/2]. Fig. 10 shows the functions in Eq. (75)f 2 for h 2 and f 3 for h 3 -for tan β 1 = 10 and tan β 2 = 2 as in Eq. (72), but keeping γ 2 free. We see that these functions are largest precisely in the approximate interval ±γ 2 ∈ [π/6, π/3]. This explains why these points are the first to be excluded by the bounds on σ(pp → h 2 ) × BR(h 2 → τ τ ), and why, going outside such bounds, some points can be preserved. 9

D. The effect of tan β's
In the last section we saw that while maintaining the main features of Eqs. (72)-(73, but enlarging the range of variation of γ 2 , we could find points allowed by all current experimental constraints. Here we exploit the variation of both tan β's in the range subject to the condition of perturbativity of the Yukawa couplings in Eq. (52). The result is shown in Fig. 11. We see that by varying the range of tan β's we can have smaller masses for the charged Higgs bosons. For tan β < 1 it is even possible to have both charged Higgs with masses below 400 GeV.

VI. GOING BEYOND EXACT ALIGNMENT
We have performed a completely uniform scan and found out that very few points survived and those were not too far away from the alignment condition of Eq. (73). So another strategy can be to scan points that differ from the perfect alignment of Eq. (73) by 1% or 10%.
In Fig. 12 we show the results for the case when we allow the parameters to differ 1% from the perfect alignment limit. Next we considered the case when the difference for perfect alignment was 10%. This is shown in Fig. 13. We see that the acceptable points which differ more from perfect alignment are less frequent, as expected. 10 Nevertheless, one can still find many points which differ from exact alignment by as much as 10%, while satisfying all experimental and theoretical constraints. And such points do allow for qualitatively different predictions, as we saw when looking at the charged scalar masses consistent with b → sγ. We conclude that imposing perfect alignment is too constraining and does not cover all the interesting features of the Z 3 3HDM.

A. Unusual signals of charged scalars
As we have seen, the contributions of the two charged scalars can exhibit large cancellations in the decays h → γγ and B → X s γ. 11 For some choices of parameter space, it is even possible that there are cancellations in both decays simultaneously. This is illustrated in Fig. 14. Figure 14: Points with significant approximate cancellation in both h → γγ (horizontal axis) and B → X s γ (vertical), which pass all theoretical and experimental bounds, including HB5. Color code: cyan is perfect alignment, red means alignment within 1%, and blue means alignment within 10%. The blue box guides the eye to those points closest to (0,0). Such charged scalars would, thus, be difficult to probe indirectly.
Notice that points with exact alignment, in cyan in Fig. 14, do not allow for cancellation in h → γγ; but alignment with 1% already does.
Most points within the blue box close to (0,0) have H + 2 decays into quarks or leptons, which are being sought at LHC. But there are points which could also be difficult to probe directly with such common searches, even tough one or both charged scalars might have relatively small masses. Indeed, one can find fine-tuned points in parameter space where the H + 2 does not decay primordially into quarks or leptons, but rather as H + 2 → H + 1 h j with h j = h 1 , h 2 , A 1 . We propose that such decays be actively searched for at LHC's next run. To aid in that experimental endeavour, we present some benchmark points (BP) in the next section.

VII. ILLUSTRATIVE BENCHMARK POINTS
This section is devoted to some benchmark points/lines, with features which may prove useful for the experimental searches.
There has been a recent interest in the literature for unusual decays of the charged Higgs [56], specially those in which the charged Higgs decays to W + h i where h i is any of the scalars or pseudo scalars in the model.
We have performed a search in our large data sets and found many points where BR(H + 1 → W + + h 125 ) was larger than 80%. From those we selected three benchmark points (BP) that we list in table I. For each of these BP we let the mass of the H + 1 vary, leaving all the other parameters fixed, obtaining benchmark lines. All these points verify all the constraints, including those from HiggsBounds-5.9.1. These BP all have the characteristic that the dominant decay of the charged H + 1 is not in the tb channel, but in W + h 125 , which makes these interesting and deserving to be searched at the LHC. Notice that, for  > M W +m A 1 explaining the decrease in our prefered branching ratio (see Fig. ??). The same happens for the channel H + 1 → W + h 2 for BP3 as can be seen in Fig. 16.

VIII. CONCLUSIONS
Multi Higgs models with N ≥ 3 allow for the possibility that all fermions of a given charge couple exclusively to one dedicated scalar. These are known as Type-Z models, and constitute a fifth alternative beyond the four natural flavour conservation models allowed in the 2HDM. We investigate the current bounds on the Type-Z 3HDM imposed by a Z 3 symmetry. We perform an up-to-date analysis including the latest data for the 125GeV Higgs [37], bounds on new scalars through the HiggsBounds-5.9.1 code [17], and the very  table I. important theoretical constraints. We use the theoretical bounds from unitarity [16] and BFB; the latter developed here for the first time. We stress the importance of using the most recent LHC bounds, which constrain severely the allowed parameter space. In particular, we show that bounds from h 2 → τ + τ − alter significantly some results in the literature [15]. This is clearly visible in our Fig. 8 and Fig. 9. Moreover, we also stress the fact that interesting physical observables may differ significantly when one considers situations close to the alignment limit, versus adopting the exact alignment limit. Indeed, current LHC bounds on the productions and branching ratios of the 125GeV neutral scalar force the measured couplings to lie close to those obtained for the SM Higgs. Nevertheless, forcing those couplings to match exactly those in the SM is too constraining on the parameter space and precludes much of the interesting new features that the Z 3 3HDM has.
We look at the constraints allowed by current data on the 125GeV Higgs decays, including a detailed look at h → γγ and its correlations with the other decays. We point out the possibility that the contributions from the two charged scalars might cancel in h → γγ. This is also possible in B → X s γ, and we explore explicitly how this allows for lower masses for the charged scalars. We provide illustrative benchmark points to aid in experimental searches. By comparing the constraints from HiggsBounds-5.7.1 and the newer HiggsBounds-5.9.1 we highlight the importance that the next LHC run will have in further constraining this model, or perhaps, finally uncovering new physics in the scalar sector.

Appendix A: Some important couplings
This appendix is devoted to some important couplings for the Z 3 3HDM used in our calculations. In our conventions these couplings include the i from the Feynman rules. These couplings were derived with the help of the software FeynMaster [30,31].

Scalar couplings to W ± bosons
We find for the neutral scalar couplings to W + W − , is to be used in Eq. (56).

Scalar couplings to charged Higgs
The couplings of the scalars h j with j = 1, 2, 3 to the charged Higgs H ∓ k−1 , H ± l−1 where k, l = 2, 3 (we do not consider here the charged Goldstone) are, where we have defined k ≡ k − 1, l ≡ l − 1 with j = 1, 2, 3 and k, l = 2, 3. Recall that v k = v k /v. The coupling λ h j ,H + k ,H − l is to be used in Eq. (57).