Heavy charged scalars from $c\bar{s}$ fusion: A generic search strategy applied to a 3HDM with $\mathrm{U}(1) \times \mathrm{U}(1)$ family symmetry

We describe a class of three Higgs doublet models (3HDMs) with a softly broken $\mathrm{U}(1) \times \mathrm{U}(1)$ family symmetry that enforces a Cabibbo-like quark mixing while forbidding tree-level flavour changing neutral currents. The hierarchy in the observed quark masses is partly explained by a softer hierarchy in the vacuum expectation values of the three Higgs doublets. As a consequence, the physical scalar spectrum contains a Standard Model (SM) like Higgs boson $h_{125}$ while exotic scalars couple the strongest to the second quark family, leading to rather unconventional discovery channels that could be probed at the Large Hadron Collider. In particular, we describe a search strategy for the lightest charged Higgs boson $H^{\pm}$, through the process $c\bar s\to H^{+}\to W^+\,h_{125}$, using a multivariate analysis that leads to an excellent discriminatory power against the SM background. Although the analysis is applied to the proposed class of 3HDMs, we employ a model-independent formulation such that it can be applied to any other model with the same discovery channel.


I. INTRODUCTION
The Standard Model (SM) remarkably stands as one of the most successful theories in physics. However, it can still be considered rather ad hoc in its nature, with unexplained features that arise from fitting the experimental data. In addition, it fails to offer an explanation to several observed natural phenomena such as dark matter, neutrino masses or baryon asymmetry in the universe. It is then natural to study extensions of the SM that, while retaining its predictive power, offer explanations or shed light into the origin of e.g. the hierarchy of fermion masses or rather specific flavour structure of the SM.
There is a plethora of such beyond the SM (BSM) theories, but not many of those offer unconventional features testable at the current experiments.
One of the simplest and most studied extensions is the class of the so-called Two-Higgs Doublet Models (2HDMs) that add a second SU(2) L doublet to the SM (an extensive review can be found in Ref. [1]). The 2HDMs offer interesting phenomenological signatures and can lead to e.g. extra sources of CP violation, dark matter candidates and stable vacua at high energies. However, they typically introduce many new free parameters, fail to address the origin of the mass hierarchy in the fermion sector of the SM and require extra discrete symmetries to avoid tree-level Flavour Changing Neutral Currents (FCNCs).
While retaining most of the features of 2HDMs, 3HDMs can offer explanations to yet unexplained features of the SM with predictions testable in the current collider measurements. In particular, the increased field content makes it possible to impose higher symmetries, which in turn can lead to interesting flavour structures.
As shown in Refs. [12,13], the most constraining realisable abelian symmetry of the scalar potential in 3HDM is U(1) × U (1). In this work, we promote the U(1) × U(1) symmetry of the scalar sector to the fermion sector, hereinafter called U(1) X × U(1) Z , in such a way that (1) no tree-level FCNCs are present, (2) a Cabibbo-like mixing is enforced, and (3) the fermion mass hierarchies are related to a hierarchy in the three vacuum expectation values (VEVs) of the doublets. This leads to a model that, although remarkably simple due to its high symmetry, is still capable of both reproducing the current experimental data and providing the exotic collider signatures. The latter is due to the fact that, as a consequence of the model symmetries, the new scalar states (both charged and neutral) couple dominantly to the second quark family.
At the LHC, the searches for charged Higgs bosons are generally categorized into two mass regions depending on whether its mass m H ± is smaller or bigger than the top quark mass m t . The motivation of this categorization comes from the properties of H ± within the various 2HDM types or supersymmetric models. Usually, for a heavy charged Higgs state (m H ± m t ), the dominant production and decay channels in the LHC context are pp → H − tb (H +t b) and H + → tb (H − →tb), respectively [14,15]. Apart from this channel, production of H ± in vector boson (W ± Z) fusion followed by the H ± → W ± Z decay is prominent in Higgs triplet models such as the Georgi-Machacek model [16].
This channel has also been searched for by the ATLAS [17] and CMS [18] collaborations recently. Conversely, a light charged Higgs boson (m H ± m t ) that decays to τν [15,19], cs [20,21] or cb [22] has also been searched for at the LHC. Previously, at LEP, pair production of H ± was considered where H ± subsequently decays to a W ± A pair [23,24] (where A is a scalar with mass m A > 12 GeV and predominantly decays to bb pairs).
Searches for heavy H ± become increasingly important with the rise of the LHC centerof-mass energy and luminosity, thus it is important to explore new production and decay modes of H ± that are predicted by various BSM theories. In this paper, we particularly focus on a new search channel where a heavy H ± resonantly decays to a W ± h 125 pair after being produced in cs (cs) fusion. This rather uncommon search channel leads to testable predictions of our model at current LHC energies. In Refs. [25][26][27], the H ± → W ± h 125 decay is considered where the H − (H + ) is produced in association with a tb (tb) pair. In our case, H + is produced in the s-channel resonance through cs fusion. In Ref. [28], the possibility of sizable cs → H + production cross section is discussed in the SUSY context where a squark mixing can circumvent the chiral suppression of the single H ± production.
In Ref. [29], cs → H + production is shown to be dominant for a 2HDM with a Yukawa sector chosen such that one doublet couples strongest to the second generation. In our model, we will see that the csH + andcsH − couplings are sizable due to the hierarchy in VEVs combined with the particular structure of the Yukawa sector.
In section II, we introduce the model, its fermion and scalar sectors and their interplay as given by the U(1) X × U(1) Z symmetry, the VEV hierarchy and the spectrum of the theory. In section III, we discuss the charged Higgs boson production and decay channels of the model and introduce a model-independent way to study H ± production in the same unconventional channels. In section IV, we show the results of a multivariate analysis for the charged Higgs boson searches and show, by using the results of a genetic algorithm scan, that our proposed theory can produce the type of signals visible with such an analysis at the LHC. Finally, we summarize and conclude our results in section VI.

