Dark matter in three-Higgs-doublet models with $S_3$ symmetry

Models with two or more scalar doublets with discrete or global symmetries can have vacua with vanishing vacuum expectation values in the bases where symmetries are imposed. If a suitable symmetry stabilises such vacua, these models may lead to interesting dark matter candidates, provided that the symmetry prevents couplings among the dark matter candidates and the fermions. We analyse three-Higgs-doublet models with an underlying $S_3$ symmetry. These models have many distinct vacua with one or two vanishing vacuum expectation values which can be stabilised by a remnant of the $S_3$ symmetry which survived spontaneous symmetry breaking. We discuss all possible vacua in the context of $S_3$-symmetric three-Higgs-doublet models, allowing also for softly broken $S_3$, and explore one of the vacuum configurations in detail. In the case we explore, only one of the three Higgs doublets is inert. The other two are active, and therefore the active sector, in many aspects, behaves like a two-Higgs-doublet model. The way the fermions couple to the scalar sector is constrained by the $S_3$ symmetry and is such that the flavour structure of the model is solely governed by the $V_\text{CKM}$ matrix which, in our framework, is not constrained by the $S_3$ symmetry. This is a key requirement for models with minimal flavour violation. In our model there is no CP violation in the scalar sector. We study this model in detail giving the masses and couplings and identifying the range of parameters that are compatible with theoretical and experimental constraints, both from accelerator physics and from astrophysics.

1 Introduction consistent with observations for the DM candidate with masses 55-74 GeV [65]. Although the high-mass region, above 500 GeV, requires severe fine-tuning it is still compatible with data. The region with the DM masses 74-500 GeV is still possible: the exclusion of a given mass region in the parameter space depends on whether one requires the IDM to provide a single component DM candidate to generate all relic density or just partially contribute to the DM relic density. More recently, Ref. [66] confirmed earlier mass-range observations of the IDM and provided a set of benchmarks for e + e − studies, for IDM detection at potential e + e − colliders (ILC, CLIC) see Refs. [67][68][69]. The aforementioned DM mass ranges are sketched in figure 1, together with those of some 3HDMs and that of the model (R-II-1a) explored here.  Figure 1: Sketch of allowed DM mass ranges up to 1 TeV in various models. Blue: IDM according to Refs. [65,66], the pale region indicates a non-saturated relic density. Red: IDM2 [70]. Ochre: 3HDM without [71][72][73] and with CP violation [74]. Green: the R-II-1a χ model presented in this paper.
The main reasons some mass regions are excluded, are qualitatively as follows: • Below m DM = O(50) GeV, it would annihilate too fast in the early Universe via offshell Higgs bosons decaying to fermions. For m DM > ∼ 5 GeV the leading contribution is from decay into b-quarks. Also, much of this range is excluded by LUX [75], by PandaX-II [76], and by XENON1T [77] data.
• In the case of small mass splitting, resonant coannihilation into gauge bosons occurs for 40-45 GeV. Another resonant coannihilation, mass splitting independent, at m DM ∼ 60 GeV is due to annihilation into the Higgs boson, dependent on λ L .
• In the region between O(80) GeV and O(m h ) it annihilates too fast via a pair of gauge or Higgs bosons. The latter channel, however, depends on λ L .
• In the region from O(m h ) to O(500) GeV it annihilates too fast via two on-or off-shell gauge bosons. Depending on the λ L value, annihilation into a pair of Higgs bosons can become significant. Also, in this region, annihilation into t-quarks becomes available.
• Above O(500) GeV an important mechanism is due to the coupling to longitudinal vector bosons [36]. With an increasing DM mass the relic density would grow. At some tens of TeV, interactions would become non-perturbative as a large value of λ L would be required to match the DM relic density value. Coannihilations between all scalars from the inert doublet become important in the limit of degeneracy.
The exact mass values where these phenomena set in, depend on details of the model and the adopted experimental constraints.
2. By having one non-inert doublet along with two inert doublets.
This latter approach has been pursued by various groups [71][72][73][74][81][82][83][84][85]. The models studied in Refs. [81,83] assume an S 3 -symmetric 3HDM potential, and take the S 3 singlet to represent the SM-like active doublet. There are then two inert doublets with zero vacuum expectation value (vev), corresponding to R-I-1 in our terminology in table 1 below. This model has three mass degeneracies: in the charged sector as well as in the CP-odd and CP-even neutral sectors of the S 3 doublet, as discussed in detail in Ref. [86]. In all these models, like in the original IDM, the lightest neutral scalar in the inert sector is prevented from decaying to SM particles by a Z 2 symmetry. In Refs. [81,83] the Z 2 symmetry is softly broken, in order to lift the degeneracy, leading to consistent models with mass, at least, in the range 40-150 GeV [83].
In Refs. [71][72][73], a Z 2 -symmetric potential is constructed, and again the vacuum with two vanishing vevs is studied, yielding a model with two inert doublets. The authors point out that having more inert fields allows for more coannihilation channels, and opens up the parameter space compared to the IDM. Consistent models are found with mass in the range 53 to 77 GeV. Also, it was found that the high-mass region, which for the IDM requires a mass above some 500-525 GeV, can give consistent models with mass down to 360 GeV [72].
In Ref. [74], the Z 2 -symmetric 3HDM is studied, with two inert doublets and a new ingredient being explicit CP violation via complex coefficients in the Z 2 -symmetric potential. This allows for more parameters, but only a modest extension of the mass ranges allowed by the CP-conserving 3HDM discussed above is found.
Finally, one can stabilise DM by a Z 3 symmetry, rather than the Z 2 symmetry [84,85]. In this case, two-component DM is considered.

The S -symmetric models
Our philosophy here is to first identify vacua with at least one vanishing vev as possible frameworks for a DM model. We shall work in the irreducible representation, where we have two SU(2) doublets, h 1 and h 2 in a doublet of S 3 , and one S 3 singlet denoted h S . The Higgs doublets with a vanishing vev are labelled as "inert", since their lightest member (assumed neutral) could be DM. They are listed in table 1 below. Some of these are stabilised by a surviving symmetry of the potential, whereas others would require an imposed Z 2 symmetry. The general S 3 -symmetric scalar potential is invariant under Z 2 : h 1 → −h 1 . We do not consider other mechanisms which would stabilise the DM candidate.
Some of these vacua are associated with massless states, hence we shall allow for soft breaking of the S 3 symmetry in the potential [86], noting that soft breaking is not possible in the Yukawa sector. When introducing soft breaking terms, the constraints will change. However, we will retain the nomenclature of the unbroken case from which they originate, thus when adding soft-breaking terms to R-I-1, we denote it r-I-1.

