Dark matter in three-Higgs-doublet models with S3 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 S3 symmetry. These models have many distinct vacua with one or two vanishing vacuum expectation values which can be stabilised by a remnant of the S3 symmetry which survived spontaneous symmetry breaking. We discuss all possible vacua in the context of S3-symmetric three-Higgs-doublet models, allowing also for softly broken S3, 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 S3 symmetry and is such that the flavour structure of the model is solely governed by the VCKM matrix which, in our framework, is not constrained by the S3 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.


JHEP01(2022)120
The micrOMEGAs code [10][11][12] is used at this stage to evaluate the cold dark matter relic density along with the decay widths and other astrophysical observables. The full impact of Cut 3 is discussed in section 6, where we also present a set of benchmarks for this model. Two different scenarios emerge right from the beginning of the numerical analysis presented in section 5, based on the ordering of the masses of the two neutral inert scalars, η and χ, which have opposite CP parities. The origin of the masses of these two fields is also different in terms of parameters of the potential. Cut 3 removes the possibility of having the η scalar lighter than χ. As a result, in this model only χ can play the rôle of DM. We also conclude that there are no good DM candidates for high mass values as is explained in our discussion. In section 7 we present our conclusions. There are several appendices where features of the model and details of our analysis are given.

Scalar dark matter
One of the simplest extensions of the SM which could accommodate DM is obtained by adding a scalar singlet. With an explicit Z 2 symmetry, to prevent specific decay channels, this extension could yield a viable DM candidate [13][14][15][16][17][18]. Direct detection constraints were studied [19][20][21][22] based on the LHC data. The parameter space of these models was further constrained after the Higgs boson discovery [23][24][25][26][27][28][29]. While being the subject of specific exclusion criteria, in general, two DM regions were identified: a low-mass region 55-63 GeV and a high-mass region above around 100-500 GeV, depending on the implementation. In alternative, the scalar singlet DM can be stabilised by a different symmetry, such as Z 3 , which was considered in refs. [30,31].
A DM candidate can also be introduced through non-trivial SU(2) L scalar n-tuples. In ref. [32] a class of so-called minimal dark matter (MDM) multiplets was proposed. The key aspect of these models is to extend the SM by additional scalar multiplets with minimal quantum numbers (spin, isospin, hypercharge) to accommodate the DM candidate. The MDM models were further on studied in refs. [33][34][35][36][37][38][39][40].

The Inert Doublet Model
One of the most popular models, capable of accommodating DM, is the IDM [4,5]. The IDM accommodates DM in the form of a neutral scalar: the lightest neutral member of an inert SU (2) Higgs doublet added to the SM. In addition to providing an economical accommodation of DM, the IDM also offers a mechanism for generating neutrino masses [41][42][43]. The IDM was studied extensively some ten years ago [36,[44][45][46][47][48][49][50][51][52][53][54][55][56][57][58]. Assuming the Higgs mass to be around 120 GeV, two DM mass regions were identified: a low and intermediate-mass region, from about 5 GeV to around 100 GeV (allowing for a heavy Higgs up to 500 GeV, this DM mass region would extend up to 160 GeV), and a high-mass region, beyond about 535 GeV. Above some 80 GeV, the annihilation in the early Universe to two gauge bosons (W + W − or ZZ) becomes very fast and the relic DM density would be too low. In addition, annihilation into two SM-like Higgs bosons is possible and is controlled by the effective λ L coupling (parameterising the Higgs-DM-DM coupling), for small λ L values the JHEP01(2022)120 overall effect is negligible. Eventually, for sufficiently heavy DM particles, above approximately 500-535 GeV, the annihilation rate drops and as a result the DM density is again compatible with the data.
When the mass splitting between the DM candidate and the other scalars in the inert sector is sufficiently small, coannihilation effects get stronger. As a result, if the Higgs-DM-DM coupling, λ L , has a suitable value, a relic density in agreement with observation can be achieved. For a fixed mass splitting with an increasing DM mass the absolute value of the λ L coupling should increase to satisfy the experimental relic density value. On the other hand, for a fixed λ L value, if the mass splittings are increased, the relic density decreases. The relic density depends also on the sign of λ L . Interplay of these parameters may result in an acceptable relic density for the DM candidate of several TeV.
With the Higgs boson discovery [8,9], it was pointed out that the decay width of the Higgs boson to two photons and invisible particles could further constrain the IDM [59][60][61][62]. It was found [63,64] that with both astroparticle and collider constraints taken into account the DM mass region below 45-50 GeV is ruled out, and thus a minimal mass threshold is set. Later, the whole low-and intermediate-mass region was shown to be 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.
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.
• 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].