II. THE MODEL
In this work, we propose a 3HDM, with features that lead to a simple yet predictive theory.
The model has a U(1) X × U(1) Z global symmetry constraining its scalar potential. This symmetry is the biggest abelian symmetry not leading to additional accidental symmetries in a 3HDM [5,12,13]. As a consequence, in the limit of one VEV being much larger than the other two, we can derive simple analytical formulas for masses and mixing matrices in the scalar sector and readily understand the features of the model and its physical consequences.
The U(1) X ×U(1) Z is also present in the fermion sector of the theory. We choose the charge assignments to constrain the Yukawa sector in a manner consistent with the experimental hierarchies in the quark mass spectrum while forbidding tree-level FCNCs arising from the scalar sector. The upside of this is that the mass hierarchy is directly connected to a VEV hierarchy, which needs not to be as strong as the hierarchy in the SM Yukawa parameters to explain the known quark masses.
A nice consequence is the opening of new search strategies for testing this model at the LHC. Due to the structure of the Yukawa sector, two physical charged Higgs bosons would be produced mainly through cs fusion, which can lead to good signal-to-background ratios as will be shown in section IV.
A. VEV hierarchy and the softly broken U(1) X × U(1) Z symmetry Besides the field content of the SM, the model has two additional scalar SU(2) L doublets for a total of three. We will denote them by H 1,2,3 , with charges shown in Table I, and expand around the vacuum as We will often focus on the case where v 3 v 1,2 . This particular limit calls for the definition of a small parameter ξ, After spontaneous symmetry breaking, in the limit ξ → 0, there remains a U(1) X × U(1) YZ symmetry where U(1) YZ is generated by a combination of the U(1) Y and U(1) Z generators. That means that all processes violating U(1) X × U(1) YZ (and in particular U(1) X ) would be suppressed by some power of ξ. As we will see, in the limit that ξ 1 it is possible to derive simple expressions for the masses and mixing matrices in the scalar sector. It is worth mentioning at this stage that, while such expressions serve as tools to understand the model's features, all scalar masses and mixing matrices will be computed fully numerically (i.e. not as expansions in ξ) when scanning the parameter space of the model.
A spontaneously broken U(1) X ×U(1) Z global symmetry would lead to massless Goldstone bosons and constrain the model significantly when considering e.g. the precise measurements of the Z-boson width. This motivates us to softly break the symmetry by adding additional mass terms in the scalar potential. The scalar potential consistent with a softly broken U(1) X × U(1) Z global symmetry group can be split into fully symmetric and soft-breaking parts as with All parameters in the scalar potential can be taken real without any loss of generality.
This is due to the fact that the parameters in V 0 are real by construction, while any phases on m 2 ij can be eliminated by field redefinitions of the three Higgs doublets. As a consequence the scalar sector of the model has no choice but to be CP-conserving.
For convenience we defineλ Finally, assuming that v 1,2,3 = 0 and requiring that the first derivative of V vanishes, we can write B. Extending the U(1) X × U(1) Z to the fermion sector We assign the quark U(1) X ×U(1) Z charges such that the neutral component of H 3 couples to only up-and down-type quarks of the third generation while the neutral components of H 1 and H 2 couple to the first and second generation down-type and up-type quarks respectively, i.e.
In this way we forbid scalar-mediated tree-level FCNCs and simultaneously enforce a Cabibbo-like quark mixing, where the gauge eigenstates of the third quark family are aligned with the corresponding flavour eigenstates. This also means that a hierarchy in the VEVs of the Higgs doublets, where v 3 v 1,2 , leads to a third quark family that is much heavier than the first two without a strong hierarchy in the Yukawa couplings. In Table I, we show the most general quark charge assignments allowing the terms in Eq. (8) once the U(1) X × U(1) Z charges of H 1,2,3 are fixed. As long as the parameters α, β, γ and δ in Table I satisfy the terms in Eq. (8) are also the only allowed quark Yukawa interactions. It is worth noting that in the mass basis, the free parameters in the quark sector are simply the quark masses and the Cabibbo angle. The reader might note that at higher orders, the Yukawa interactions only allow for a mixing between the first and second quark generations, thus opening the question of how to reproduce the observed full CKM mixing in the quark sector. As this model is thought as an effective theory, one can write the following dimension-6 operators consistent with the imposed symmetries Such terms will induce naturally small (suppressed by a scale of new physics) mixing terms with the third quark family once Higgs VEVs appear. The operators can in principle be generatedà la Frogatt-Nielsen mechanism [30] by integrating out the heavy fields of a high-energy theory. A deeper analysis of this is beyond the scope of this paper.
Finally, we note that the lepton Yukawa sector can be made very SM-like by assigning the lepton U(1) X × U(1) Z charges such that they only couple to H 3 . We will assume that this is the case throughout this work, and will not discuss the implications on lepton phenomenology any further. However, we want to point out that there are also other interesting scenarios, e.g. where the leptons couple to H 1,2,3 such that the lepton mass hierarchies are also related to v 1,2 v 3 .