The scalar potential
In terms of the S 3 singlet (1 : h S ) and doublet (2 : (h 1 h 2 ) T ) fields, the S 3 -symmetric potential can be written as [87][88][89]: Note that there are two coefficients in the potential that could be complex, thus CP can be broken explicitly. For simplicity, we have chosen all coefficients to be real. In spite of this choice there remains the possibility of breaking CP spontaneously.
In the irreducible representation, the S 3 doublet and singlet fields will be decomposed as where the w i and w S parameters can be complex. The symmetry of the potential can be softly broken by the following terms [86]: In accordance with the previous simplification of couplings we assume that the soft terms are real. Soft symmetry-breaking terms are required whenever we work with solutions with λ 4 = 0 since in this case most vacua lead to massless scalar states, Goldstone bosons of an O(2) symmetry resulting from this choice. In our analysis we shall not make use of these terms. We recall that in order for a doublet to accommodate a DM candidate it must have a vanishing vev, since otherwise it would decay via its gauge couplings, e.g., the SW + W − and SZZ couplings.

The Yukawa interaction
Whenever the singlet vev, w S , is different from zero we can construct a trivial Yukawa sector, L Y ∼ 1 f ⊗ 1 h . In this case, the fermion mass matrices are: where y's are the Yukawa couplings of appropriate fermions. Another possibility is when fermions transform non-trivially under S 3 , with a Yukawa Lagrangian written schematically as Such structure yields the mass matrix for each quark sector (d and u) of the form If, for simplicity we assume all y's to be real, then a complex Cabibbo-Kobayashi-Maskawa (CKM) matrix would have to be generated from complex vacua. As we point out in the sequel, in some cases realistic quark masses and mixing can only be generated if the quarks are taken to be S 3 singlets and only couple to h S , which requires complex Yukawa couplings as in the SM. On the other hand, even in some of the cases when the vevs are complex, the CKM matrix is not always complex. The hermitian quantity, is going to be useful in our discussion. The fermion mass matrices M f are diagonalised in terms of the left-handed and the right-handed fermion rotation matrices, which, in general, are not equal. However, the quantity M f M † f is diagonalised in terms of the left-handed rotation matrix only, and M † f M f accordingly in terms of the right-handed rotation matrix. The eigenvalues of H f will be the squared fermion masses.
It is instructive to count the number of parameters available in these models, and compare with the number of physical quantities to be fitted. The full Yukawa sector, eq. (3.5), has 10 parameters (5 y d i and 5 y u i ). If w S = 0, this number is reduced to 6 (y = 0). If either w 1 = 0 or w 2 = 0, we still have 10 parameters, but if both are zero, w 1 = w 2 = 0, the number is reduced to 4. The available parameters will be "used" to fit 6 quark masses plus 4 parameters of the CKM matrix. Apart from that, the most general vacuum configuration is given by 3 absolute values and 2 complex phases.
As we are interested in a DM candidate, one of the requirements is to have a field with a vanishing vacuum expectation value. Let us briefly consider what happens with the Yukawa sector in this case. When the DM candidate resides in the scalar S 3 singlet, w S = 0, we need the fermions to couple to the S 3 doublet, schematically represented by L Y ∼ 2 f ⊗ 2 h . In some cases, the DM candidate can reside in the scalar S 3 doublet. To keep the notation simple, we still write the Yukawa sector as L Y ∼ (2 ⊕ 1) f ⊗ (2 ⊕ 1) h , as the general form of the fermion mass matrices persists. However, in order to stabilise the DM candidate one needs to introduce an additional Z 2 symmetry in the Yukawa sector to decouple a specific inert doublet from the fermionic sector.