JHEP01(2022)120
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.

JHEP01(2022)120
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 L Y ∼ (2 ⊕ 1) f ⊗ (2 ⊕ 1) h , one doublet and one singlet of S 3 , 2 : Such structure yields the mass matrix for each quark sector (d and u) of the form

JHEP01(2022)120
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   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]. 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).
The Yukawa sector is unrealistic.
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: • R-II-1a: (0, w 2 , w S ).
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.
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 JHEP01(2022)120 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 .
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.
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 .
Both the scalar and Yukawa sectors are unrealistic.
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.
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 This case might result in a viable DM candidate provided that the S 3 symmetry of the potential is softly broken.

JHEP01(2022)120
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 .
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 non-negligible 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.

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)

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: where M 2 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.

JHEP01(2022)120
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:

JHEP01(2022)120
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: 1 1 In our study neutrino masses are of no particular interest.

JHEP01(2022)120
where the Yukawa couplings y (u, d) ij 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): (tan β) R-II-1a = 1 tan β 2HDM, Type I . (4.42) 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 JHEP01(2022)120 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; 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 theory constraints
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 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 perturbativity constraint restricts large values of the λ i couplings, see appendix B.3, and sets a limit 0 < λ 1 + λ 3 ≤ 2π/3 (A.4a). For m 2 η 3m 2 h , we arrive at the bound In our model fermions couple only to the h S doublet. This can be used to put a limit on the tan β value. From the definition of the Yukawa couplings, Y f = √ 2m f 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π. . 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 JHEP01(2022)120 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 3 -based 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 JHEP01(2022)120 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.
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 twogluon 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].
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,

Invisible decays, h → inv.
The SM-like Higgs boson can decay to the lighter scalars h → ϕ i ϕ j provided 2m ϕ i < m h . If the decays are kinematically allowed, these processes can enhance the total width of the SM-like Higgs state sizeably. 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. 2 Due to CP conservation there is no coupling g (h η χ), therefore the Higgs-like particle can decay only into pairs of inert neutral scalars,

JHEP01(2022)120
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 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 η , After imposing α + β = π/2 we arrive at the SM-like Higgs couplings.
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].
We shall adopt the cold dark matter relic density value of 0.1200 ± 0.0012 taken from PDG [95]. The relic density parameter will be evaluated using a 3-σ tolerance and assuming an additional 10 per cent computational uncertainty, corrponding to the [0.1075; 0.1325] region. 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; Figure 8. Feynman diagrams contributing to XX (particles of the inert sector) annihilation channels at high DM masses.
• 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 activeinert scalar channels. This explains why the DM relic density is saturated at   . 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.
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 JHEP01(2022)120 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. 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 low-mass 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 JHEP01(2022)120 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.

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. 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.  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.

JHEP01(2022)120
Finally, in  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.

JHEP01(2022)120 7 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.

JHEP01(2022)120
Acknowledgments 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

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: + λ a c 3 α s β + 2c 2 α s α (3λ 4 s β + (3λ 8 − λ a ) c β ) , The trilinear couplings involving the charged fields are: one would arrive at a value of the potential which would be lower than what actually is possible to achieve within the space available. In our case, due to the freedom of the λ 4 coupling, which breaks the O(2) symmetry, the stability conditions are rather involved. We parameterise the SU(2) doublets as where the norms of the spinors ||h i || are parameterised in terms of the spherical coordinates
First, we check if the necessary stability constraints are satisfied, see table 3. Next, with the help of the Mathematica function NMinimize, using different algorithms, a further numerical minimisation of the potential is performed.

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 scattering matrix to be within a specific limit. In our scan we assume that one is given by the value |Λ i | ≤ 16π [130]. Some authors prefer a more severe bound |Λ i | ≤ 8π [131,132]. We compare the impact of both in figures 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.

C.2 V and U matrices
From refs. [133,134] we determine the V and U matrices 4 for R-II-1a in the Higgs basis (4.33)