C. The spectrum, mixing matrices and interactions of the scalar sector
After spontaneous symmetry breaking, the mass terms in the scalar potential V in Eq. (3), can be neatly expressed using We note that both M 2 C,P have an eigenvector ∝ v i with vanishing eigenvalue. The corresponding Goldstone states become the longitudinal polarization states of the massive electroweak gauge bosons.
The electrically neutral scalar, pseudo-scalar and charged scalar mass eigenstates, are related to the interaction eigenstates as The states A G and H ± G in Eq. (13) denote the Goldstone bosons. Working in the ξ 1 limit, the mixing matrices S, P and C are identical up to O(ξ) but differ at O(ξ 2 ). It is here convenient to define an angle β ∈ [0, π 2 ] as To the second order in ξ, we have with For the O(ξ 2 ) pieces, we have Here, Z 1,2,3 parametrize the leading order difference between the mixing matrices, which will be important as these parameters determine the off-diagonal scalar-scalar interactions with the electroweak gauge bosons. We also note that as X, Y , Z 1 , Z 2 , Z 3 get larger, the expansion in ξ becomes less reliable.
The state h 125 contains mostly h 3 , meaning that it couples substantially to the third quark family. It also receives a mass of the order of v 3 ∼ v, making this state our candidate for the observed SM Higgs-like 125 GeV state. The exotic scalars h a,b , A a,b and H ± a,b can all be made heavy as the leading order contribution to their masses is inversely proportional to ξ. To this leading order, {h a , A a , H ± a } are degenerate in mass. This is also the case for More accurately, the masses are given by to O(ξ). This means that the exotic scalars and pseudo-scalars are typically very close . Note also that the couplings λ ij can either be positive or negative, such that the charged scalars H + a,b can both be heavier and lighter than the neutral scalars in the respective family.
We conclude this section by listing the trilinear interactions between the physical scalars and the electroweak gauge bosons, as they are relevant for the collider phenomenology of the charged Higgs boson discussed in section III. The interactions between the neutral scalars, charged scalars and the W boson are given by with The top line in Eq. (22) is written in the interaction eigenbasis of the scalars, while the bottom line is the same expression in terms of the mass eigenstates. The W boson also couples to pairs of charged scalars and pseudo-scalars as with Similarly, for the trilinear interactions with the Z boson, we have where c W is the cosine of the Weinberg angle and