Vacua with at least one vanishing vev
In table 1, we list the vacua of the S 3 -symmetric potential that could accommodate DM [6]. The vacua will be represented in the form Whenever h 0 S = 0 we indicate that the fermions may transform trivially under S 3 . This is the simplest choice. In all other cases we need them to transform non-trivially in order to acquire masses. Whenever λ 4 = 0, the scalar potential acquires an additional O(2) symmetry between h 1 and h 2 , which could be spontaneously broken. This breaking leads to one massless neutral scalar, see Ref. [86] for a classification of massless states. Another interesting feature is that the scalar potential with the λ 4 = 0 constraint will have two additional Z 2 symmetries beyond h 1 → −h 1 , involving h 2 and h S . The symmetry h 1 → −h 1 is a symmetry of the potential in the irreducible representation. If h 1 does not acquire a vev it remains as an unbroken symmetry. Column 2 lists the vevs in the irreducible representation. S 2 and S 3 symmetry in the fourth column refer to remnant symmetries explicit in the defining representation, as does "cyclic Z 3 ". For λ 4 = 0 real vacua can at most break S 3 → S 2 [90]. Whenever λ 4 = 0 the S 3 symmetry can be fully broken by real vacua.
Below, we indicate some pros and cons of different vacua of S 3 -3HDM, following the nomenclature of Ref. [6].
• R-I-1: (0, 0, w S ) This case might result in a viable DM candidate. A related case was studied in Refs. [81,83]. In order to stabilise h 2 they forced λ 4 = 0 and found that this model may result in a viable DM candidate.
Scalar sector: The DM candidate resides in h 1 and thus is automatically stabilised.
There are three pairs of mass-degenerate states: a charged pair and two neutral pairs. In principle, one could lift this degeneracy by softly breaking the S 3 symmetry of the scalar potential. The only possible soft breaking term is µ 2 2 . If this term is present, there are no mass degeneracies. In addition, it is possible to stabilise the h 2 doublet by forcing λ 4 = 0. There is no spontaneous symmetry breaking associated with the inert doublets and therefore no Goldstone states would arise due to spontaneously broken O(2).
Yukawa sector: L Y ∼ 1 f ⊗ 1 h can give realistic fermion masses. Due to freedom of parameters, this can give a realistic CKM matrix and no flavour-changing neutral currents (FCNC).
Scalar sector: The Z 2 symmetry is preserved for (h 2 , h S ) → −(h 2 , h S ), or equivalently this translates into h 1 → −h 1 , and thus we have two stabilised inert doublets.
Yukawa sector: This indicates that one of the fermion mass eigenvalues vanishes, i.e., there will be a massless fermion.
• R-I-2b,2c: (w, ± √ 3w, 0) The S 3 symmetry of the scalar potential needs to be softly broken. The Yukawa sector is unrealistic. Scalar sector: There is mixing present in the mass-squared matrix as λ 4 = 0 and thus DM is not stabilised. If we artificially put λ 4 = 0, one of the neutral states would become massless. The DM is left stabilised only if the soft breaking term µ 2 2 together with ν 2 12 are introduced. Yukawa sector: This case results in a viable DM candidate provided that the Yukawa sector is trivial. Scalar sector: The DM candidate resides in h 1 and thus is automatically stabilised.
Yukawa sector: L Y ∼ (2 ⊕ 1) f ⊗ (2 ⊕ 1) h can give realistic fermion masses. The CKM matrix splits into a block-diagonal form and thus is unrealistic. However, there is another possibility to construct the Yukawa Lagrangian of the form L Y ∼ 1 f ⊗ 1 h .
• R-II-2: (0, w, 0) The S 3 symmetry of the scalar potential needs to be softly broken. The Yukawa sector is unrealistic.
Scalar sector: Due to λ 4 = 0, Z 2 is preserved for both h 1 and h S . The λ 4 = 0 constraint results in an additional Goldstone state. The only soft breaking term which does not survive minimisation is ν 2 12 . Also, ν 2 02 cannot be the only soft breaking term as then the massless state would survive. Then, the µ 2 2 and ν 2 01 couplings are free parameter as those do not depend on the minimisation conditions. The coupling ν 2 02 would require λ 4 = 0 and therefore DM is only stabilised in h 1 . Yukawa sector: L Y ∼ 2 f ⊗ 2 h can give realistic masses. However, the CKM matrix is split into a block-diagonal form.
• R-II-3: (w 1 , w 2 , 0) The S 3 symmetry of the scalar potential needs to be softly broken. The Yukawa sector is unrealistic.
Scalar sector: Due to λ 4 = 0, Z 2 is preserved for h S . An additional Goldstone boson is present. The possible softly broken couplings are µ 2 2 and ν 2 12 . If µ 2 2 is the only term present, for consistency ν 2 12 = 0, this would require to impose either w 1 = 0 or w 2 = 0. When ν 2 12 is considered, it results in w 1 = w 2 . It is also possible to have both terms present simultaneously.
Yukawa sector: L Y ∼ 2 f ⊗ 2 h can give realistic fermion mass eigenvalues. However, the CKM matrix is unrealistic and there are no free parameters to control FCNC. The Yukawa sector results in six degrees of freedom: three from the u-type couplings and three from the d-type couplings. In case of the model with both couplings, r-II-3-µ 2 2 -ν 2 12 , there is an additional degree of freedom in terms of the ratio between the vevs w 1 and w 2 .
• R-III-s: (w 1 , 0, w S ) This case might result in a viable DM candidate provided that the S 3 symmetry of the scalar potential is softly broken.
Scalar sector: Due to λ 4 = 0, Z 2 is preserved for h 2 . An additional Goldstone boson is present. Possible soft symmetry breaking terms are µ 2 2 and ν 2 01 . Both of these couplings are unconstrained by the potential. In Ref. [6] this vacuum with w 1 complex was denoted as C-IV-a and it was pointed out that the constraints made it real, therefore there is no need for λ 7 to be zero.
Yukawa sector: L Y ∼ (2 ⊕ 1) f ⊗ (2 ⊕ 1) h can give realistic fermionic masses and the CKM matrix. However, there are FCNC present in this case. In total, there are ten Yukawa couplings and a ratio between vacuum values. Due to such high number of free parameters it might be possible to control the overall effect of FCNC. There is also a possibility to construct the Yukawa Lagrangian of the form L Y ∼ 1 f ⊗ 1 h .
Scalar sector: In order to stabilise DM in h S , we are forced to impose λ 4 = 0. In general, this model results in two neutral mass-degenerate pairs. No soft symmetry breaking terms survive and thus there are no c-I-a models.
Yukawa sector: L Y ∼ 2 f ⊗ 2 h can give realistic fermion masses. However, the CKM matrix is split into a block-diagonal form.
• C-III-a: (0,ŵ 2 e iσ 2 ,ŵ S ) This case might result in a viable DM candidate provided that the Yukawa sector is trivial. Scalar sector: The DM candidate resides in h 1 and thus is automatically stabilised.
Yukawa sector: L Y ∼ (2 ⊕ 1) f ⊗ (2 ⊕ 1) h can give realistic fermion masses. However, the CKM matrix is split into a block-diagonal form. Another possibility is to construct the Yukawa Lagrangian of the form L Y ∼ 1 f ⊗ 1 h .
• C-III-b: (±iŵ 1 , 0,ŵ S ) This case might result in a viable DM candidate provided that the S 3 symmetry of the potential is softly broken.
Scalar sector: Due to λ 4 = 0, DM is stabilised in h 2 . An additional Goldstone boson is present. The only possible soft symmetry breaking term is µ 2 2 . Yukawa sector: L Y ∼ (2 ⊕ 1) f ⊗ (2 ⊕ 1) h results in realistic fermion masses and CKM matrix. This model has eleven free parameters. FCNC are present. Another possibility is to construct the Yukawa Lagrangian of the form L Y ∼ 1 f ⊗ 1 h .
• C-III-c: (ŵ 1 e iσ 1 ,ŵ 2 e iσ 2 , 0) The S 3 symmetry of the scalar potential needs to be softly broken. The Yukawa sector is most likely unrealistic.
Yukawa sector: L Y ∼ 2 f ⊗ 2 h can give realistic fermion mass eigenvalues. When the µ 2 2 or ν 2 12 terms are considered, there are seven free parameters. Preliminary check of the model with just ν 2 12 resulted in unrealistic CKM and a non-negligible FCNC contribution [91]. Most likely, exactly the same situation arises when µ 2 2 is added. However, the c-III-c-µ 2 2 -ν 2 12 model has one additional free parameter due to unfixed vevs. Nevertheless, even if the CKM values can be fitted, there are still FCNCs that have to be controlled.
This case might result in a viable DM candidate provided that the S 3 symmetry of the potential is softly broken.
Scalar sector: Due to λ 4 = 0, DM is stabilised in h 2 . There are two massless states present in this model. Possible soft symmetry breaking terms are µ 2 2 and ν 2 01 . In case of the c-IV-a-µ 2 2 model, the overall phase gets fixed (iŵ 1 , 0,ŵ S ). Yukawa sector: L Y ∼ (2 ⊕ 1) f ⊗ (2 ⊕ 1) h can give realistic fermion masses and CKM matrix. In total, there are twelve (eleven when only µ 2 2 is present) free parameters. The FCNCs are present. Another possibility is to construct the Yukawa Lagrangian of the form All in all, there are several models with a potential DM candidate. In our classification, which is summarised below, whenever unwanted Goldstone bosons are present, we refer to the need to include soft symmetry-breaking terms. Most of the models are ruled out due to unrealistic Yukawa sector. Full analysis of the Yukawa sector is out of scope of this paper. It should be noted that the non-trivial Yukawa sector might result in nonnegligible FCNC and thus some of the models would be ruled out. Whenever the quarks transform trivially under S 3 they can only couple to one Higgs doublet, the S 3 singlet, and thus FCNC are not present. We have the following possible models, indicating where the DM candidate could reside and the Yukawa sector: Despite the variety of models presented above, there are only three models with an S 3symmetric scalar potential which is not softly broken and with a realistic Yukawa sector, that could result in a viable DM candidate. These models are R-I-1 (0, 0, w S ), which was covered in Refs. [81,83] assuming a specific limit, R-II-1a (0, w 2 , w S ), and C-III-a (0,ŵ 2 e iσ 2 ,ŵ S ). Further on, we focus on the R-II-1a model and show that it could result in a viable DM candidate. The C-III-a model, which has spontaneous CP violation, will be presented elsewhere. 4 The R-II-1a model 4