D. Scalar-fermion couplings
Knowing the mixing matrices S, P and C for the neutral scalars, pseudo-scalars and charged scalars, respectively, to the first orders in ξ, it is straightforward to obtain the Yukawa interactions between the physical scalars and the quarks. Using Eq. (14) we find that h 125 couples to quarks in a way similar to the SM, For the third quark family, this is an obvious consequence of the model's symmetries, as t and b quarks receive their masses from H 3 with v 3 v, and h 125 is mostly made of h 3 .
On the other hand, the first and second family get their masses from H 1,2 with v 1,2 v 3 , so the corresponding Yukawa couplings with the gauge eigenstates H 1,2 are quite large as . When shifting to the mass eigenbasis, h 1,2 contribute to h 125 only at O(ξ) thus giving an overall coupling of O(m q /v 3 ).
In the same process, we also find the interaction terms between the quarks and the exotic scalar states h a,b , A a,b and H ± a,b . Couplings to the third quark family are generally quite small ∼ m t,b ξ/v 3 . In our model, phenomenologically the most relevant couplings are with the second quark family instead, which to the leading order in ξ read where θ C is the Cabbibo angle. When the masses of the scalars are in the appropriate range, we can expect that the charged scalars H + a,b would be produced in collider experiments through cs fusion while h a and A a (h b and A b ) would mainly be produced by the ss (cc) fusion.

III. A MODEL INDEPENDENT APPROACH
One of the interesting features of our model is the existence of heavy charged scalars H + (H − ) that mostly couple to a cs (cs) pair as their interactions with tb (tb) are small due to the model symmetries. Furthermore, we find that H ± can decay to a W ± h 125 pair with a sizable branching ratio (BR) which is still allowed by the current experimental data. It turns out that this unconventional channel, while not explored in the literature before for m H ± > 200 GeV, can be a rather clean way to search for charged scalars at the LHC.
In the following, we adopt a model independent approach in searching for charged scalars exhibiting those features. In section V, we will show how the analysis can be used to find discovery regions in the parameter space of the 3HDM we have proposed above. We take a model independent approach to not only test the predictions of our model, but also to offer a guideline for our experimental colleagues to implement this new search channel in the experimental analyses.
We start with the following model independent Lagrangian for H ± including its kinetic (L kin ) and interaction (L int ) terms There are four free parameters in the above Lagrangian viz. the charged Higgs mass m H ± , and the three couplings κ p cs , κ m cs and κ W h 125 . In general, κ p cs and κ m cs both could be non-zero. In that case, the production cross section, σ(pp → H ± ) is proportional to the combination (κ p cs ) 2 + (κ m cs ) 2 . Therefore, instead of two free couplings, we introduce a single free parameter κ cs which is, κ 2 cs = (κ p cs ) 2 + (κ m cs ) 2 . From the above model independent Lagrangian, we see that H + has only two decay modes: W + h 125 and cs. The corresponding tree-level partial widths are given by where m h 125 = 125 GeV. The expression Γ(H + → cs) is given in the limit of massless c and s quarks. In general, H ± can have other decay modes too. We, therefore, take the BR of the decay mode H ± → W ± h 125 denoted by BR W h 125 as a free parameter instead of κ W h 125 . So, one can write the following in the narrow-width approximation, where σ 0 (m H ± ) is the cross section of pp → H ± for κ cs = 1. We show σ 0 (m H ± ) at the LHC ( √ s = 13 TeV) as a function of m H ± in Fig. 1.

IV. SEARCH FOR CHARGED SCALARS PRODUCED BY cs FUSION
We implement the model independent Lagrangian of H ± as shown in Eqs. (30) and (31) in FeynRules [31]

A. Signal
We focus the H + (H − ) production from the cs (cs) initial state followed by the decay We consider a semileptonic final state where W ± decays leptonically and h 125 decays to bb. Therefore, the chain of the signal process in our case is Here, = {e, µ}. We then have one charged lepton, two b-jets and missing transverse energy in the final state and our event selection criteria is exactly one charged lepton (either an electron or a muon including their anti-particles), at least two jets and missing transverse energy that pass the following basic selection cuts:  • 1b-tag: In this category, we demand at least one b-tagged jet among the two leading p T jets.
• 2b-tag: In this category, we demand that both the two leading p T jets are b-tagged.
This category is a subset of the 1b-tag category.
To reconstruct the Higgs boson, we apply an invariant mass cut |m H ± − m h 125 | < 20 GeV around the Higgs boson mass m h 125 = 125 GeV. However, the full event is not totally reconstructible due to the presence of the missing transverse energy.

B. Background
The main background for the signal with one lepton, at least one or two b-tagged jets and missing energy can come from the following SM processes: 1. W ± + jets: The definition of our inclusive W ± + jets background includes up to two jets and we include the b parton in the jet definition i.e. j = {g, u, d, c, s, b}.
We generate these background events in two separate parts. In one sample, we only consider light jets i.e.j = {g, u, d, c, s} and combine pp → W ± + (0, 1, 2) j processes where we set the matching scale Q cut = 25 GeV. This background is the largest (the cross section is about 1.53 × 10 5 pb at the LHC, with √ s = 13 TeV, without any cut) among all the dominant SM backgrounds we have considered. Although the bare cross section is large, it will reduce drastically after b-tagging due to a small mistagging (light jet is tagged as b-jet) rate. We find that its contribution in the 1b-tag category is substantial but in the 2b-tag category is very small. In the other sample, we consider at least one b parton in the final state where we combine pp → W ± bj and pp → W ± bb processes (no SM pp → W ± b process exists). This background will contribute significantly in both the categories. We include pp → W ± h 125 → W ± bb and pp → W ± Z → W ± bb processes in the pp → W ± bb channel.
2. tt + jets: The definition of our inclusive tt + jets background includes up to two jets containing also b partons. We generate this background by combining pp → tt+(0, 1, 2) j processes using the matching scale Q cut = 25 GeV. The matched cross section is about 431 pb before the top decay and without any selection cut applied.
We find that this background is the dominant one after the strong basic selection cuts (applied before passing the events to the MVA).

Single top: This background includes three types of single top processes -s-channel
single top (such pp → tb), t-channel single top (i.e. pp → tj) and single-top associated with W (such as pp → tW ± ) processes. Note that for the pp → tW ± process, the selected lepton can come from two possible ways, either from the decay of the associated W ± or from the W ± coming from the top decay. These two possibilities are properly included in our event sample. The single top background also contributes significantly to the total background. 4. Diboson: This background includes pp → W ± W ∓ → W ± + jj and pp → W ± Z → W ± + jj processes where two light jets come from the decay of W or Z bosons. In this background, we have also included pp → W ± Z → W ± νν processes where two selected jets come from the parton showers. This background reduces drastically due to the small mistagging efficiency of light jets that are misidentified as b-jets.
Finally, in the MVA this background contributes negligibly to the total background.
Note that two diboson production processes viz. pp → W ± h 125 → W ± bb and pp → W ± Z → W ± bb are already considered in the W + jets background.
5. QCD multijets: The multijet background arises due to QCD interactions at the LHC and has a very large production cross section, especially in the soft region.
The QCD-induced multijet production processes can potentially contribute to the total background for our signal by faking the lepton, / E T and b-tagged jets. It is impractical to study this part of the background using a Monte-Carlo simulation since it is computationally challenging to generate enough events due to very low fake rates. In experimental analyses, this contribution is usually estimated from the data. In our analysis, we do not consider this background since it will be largely diminished after strong preselection cuts and will be further reduced due to small fake rates of the considered final states.
The SM background, especially the W + jets component, is large and therefore one has to design a clever set of cuts which would notably reduce such a background but would not notably affect the signal. This implies that the cut efficiency for the background is very small and hence, a large number of background events has to be generated. In order to avoid the generation of a large event sample, we apply a strong cut on the partonic center-of-mass energy, √ŝ > 200 GeV at the generation level of all background processes.
This cut can reduce the W + jets background by two orders of magnitude. However, this cut has no or very little effect on the other backgrounds viz. tt + jets, single top and diboson ones since the threshold energy for them is either above or slightly below 200 GeV. In the case of a signal, √ŝ is always above 200 GeV since we are interested in the parameter space regions where m H ± > m W + m h 125 205 GeV.
One should note that, in reality, the full reconstruction of √ŝ of an event is not possible if there is missing energy present in that event. In this case, one can construct an inclusive global variable √ŝ min defined in Ref. [40] which is closest to the actual √ŝ of the event.
One can roughly approximate √ŝ ≈ √ŝ min if there is only one missing neutrino in the event but this approximation gets poorer with the increase of the number of neutrinos in the final state. For simplicity, we have used the cut √ŝ > 200 GeV at the generation level. But in reality, one can use a cut on √ŝ min to trim the background before passing it for further analysis.