.1 Generalities
The R-II-1a vacuum is defined by [6] {0, w 2 , w S }, (4.1) and the minimisation conditions are: The Z 2 symmetry is preserved for: Hence, the inert doublet is associated with h 1 , as h 1 = 0. A trivial Yukawa sector is assumed, L Y ∼ 1 f ⊗ 1 h , and thus the S 3 singlet is solely responsible for masses of fermions, making w S a reference point. Therefore, we define the Higgs-basis rotation angle as: After a suitable rephasing of the scalar doublets we chose w S > 0. With w 2 possibly negative, the Higgs basis rotation angle will be in the range β ∈ [− π 2 , π 2 ]. Therefore, the vevs can be parameterised as: The Higgs basis rotation is given by: (4.8)

R-II-1a masses 4.2.1 Charged mass-squared matrix
The charged mass-squared matrix in the {h + 1 , h + 2 , h + S } basis is given by: The charged mass-squared matrix is diagonalisable by the rotation (4.7). Therefore, the mass eigenstates can be expressed as: with masses: Positivity of the squared masses requires the following constraints to be satisfied:

Inert-sector neutral mass-squared matrix
The mass terms of the neutral components of the h 1 doublet are already diagonal. The masses of the two neutral states are given by: Positivity of the masses squared requires the following constraints to be satisfied:

Non-inert-sector neutral mass-squared matrix
The neutral mass-squared matrix is block-diagonal in the basis {η 2 , η S , χ 2 , χ S }. Therefore the mass-squared matrix can be split into two blocks: where "2S" refers to the mixing of h 2 and h S . The mass-squared matrix of the CP-odd sector is: It can be diagonalised by performing the R β rotation (4.7). The two CP-odd states are: Positivity of the squared masses requires The CP-even mass-squared matrix is: The mass-squared matrix is diagonalisable by The CP-even states are: with masses: (4.28) We identify the lighter state, h, as the SM-like Higgs boson.
It should be noted that we identified the mass eigenstates of the CP-even sector without performing a rotation to the Higgs basis. Had we done that according to the CP-even sector would have been diagonalised by the additional rotation (4.31)

Mass eigenstates
In terms of the mass eigenstates, the SU(2) doublets can be written as: whereas in the Higgs basis the SU(2) doublets can be written as: . (4.33c) The expressions for the squared masses can be inverted to yield the scalar potential couplings: The model is invariant under a simultaneous transformation of which could have been adopted, but is not, in order to reduce ranges of parameters.

The R-II-1a couplings
Below, we quote the gauge and Yukawa couplings of the R-II-1a model. The scalar-sector couplings are collected in appendix A.

Gauge couplings
After substituting the doublets in terms of the mass eigenstates (4.32) into the kinetic Lagrangian, the resulting terms are: where we have left out couplings involving Goldstone fields. From the interaction terms ZZh and ZZH it follows that the states h and H are CP-even and therefore the state A is CP-odd. Provided that the h scalar is associated with the SM-like Higgs boson, from the interactions hZZ and hW ± W ∓ it follows that the SM-like limit is reached for sin(α + β) = 1. (4.37)

Yukawa couplings
As noted earlier, there are two possibilities to construct the Yukawa Lagrangian: Although the first option can give realistic fermion masses, the CKM matrix splits into a block-diagonal form. We consider the trivial representation for fermions 6 : where the Yukawa couplings y are assumed to be real, as mentioned earlier, andh S is the charge conjugated of h S , i.e.,h S = iσ 2 h * S . The fermion mass matrices are to be transformed as the superscripts "0" on the fermion fields indicate a weak-basis field.
When considering the trivial Yukawa sector, the CKM matrix, V CKM = V † u V d , can be easily fixed to match the experimental value. Moreover, there is no naturally, at tree-level, occurring source responsible for FCNC. The scalar-fermion couplings can be extracted from eq. (4.38) by transforming into the fermion mass-eigenstates basis and multiplying the appropriate coefficients, next to the fields, by −i: and for the leptonic sector, the Dirac mass terms would lead to similar relations. The SM-like limit for the scalar h, Finally, the charged scalar-fermion couplings are: The structure of the Yukawa couplings is the same as that of the 2HDM, Type I, except that our definition of tan β is the inverse, since the singlet vev is here taken as the reference (denominator): Note also that α is defined differently.

Model analysis
For simplicity, we introduce a generic notation for different scalars, Precise measurements of the W ± and Z widths at LEP [92] forbid decays of the gauge bosons into a pair of scalars. In our case, the lower limits on the scalar masses are set by the following constraints: m ϕ ± i > 1 2 m Z , and m ϕ i + m h ± > m W ± , and m η + m χ > m Z . Usually, a conservative lower bound for the charged masses m ϕ ± i ≥ 80 GeV is adopted [93,94]. We assume a slightly more generous lower bound of m ϕ ± i ≥ 70 GeV. The model is analysed using the following input: • Mass of the SM-like Higgs is fixed at m h = 125.25 GeV [95]; • The Higgs basis rotation angle β ∈ [− π 2 , π 2 ] and the h -H diagonalisation angle α ∈ [0, π]; • The charged scalar masses m ϕ ± i ∈ [0.07, 1] TeV; • The inert sector masses m ϕ i ∈ [0, 1] TeV. Either η or χ could be a DM candidate, whichever is lighter; • The active sector masses {m H , m A } ∈ [m h , 1 TeV]; For the numerical parameter scan, both theoretical and experimental constraints are evaluated. Based on the constraints, several cuts are defined and applied: • Cut 1: perturbativity, stability, unitarity checks, LEP constraints; • Cut 2: SM-like gauge and Yukawa sector, electroweak precision observables and B physics; • Cut 3: h → {invisible, γγ} decays, DM relic density, direct searches; with each of the subsequent constraint being superimposed over the previous ones. Imposing the theory constraints discussed in appendix B, we can exclude parts of the parameter space, as illustrated in figure 2. Low values of w 2 , and hence of β (see eq. (4.6)), are disfavoured by the R-II-1a model. This can be seen by inspecting the λ 1 , λ 2 , and λ 3 couplings (4.34), these couplings are proportional to 1/w 2 2 . A particularly instructive combination, expanded for small β, is