C. Multivariate analysis
A W h 125 resonance, similar to our case, can also appear from the decay of a heavy charged gauge boson, W . The search for W in the ± + / E T + bb channel (same final state that we are interested in) has been carried out at the LHC [41,42]. In these searches, they mainly focus in the TeV-scale W mass and the analyses are done using the cut-based techniques. A cut-based analysis may not perform well in our case, especially for low m H ± region due to the presence of a large SM background [43,44]. Therefore, we choose to use a MVA to obtain a better signal-to-background discrimination which usually leads to a better significance than a cut-based analysis. See Ref. [45] for a brief review on various multivariate methods and their use in collider searches. In this paper, we only use multivariate techniques and do not compare our achieved sensitivity with the cutbased techniques.
We choose the following twelve simple kinematic variables that are also listed in Table III for our MVA.
• Transverse momenta of lepton, p T ( ) and two leading-p T jets, p T (J 1 ) and p T (J 2 ).
• Missing transverse energy / E T and pseudorapidity of / E T vector denoted by η( / E T ).
• Scalar sum of transverse momenta of all visible particles denoted by H T .
• Invariant mass of two leading-p T jets denoted by M (J 1 , J 2 ).
These variables are chosen by comparing their distributions for the signal generated for m H ± = 300 GeV with the total background distributions. They are selected from a bigger set of variables based on their discriminating power and less correlation. In Fig. 2, we show the normalized distributions of these variables for the signal with m H ± = 300 GeV and the total background. Similar distributions for m H ± = 500 GeV are shown in Fig. 3. From these figures, one can see that each of these distributions has reasonable discriminating power between the signal and the background. We use these kinematic variables simultaneously in a MVA whose output shows large differences in their shapes for the signal and the background. One should notice that the signal distributions deviate more from the background ones as we increase m H ± . Therefore, isolation of the signal from the background becomes easier for heavier resonances. We, therefore, tune our MVA for lower masses and use the same optimized analysis for larger masses. These numbers are obtained for m H ± = 300 GeV for the 2b-tag category. These numbers can vary for other choices of parameters.
In Table III The corresponding cut efficiencies are shown as functions of C in Fig. 4b (Fig. 4d) for m H ± = 300 GeV (m H ± = 500 GeV). The optimal BDT cut (BDT opt ) is defined for which fb −1 ) before the BDT analysis, it is possible to achieve a maximum 5σ significance for BDT opt 0.26. After this cut, the number of signal events is reduced to 118 from 222 but the background events are drastically reduced to 436 from 33031. In Table IV, we show N S and N B along with N bc S , the minimum number of signal events before the BDT cut that is required to achieve 5σ significance, for different m H ± values and for the two selection categories. events that can be discovered with 5σ significance using our MVA is denoted by N bc S (this is before the optimal BDT cut as shown in the third column). The signal and background events that survived after the optimal BDT cut are denoted by N S and N B , respectively, and they lead to 5σ significance.
the parameter space where that is the case, which shows that if limits are set by the experimental collaborations, the theory can be further constrained using the current experimental data.
The first task is to match our model to the Lagrangian in Eqs. (30) and (31). For each parameter space point, we choose the lightest charged scalar for the analysis. Although we concentrated our search in the parameter space region with v 1,2 v 3 , as to exploit the SM-like h 125 state in that limit, we do not rely on the validity of the expansion in small ξ in this analysis. The couplings κ cs and κ W h 125 are found in Eq. (22) after a numerical calculation of the spectrum and mixing matrices. To find the discovery reach of our parameter space, we translate N bc S in terms of the model parameters by using the following relation where σ is the cross section after showering and hadronization, S is the signal cut efficiency and L is the integrated luminosity.
For the calculation of BR W h 125 , it is important to note that although in general H ± can decay to W ± h a,b,125 , we are interested only in the decay mode involving h 125 , as in our model this is the state that couples the strongest to bb (see Eqs. (28) and (29)).
In addition, our model must be able to pass several consistency tests in order to be phenomenologically viable, such as reproducing the electroweak precision measurements.
The original formulation [46] for BSM contributions to the electroweak precision observables in terms of the S, T and U parameters assumes that the scale of new physics is 1 TeV. As our model allows for new exotic scalars to have masses around the electroweak scale, we must employ the more general formalism introduced in Refs. [47,48] with an extended set of oblique parameters S, T , U , V , W and X. These can then be used to calculate S , T and U for which the standard Z-pole constraints on S, T and U apply.
To compute S , T and U , we have applied the results in Ref. [49], in which S, T , U , V , W and X are computed for a general N -Higgs Doublet Model with the inclusion of arbitrary numbers of electrically charged and neutral SU(2) L singlets. To summarize, when scanning the model parameter space for phenomenologically interesting regions, we look for points for which the following constraints are satisfied: • There are no tachyonic scalar masses and the scalar potential is bounded from below (the corresponding constraints on the quartic couplings can be found in Ref. [12] taking into account that our λ ii differ by a factor two from theirs).
• The SM Higgs-like scalar has a mass no more than 5 GeV away from the observed 125 GeV value, and has a Yukawa coupling to the top quark satisfying |y tth 125 | ∈ [0.9, 1.1].
• The exotic decays Z → h a,b A a,b are kinematically forbidden, as to not be in conflict with the precision measurements of the Z width.
• The lightest charged Higgs has a mass in the range [m • The computed values of S , T and U fall within the error bars on S, T and U as reported in Ref. [50].
• The value of κ 2 cs × BR W h 125 is at least 0.5 above the 100 fb −1 discovery threshold for the 1b-tag category set by the MVA.

A. Scanning the parameter space
A random scan over the parameter space of the theory is both computationally expensive and not efficient. A good alternative, without the need for sophisticated statistical methods but still very powerful is the use of Genetic Algorithms (GA).
Following the guidelines set in Ref.
[51], we wrote a GA in Mathematica for finding the parameter points in the discovery region, with a fitness function taking into account all the constraints listed above and including the so-called biodiversity enhancement to explore the parameter space more thoroughly.
GAs start from a randomly generated initial population, with each full cycle resulting in a new generation of candidates. The fittest parameter points are selected for every generation and their parameters are modified (by crossover and/or mutations) leading to a new generation. The new candidate points are then used in the next iteration of the GA.
The GA finishes when either a maximum number of generations or a satisfactory fitness level is reached. We decided to build the GA relying on mutations only as it usually performs comparably to GAs including a crossover but it is simpler to implement, and it was stopped once a given number of valid parameter points was reached.

B. Results of the GA parameter scan
We performed five independent scans with different initial population sizes ranging from 50 to 1000, with varying mutation rates and different lower limits on m H ± . We found 2116 parameter space points of the proposed model satisfying all the constraints within the discovery region of our analysis. In Figs. 5a and 5b, we show the 5σ discovery contours of κ 2 cs × BR W h 125 corresponding to 1b-and 2b-tag categories, respectively, as functions of m H ± for L = 50, 100 fb −1 at the LHC ( √ s = 13 TeV). Here, these functions are overlayed with the corresponding values for the parameter points found by the GA scanning procedure. We find that both selection categories are almost equally sensitive in probing the parameter space of our model. However, the 2b-tag category is slightly more sensitive than the 1b-tag category since the background reduction is better for the former. The irregularities in the charged Higgs mass dependence seen in Figs. 5a and 5b are due to a combination of points from scans with different lower limits on m H ± .
As discussed before, since a W can also produce a W h 125 resonance, we compare our reach with the W h 125 resonance search data. In Fig. 5b, the shaded region is excluded by the ATLAS W h 125 resonance search data [41] in the + / E T + bb channel. To obtain this, we translate the 95% confidence level (CL) upper limit (UL) on the cross section set by ATLAS in terms of our model parameters by using the following relation, where W and H ± are the cut-efficiencies for the W and H ± respectively and they are different, in general. For simplicity, we assume W = H ± while obtaining the exclusion region on our model parameters. For instance, for m H ± = 800 GeV, κ 2 × BR W h 125 2 × 10 −3 is excluded with 2σ CL using L ≈ 36 fb −1 data but κ 2 × BR W h 125 2 × 10 −3 region can be discovered with 5σ significance if we go to a higher luminosity. The exclusion region starts from m H ± = 500 GeV since the latest data used here are available for W masses above 500 GeV.  It is worth noting that although the GA did not rely on the validity of the ξ 1 expansion, it often found points where that is the case. Although the initial populations had a hierarchy in the VEVs of the scalar fields (v 1,2 v 3 ), the GA had no inherent constraints stopping it from exploring the regions without it. Notably, the majority of the found points did show that feature and therefore a diminished hierarchy in the Yukawa couplings of the quark sector together with all the features described in section II. That can also be seen in Figs. 5a-5b where we have indicated the specific values of ξ for the valid parameter points. We have checked the difference between the O(ξ) expressions for masses and the full numerical calculation, and find that for a vast majority of the valid parameter points the ξ 1 expansion is reliable.
In Fig. 7, we show distributions of quartic couplings, values of ξ and mass parameters in the scalar potential, as well as the Higgs VEVs, for points in the discovery region found by the GA. Indeed, the typical values of v 1,2 are likely to be below 100 GeV, with v 2 extending over a larger domain than v 1 , while v 3 values are mostly concentrated close to the maximal 246 GeV limit. There is still a small number of valid points incompatible with the ξ-expansion due to the smallness of one of m 2 ij . Such points would still be in the discovery region, although without the features that assume ξ 1. As can be seen from Fig. 6b, the masses of the exotic scalars and pseudo-scalars tend to align as predicted by the ξ 1 expansion (see Eq. (21)).

VI. SUMMARY AND CONCLUSIONS
We have, in this article, introduced a class of 3HDMs with a global U(1) X × U(1) Z family symmetry that is softly broken by bi-linear terms in the scalar potential. We have shown how to assign the X and Z charges of the quarks such that no tree-level FCNCs are present, while enforcing the Cabibbo-like structure of V CKM . We described how a mixing with the third quark family can be induced from dim-6 operators, which would explain the smallness of the corresponding entries in V CKM . Moreover, we showed that a hierarchy in the VEVs of the three Higgs doublets, v 1,2 v 3 , leads to a heavy third quark family without the need for a strong hierarchy in the Yukawa couplings (contrary to what happens in the SM where e.g. y up /y top ∼ 10 −5 ). The same hierarchy has been exploited to derive simple closed expressions for the scalar masses and mixing matrices by expansions in the small parameter ξ ≡ v 2 1 + v 2 2 /v 3 1.
A generic prediction of the model is that the new scalars h a,b , A a,b and H ± a,b are likely to couple strongly to the s and c quarks, yielding different signatures in colliders at variance with the standard searches focusing on the third quark family. As an example, we studied collider phenomenology of the lightest charged Higgs when its mass is in the 250 -1000 GeV range, under the assumption that the other charged Higgs is sufficiently heavy to be dropped out of the analysis. In that case, the lighter charged Higgs would be resonantly produced through a cs fusion, and, for certain regions of the parameter space, subsequently decay to W h 125 . All other decay channels are assumed to only contribute to its total width.
We particularly focused on one of the possible channels -the cs → H + → W + h 125 channel, which has not been explored before in the context of heavier charged Higgs searches. This channel is specific to our class of 3HDMs and is particularly sensitive to the sub-TeV charged Higgs mass and small-ξ regions. We showed that this unconventional channel, when combined with the power of a multivariate analysis, leads to good signal-tobackground ratios even for masses below 500 GeV and thus can be used to probe models with that particular feature at the LHC. We employed a model independent formulation so that our approach can be applied to any model which predicts a sufficiently large cross section for the cs → H + → W + h 125 process to be observed in the future LHC runs.
Our analysis can also be applied to improve sensitivity for W searches especially for the sub-TeV masses.
We then used a genetic algorithm to find parameter space points in our 3HDM which would yield signals with > 5σ significance, while still satisfying the standard phenomenological constraints. Although the scan did not rely on ξ 1, a vast majority of the points were consistent with that limit and thus showed all the features mentioned above and described in section II. This shows that the described unconventional search strategy can effectively probe realistic multi-Higgs theories with the current LHC data, and so we think it should be seriously considered by our experimental colleagues.