Imposing theory constraints
with λ 1 minimised for α = π/2. With a decreasing denominator ∼ β 2 , we want to control the overall value of |λ i |, and therefore the value of the numerator must also decrease. The v cos β , and the perturbativity requirement, |Y f | ≤ 4π, it follows that the most stringent bound comes from the heaviest state, which is m f = m t . For these values we find that | tan β| ≤ 12.6. While there would be Landau poles in the vicinity of this non-perturbative region, we note that other constraints prevent parameters from getting close.
Finally, we shall assume that the unitarity constraint is satisfied for the value of 16π rather than 8π, see appendix B.2. It turns out that points surviving all of the checks satisfy unitarity constraint at the value of 8π. Figure 3: Constraints on α and β from the gauge and Yukawa couplings. We take β ∈ [− π 2 , π 2 ] and α ∈ [0, π], see also eq. (4.35). The white region is excluded. The black diagonal line identifies the SM-like limit, which is α + β = π/2. The coloured regions are compliant with Cut 1, whereas the yellow regions are also compliant with Cut 2 constraints. Values of β for which w 2 or w S vanish are identified.

The SM-like limit
The Standard Model like limit refers to the limit in which the CP even boson, which is in the only doublet that acquires vev when going to the Higgs basis, is already a physical boson, i.e., a mass eigenstate and therefore, does not mix with the other neutral scalars. This doublet is the one that contains the would-be Goldstone bosons G ± and G 0 . This limit is special and is referred to as the SM-like limit because this CP even neutral boson behaves in many aspects as the SM boson. This limit is reached for sin(α + β) = 1 (α = −β + π/2) as stated before. In this case we see from eq. (4.36a) that only h has couplings of the type V V H and their strengths coincide with those of the SM. Furthermore, from the fact that h is in the only doublet that acquires vev we see that in this limit h couples to the fermions with the same strength as the SM Higgs boson. In fact, this can be seen from eq. (4.39), since in this limit sin α = cos β.
We shall adopt the following 3-σ bounds from the PDG [95]: Note that we must impose the same sign for these two couplings of eq. (5.3), in order not to spoil the interference required for h SM → γγ.

Electroweak precision observables
The electroweak oblique parameters are specified by the S, T , and U functions. Sufficient mass splittings of the extended electroweak sector can account for a non-negligible contribution. The S and T parameters get the most sizeable contributions. Results are compared against the experimental constraints provided by the PDG [95], assuming that U = 0. The model-dependent rotation matrices, needed to evaluate the set of S and T , are presented in appendix C.2.

B physics constraints
The importance of a charged scalar exchange for theB → X(s)γ rate has been known since the late 1980's [96][97][98]. The rate is determined from an expansion of the relevant Wilson coefficients in powers of α s /(4π), starting with (1) the matching of these coefficients to the full theory at some high scale (µ 0 ∼ m W or m t ), then (2) evolving them down to the low scale µ b ∼ m b (in this process the operators mix), and (3) determine the matrix elements at the low scale [99][100][101][102][103][104][105][106][107][108][109][110][111][112]. We here follow the approach of Misiak and Steinhauser [113]. While the considered S 3based models have two charged Higgs bosons, only one couples to fermions. This implies that we may adopt the approach used for the 2HDM with relative Yukawa couplings of the active charged scalar, eq. (4.41) (in the notation of Ref. [113]), We show in figure 4 the regions in the tan β -m H + parameter plane that are not excluded by this constraint. The situation is quite different from that of the more familiar 2HDM with Type II Yukawa couplings. According to eq. (5.4) the relevant couplings are the same as those of the 2HDM Type I model, with the exception that we are here interested in small values of tan β. While theB → X(s)γ constraint excludes large values of tan β, low values are for R-II-1a also cut off due to the perturbativity constraint, we have | tan β| > ∼ 0.26, as discussed in section 5.1.
Since A u A d > 0, there is a region of negative interference between the SM-type contributions and the loop with the charged Higgs: as we increase the value of tan β for fixed m H + , the branching ratio will first diminish, and then at some point come back up, as illustrated in the lower right-hand corner of the left panel of figure 4. This interference region is ruled out by Cut 3. For any fixed value of tan β, on the other hand, for sufficiently high mass m H + , the rate approaches the SM value.

LHC Higgs constraints
We require that the Higgs-like particle, h, full width is within Γ h = 3.2 +2.8 −2.2 MeV, which is an experimental bound adopted from [95]. In the SM the total width of the Higgs boson is around 4 MeV. The upper value, i.e., Γ h = 6 MeV is used in preliminary checks within the spectrum generator.

Decays h → γγ and h → gg
The di-photon partial decay width is modified by the charged-scalar loop in comparison to the SM case. The one-loop width is known [2,114,115]: where α is the fine-structure constant, Q f is the electric charge of the fermion, N c = 3 (1) for quarks (leptons), and the C's are the couplings normalised to those of the SM, The spin-dependent functions F S 1/2 , F 1 , and F 0 can be found in appendix C. Contributions from the fermionic loop and gauge loop versus the charged scalar loop are presented in figure 5. Figure 5: Scatter plots of additional contributions to the di-photon decay amplitudes, normalised to the SM value, expressed in per cent. Complex numbers arise from the spin-dependent function F S 1/2 of equation (C.1b) for fermions lighter than the Higgs-like particle h. The A S value represents the normalised contribution from the charged scalar loop, A S /|A SM |, while the ∆A W+F value stands for an additional contribution to the SMlike part due to the W and fermion loops, ∆A W+F = (|A W+F | − |A SM |) /|A SM |. Yellow dots represent parameters surviving Cut 2, while the red ones represent those surviving also LHC Higgs-particle constraints: full width, the invisible branching ratio, and the di-photon constraint.
In the SM case, the dominant Higgs production mechanism is through gluon fusion. However, due to experimental limitations we do not explicitly consider constraints on this channel. The rate for the two-gluon decay at the leading order is [116][117][118][119] where α S is the strong coupling constant. The decay width of this process can be enhanced or diminished with respect to the SM case. Such behaviour is caused by an additional factor for the amplitude, g R-II-1ā The normalised two-gluon branching ratio to the SM case is depicted in figure 6. The gluon branching ratio for DM mass below m h /2 can become low due to the opening of the invisible channel, h → ϕ i ϕ i . However, such cases are partially excluded by other LHC Higgs-particle constraints of Cut 3. Even with more experimental data collected, the two-gluon constraint will not play a very significant role in terms of constraining the R-II-1a model. Most of the two-gluon points are within the range µ gg ∈ [0.9, 1.1]. Figure 6: Two-gluon versus di-photon Higgs-like particle branching ratios, normalised to the SM value, R ii ≡ Br(h → ii)/Br SM (h → ii), for different R-II-1a mass orderings. The black lines represent the SM case. A proper evaluation of the di-photon channel (5.9) would require an additional factor from the two-gluon channel arising from the h production, i.e., each point would have to be scaled by multiplying both of the values presented in the plot, µ γγ = R gg R γγ . Points above the R gg = 1 line would shift right, while those below would shift left.
In light of the above discussion, we do not aim to account for the correct h two-gluon production factor and approximate the di-photon channel strength to be with µ γγ = 1.11 ± 0.10 [95]. We evaluate this constraint allowing for an additional 10 per cent computational uncertainty, and impose an (n = 3)-σ tolerance, The decay width of h into a pair of scalars ϕ i is given by with a symmetry factor (2 − δ ij ), where δ ij is the Kronecker delta. The appropriate trilinear couplings are given by equations (A.2a) and (A.2c). In appendix A the overall factor of "−i" is left out 7 . Due to CP conservation there is no coupling g (h η χ), therefore the Higgs-like particle can decay only into pairs of inert neutral scalars, These are the only two possible channels (h → ηη and h → χχ) which can contribute to the invisible decay. Channels with other scalars are kinematically inaccessible due to the assumed limits, i.e., mass ordering and LEP constraints. The experimental bounds can be applied directly if there is only a single channel open. However, there is a possibility that both inert neutral states, η and χ, can be kinematically accessible. In a simplistic approximation a particle escapes a detector of 30 meters, assuming no time dilation factor, if its lifetime exceeds a value of τ 10 -7 seconds. The value can be expressed in terms of the total decay width, Γ tot 6.6 × 10 -18 GeV. Based on the total width of the next-to-lightest state, ϕ j , two situations are possible: • If the particle is not long-lived, Γ tot (ϕ j ) > 6.6 × 10 -18 GeV, it will decay within the detector through h → ϕ j ϕ j → ϕ i ϕ i Z * Z * , with the Z * subsequently also decaying.
In this case only the h → ϕ i ϕ i channel will contribute to Br (h → inv.).
• When Γ tot (ϕ j ) < 6.6 × 10 -18 GeV, the invisible branching ratio will be given by (5.14) In the R-II-1a model we found that the decay rate of the next-to-lightest inert neutral particle is way above Γ = O(10 -18 ) GeV. In our calculations we shall adopt the PDG [95] constraint, which is Br exp (h → inv.) < 0.19.

The h scalar self interactions
Let us next consider the trilinear and quadrilinear self interactions, given by eqs. (A.1a) and (A.4b). The SM Higgs self-interactions are [122] In the R-II-1a model, these couplings can be expanded in terms of m 2 h , m 2 H and m 2 η ,
In the future, the trilinear Higgs self interactions may become a crucial test for new physics. For this purpose, we show in figure 7 the h trilinear coupling, relative to the SM value.

Astrophysical Observables
We consider a standard cosmological model with a freeze-out scenario. The cold dark matter relic density along with the decay widths discussed above and other astrophysical observables are evaluated using micrOMEGAs 5.2.7. The 't Hooft-Feynman gauge is adopted, and switches are set to default values VZdecay = VWdecay = 1, identifying that 3-body final states will be computed for annihilation processes only. The fast = −1 switch identifies that very accurate calculation is used. The steering CalcHEP [123] model files are produced with the help of SARAH [124,125].
Let us recall results of figure 1. In the IDM, two regions compatible with the cold dark matter relic density were identified. These correspond to the high-mass region and the intermediate-mass region. The high-mass region is in agreement with the relic density due to two factors: • Near-mass-degeneracy among the scalars of the inert sector. Small mass splittings correspond to tiny couplings, and different inert-scalar contributions to the annihilation will be suppressed and result in an acceptable relic density; • Freedom to choose the Higgs boson portal coupling λ L . This parameter controls the trilinear XXh and quartic XXhh coupling, and must be sufficiently small. Here, X refers to scalars of the inert sector, both charged and neutral.
In this region, the main DM annihilation is into W + W − . However, the relic density can be maintained at an acceptable level by suppressing annihilation via an intermediate h boson and into a pair of h bosons, illustrated in figure 8. The desired relic abundance can be achieved by adjusting the mass splittings. Whereas the IDM and the 3HDMs considered in figure 1 have an adjustable portal coupling (often referred to as λ L for the IDM), the present model is constrained by the underlying S 3 symmetry. Here, there is not a single portal coupling, but two: a trilinear XXh and a quartic XXhh, plus additional ones involving other scalars. Furthermore, these are not "free", but correlated with other features of the model. In particular, they are constrained by the scalar masses and two angular parameters, α and β.
Let us consider a simplified picture with heavy active scalar bosons. At high DM masses the main DM annihilation channel is into W + W − . This channel is controlled by the gauge coupling. In addition, channels leading to h scalars will be accessible, as illustrated in figure 8. These amplitudes are controlled by the portal couplings which should be constrained, since otherwise the DM relic density becomes too low.
To first order in δ (in the neighbourhood of the SM-like limit), where α = π/2 − β + δ, the behaviour of the portal couplings is quite simple. In the limit of δ → 0 we arrive at This relation shows that the portal couplings will grow with increasing DM mass. Values of the trilinear and quartic couplings are shown in figure 9. As shown in this figure, the correlation with DM mass (5.18) is qualitatively borne out by the parameter points surviving Cut 2.
In the aforementioned simplification we argued that there are no good DM candidates for high mass values. In the full R-II-1a model other annihilation channels involving active scalars are also accessible. For example, there is a contribution from the H scalar to the X X → h h process through the s-channel. Furthermore, there are other accessible active-inert scalar channels. This explains why the DM relic density is saturated at Ωh 2 = O(10 −4 ) at high DM values, see figure 10. Figure 9: Absolute value of the trilinear portal coupling |g(XXh)/v| (top) and the quartic portal coupling |g(XXhh)| (bottom) versus the lightest inert particle mass. Figure 10: Dark matter relic density for the R-II-1a model. The region compatible with the observed DM relic density does not allow for masses above around 120 GeV. In the high-mass region, m ϕ i > ∼ 500 GeV, the DM relic density is shown to be too low.
Going down in the DM mass, the situation changes as follows. At around m ϕ i ≈ 200 GeV there is a kink and the maximal relic density increases from Ωh 2 = O(10 −4 ) up to Ωh 2 = O(10 −2 ). In this range, the main annihilation (or loss) mechanisms are via the channels ϕ i ϕ i → hh and ϕ i ϕ i → W ± W ∓ . In this same mass region many parameter points also yield Ωh 2 < O(10 −4 ). This happens when the dominant annihilation channels are ϕ i ϕ i → AZ and ϕ i ϕ i → H ± W ∓ , through either h or H in the s-channel, or via ϕ j or h ± (based on quantum numbers) in the t-channel.
Then, in the mass region m ϕ i ∈ [m h /2, 120 GeV], the relic density ranges from above unity down to below Ωh 2 = O (10 −5 ). The main annihilation channels are into a pair of (virtual) W ± bosons or b-quarks, with the latter becoming increasingly important as the DM mass decreases. However, there are also cases when the dominant annihilation channel is ϕ i ϕ i → gg, which can contribute more than 50%.
Finally, in the DM mass region corresponding to values below m h /2, the primary annihilation or loss mechanism is ϕ i ϕ i → bb trough a virtual h. This channel depends critically on the portal, i.e., the trilinear coupling g(ϕ i ϕ i h).

The η case
For the case when the η scalar is the lightest ("η case"), from figure 10 it looks as if there could be solutions for the range of m η ∈ [2, 120] GeV. However, in the lowmass range of this interval, the h invisible branching ratio, together with the relic DM density constraint becomes incompatible with the experimental data. In light of this fact, in the remainder of the discussion presented in this paragraph we shall focus on η masses up to 120 GeV since we already know that criterion (a ) excludes higher masses. Both constraints, i.e., (a ) and (b ), alongside with Cut 1 and Cut 2, are satisfied within the region m η ∈ [40,120] GeV. The final checks are then the di-photon (c ) and the direct detection constraints (d ). These two constraints are very severe and eliminate a large region of the parameter space. When imposed separately, the strongest constraint comes from the direct detection criteria, which are satisfied in the mass region m η ∈ [43,120] GeV and at values below m η 10 GeV.
We list different Cut 3 paired constraints in figure 11 after imposing cuts 1 and 2. We work with 8 input parameters, 6 masses and 2 angles. After imposing Cut 3, there will be different domains allowed by each of the three checks, Ωh 2 , DD, or LHC, imposed separately. The intersection of all these domains would correspond to Cut 3 being satisfied. Figure 11 shows the allowed mass regions of the DM candidate after imposing two of the different checks at a time. Overlapping lines do not guarantee that there are regions of parameters satisfying all of the constraints simultaneously since there are seven more parameters to consider. There is no region for either (Ωh 2 +DD) or (LHC+DD) satisfied for m η ≈ m h /2. In fact, for the η case we found no parameter point satisfying all of the Cut 3 constraints simultaneously. Figure 11: Allowed mass regions of the DM candidate involving different Cut 3 constraints. Blue: relic density satisfied together with direct detection. Purple: LHC Higgs constraints along with direct detection constraints. Red: relic density and LHC Higgs constraints. Grey: all of the Cut 3 constraints are satisfied. Note that additional input parameters are not shown.

The χ case
For the case when the χ scalar is the lightest ("χ case"), the Ωh 2 distribution is slightly shifted towards lower relic density values, as shown in figure 10. As a result, the region compatible with the relic density is m χ ∈ [2, 105] GeV. For the η case with masses below 40 GeV, when the (b ) constraint is imposed, all Ωh 2 are above 0.22. This does not apply to the χ case since in this case Ωh 2 can go as low as ≈ 0.07. Nevertheless, the sub-40 GeV region is not compatible with Cut 3. However, the region m χ ∈ [52.5, 89] GeV survives cuts 1 to 3 when applied simultaneously. The lower DM mass range is compatible with other models presented in figure 1, while slightly heavier DM candidates are also allowed within the R-II-1a framework. When applied simultaneously, conditions (a ) and (d ) are satisfied for a broader range of m χ ∈ [45.5, 92] GeV. Bounds from pairs of different Cut 3 checks are shown in figure 11.
It is instructive to see which points within the parameter range survive all the cuts. The mass scatter plots of Cut 3 superimposed on the Cut 2 points are given in figure 12. There are no solutions with active neutral states being degenerate. Moreover, these masses reach at most m max A ∼ 340 GeV and m max H ∼ 450 GeV. The CP-odd state, A, can be as light as the observed SM-like Higgs boson, m A m h . On the other hand, such low masses for the H boson are disfavoured by Cut 3. The charged bosons that survive Cut 3 are also light, with m max H ± ∼ 460 GeV and m max h ± ∼ 310 GeV. It is interesting to note that the charged active scalar can be as light as m min H ± ∼ 177 GeV. Future experimental data on the decays of the charged bosons will be useful to test the model. For the B physics constraints we applied only the indirect experimental bounds on theB → X(s)γ rate, however other channels might lead to stronger constraints on the mass of the H ± scalar. Finally, as noted earlier we have for the DM candidate m χ ∈ [52.5, 89] GeV and the other inert neutral state, η, can be as heavy as m max η ∼ 310 GeV. It is also interesting to note that either η or h ± can be the next-lightest member of the inert doublet.
In addition to the mass parameters, angles are also used as input. The allowed ranges in the α -β plane are shown in figure 13. The solutions satisfying all cuts are asymmetrically distributed along the black diagonal line, which represents the SM-like limit. In the right panel of figure 13 we also show the gauge (κ V V ) and Yukawa (κ f f ) couplings, relative to the SM values, see eqs. (5.3). Figure 13 shows that the experimental data are more stringent for κ f f than for κ V V . The lack of symmetry in the distribution of the grey points with respect to the diagonal, corrresponding to the SM-like limit, in the left panel translates into a significant population of κ f f -values below unity in the right panel. Figure 12: Scatter plots of masses that satisfy different cuts, for both orderings of η and χ. Identical notation as in figure 2. The light-blue region satisfies Cut 1 and accommodates the 16π unitarity constraint. The yellow region accommodates a 3-σ tolerance with respect to Cut 2, whereas in the green regions, the model is within the 2-σ bound of these values. The grey points are compatible with all cuts and are only present in the right-hand panels.
We present direct detection constraints in figure 14. We note that in practically the whole mass range there are parameter points at lower cross sections. Thus, a future improvement on this direct detection constraint is not obviously going to reduce the range of masses allowed by the model.
Finally, in table 2 we show some benchmarks.  The spin-independent DM-nucleon cross section compatible with XENON1T [77] data at 90% C.L. The points represent Cut 3 satisfied for the χ DM case. The red line corresponds to an approximate neutrino floor, which can be defined in various ways. Table 2: Benchmark points and dominant decay modes. The "q" notation refers to a sum over the light quarks, d, u, s and c, "l" refers to all leptons, and "ν" to all neutrinos.

Concluding remarks
It is possible to incorporate a DM candidate within the S 3 -symmetric scalar model. The DM candidate requires one of the Higgs doublets to be inert. As S 3 symmetry is assumed and a DM candidate is sought, this requirement imposes constraints on the structure of the Yukawa Lagrangian. As a matter of fact, we found no possible combinations of the exact S 3 -symmetric scalar potential and a non-trivial Yukawa Lagrangian, which could accommodate a DM candidate. Hence we require fermions to transform trivially under S 3 . When soft-symmetry breaking terms are present, it is possible to construct the Yukawa Lagrangian with a non-trivial S 3 structure. Due to such behaviour an ad-hoc S 3 is not appealing in the context of being simultaneously applied to the scalar sector and generating a non-trivial Yukawa sector, while trying to explain DM.
In this work we focused on a specific S 3 -symmetric scalar model R-II-1a. As was shown in the paper, this model is compatible with several constraints, both theoretical and experimental. In the IDM and the literature on 3HDMs there is a viable DM high-mass region present. This is not true in our case. The main difference is that the inert-active scalar couplings in the R-II-1a are constrained by the underlying S 3 symmetry and hence the portal couplings are harder to adjust.
We analysed the R-II-1a model numerically. Within the model there are two possible DM candidates present, η and χ. After performing the analysis we found no solutions satisfying all of the constraints with η being the lightest. However, for the χ case there is a range compatible with all constraints, m χ ∈ [52.5, 89] GeV. Constraints in this mass range are compatible with data at the 3-σ level. As compared with the IDM, the R-II-1a model allows solutions in the intermediate-mass range up to somewhat higher values, but can not satisfy all constraints in the high-mass region.
The model differs from the IDM in having two non-inert doublets. The corresponding scalars must be rather light, as shown in figure 12. If these were to be produced at the LHC, they could decay to a pair of scalars from the inert doublet, as well as to final states familiar from the 2HDM. Furthermore, one could imagine direct production of two scalars of the inert doublet, for example via W W or W Z fusion, at rates controlled by the gauge couplings. For all of these cases, the charged h ± and neutral η would decay via gauge boson emission to the DM, h ± → W ± χ and η → Zχ, leading to mono-vector events [126,127].
The present model looks very promising. There are some similarities between this model and the 2HDM. The present framework preserves CP both at the Lagrangian level and by the vacuum. The presence of both g (ηh ± H ∓ ) and g (χh ± H ∓ ) couplings suggests there might be mixing at the one-loop level, but the two diagrams associated with the different charge assignments cancel. Since there is no CP violation in the scalar sector, the Yukawa couplings must be complex and CP is only violated via the CKM matrix. The fermions only couple to one of the active scalars, the one that is a singlet under S 3 since the fermions also transform trivially under S 3 . In the S 3 -symmetric 3HDM there are also regions of parameter space leading to vacua that violate CP spontaneously, as discussed in Refs. [6,86], showing that the S 3 -symmetric 3HDM has a very rich structure.

Acknowledgements
It is a pleasure to thank Igor Ivanov, Mikolaj Misiak and Alexander Pukhov for very useful discussions. We thank the referee for pointing out a mistake in the earlier version. PO is supported in part by the Research Council of Norway. The work of AK and MNR was partially supported by Fundação para a Ciência e a Tecnologia (FCT, Portugal) through the projects CFTP-FCT Unit UIDB/00777/2020 and UIDP/00777/2020,

A Scalar-scalar couplings
For simplicity, the scalar-scalar couplings are presented with the symmetry factor, but without the overall coefficient "−i". We denote the "correct" couplings as g ... = −ig (. . .). We shall abbreviate c θ ≡ cos θ, and s θ ≡ sin θ, and t θ ≡ tan θ for any argument θ.
The trilinear scalar-scalar couplings involving the same species are: The trilinear couplings involving the neutral fields are: The trilinear couplings involving the charged fields are: The quartic couplings involving the same species are: The quartic couplings involving only the neutral fields are: The quartic couplings involving both neutral and charged fields are: The quartic couplings involving only the charged fields are:

B Theory constraints
We impose certain data-independent theory constraints on the models.

B.1 Stability
Necessary, but not sufficient, conditions for the stability of an S 3 -symmetric 3HDM were provided in Ref. [89]. In Ref. [6], based on the approach of Refs. [78,128], necessary and sufficient conditions for models with λ 4 = 0 were discussed. It was later pointed out in Ref. [129] that parameterisation used in Ref. [78], and hence in Ref. [6], is not correct 8 , and one would arrive at a value of the potential which would be lower than what actually is possible to achieve within the space available.

B.2 Unitarity
The tree-level unitarity conditions for the S 3 -symemtric 3HDM were presented in Ref. [89]. The unitarity limit is evaluated enforcing the absolute values of the eigenvalues Λ i of the Table 3: Reproduction of the necessary stability conditions of Ref. [89] in terms of the parameterisation given by (B.2).

B.3 Perturbativity
The perturbativity check is split into two parts: couplings are assumed to be within the limit |λ i | ≤ 4π, the overall strength of the quartic scalar-scalar interactions is limited by |g ϕ i ϕ j ϕ k ϕ l | ≤ 4π. For the R-II-1a model, the list of the quartic scalar interactions g ϕ i ϕ j ϕ k ϕ l can be found in appendix A. From the interactions ηηηη and χχχχ (A.4a), it follows that 0 < λ 1 +λ 3 ≤ 2π/3. Evaluation of other couplings is more involved.  and