Multi-Higgs models with CP symmetries of increasingly high order

When building CP-symmetric models beyond the Standard Model, one can impose CP-symmetry of higher order. This means that one needs to apply the CP-transformation more than two times to get the identity transformation, but still the model is perfectly CP-conserving. A multi-Higgs-doublet model based on CP-symmetry of order 4, dubbed CP4, was recently proposed and its phenomenology is being explored. Here, we show that the construction does not stop at CP4. We build examples of renormalizable multi-Higgs-doublet potentials which are symmetric under CP8 or CP16, without leading to any accidental symmetry. If the vacuum conserves CP-symmetry of order 2k, then the neutral scalars become CP-eigenstates, which are characterized not by CP-parities but by CP-charges defined modulo 2k. One or more lightest states can be the dark matter candidates, which are protected against decay not by the internal symmetry but by the exotic CP. We briefly discuss their mass spectra and interaction patterns for CP8 and CP16.


I. INTRODUCTION
The Standard Model (SM) is agnostic about the origin of the CP -violation which we observe in weak interactions [1]. SM simply postulates it and describes it via the complex Cabibbo-Kobayashi-Maskawa (CKM) matrix [2], but it provides no answer why the CP -symmetry should be broken at all. The search for a dynamical reason of why CP is broken is one of the motivations for building models beyond the SM (bSM), especially those with non-minimal Higgs sectors [3]. In fact, in the same year as Kobayashi and Maskawa put forward the idea that three quark generations can accommodate all CP -violating phenomena [2], T. D. Lee proposed in [4] the two-Higgs-doublet model (2HDM), where the CP -symmetry is broken spontaneously, as a result of the minimization of a CP -symmetric Higgs sector. At present we know that the CKM paradigm is indeed at work, while the Higgs boson properties revealed by the LHC are compatible with the SM Higgs [5]. Still, since the origin of the complex CKM matrix remains unexplained and since baryogenesis calls for yet additional sources of CP -violation, the intensive exploration of 2HDM [6] and more sophisticated multi-Higgs models [3,7] continues at full speed.
Investigation of how CP can be violated in multi-Higgs models recently led to a model with a peculiar form of CP -conservation [8]. Dubbed CP4 3HDM, this model is based on three Higgs doublets and incorporates a generalized CP -symmetry of order 4 denoted CP4. 1 Although it is known since long ago that models with several scalar fields φ i , i = 1, . . . , N with identical quantum numbers allow for unconventional definitions of CP -symmetry [9][10][11], in the vast majority of cases these definitions can be reduced to the standard one, with X ij = δ ij , by a basis change. For example, in 2HDM, one can define the CP -symmetry in the scalar sector in a variety of ways [12][13][14][15], but whatever definition one takes, the scalar sector of the model contains, in an appropriate basis, the conventional CP -symmetry [16]. 2 However, CP4 being a symmetry transformation of order 4, is markedly different. One needs to apply it four times, not twice, to obtain an identity transformation on fields. Thus, it cannot be reduced to the ordinary CP -symmetry by any basis change. In other words, the matrix X ij in (1) cannot be linked to δ ij by any basis change. This feature has clearly visible consequences in the scalar potential: despite the model is CP -conserving, it is impossible to find a basis in which all coefficients would be real. Technically, this is due to the existence of gauge-invariant non-hermitian combinations of Higgs fields which are invariant under CP4 instead of being mapped to their hermitian conjugates.
This property may also be linked to an interesting group-theoretical observation made in [18,19]. If one starts with a certain symmetry group G (which may include the Lorentz and gauge groups) and enlarges it with a CP -type symmetry, then this CP -transformation acts on G by an outer automorphism [10]. It turns out that the structure of the group G and the properties of its Clebsch-Gordan coefficients may influence this construction. In particular, for certain groups, this construction is possible but leads to a higher-order CPtransformation due to the complex Clebsch-Gordan coefficients. This offers another look at how irremovable complex coefficients may arise in CP -conserving models.
CP4 3HDM is the minimal multi-Higgs-doublet model whose scalar sector incorporates only CP4 without any accidental symmetries [8,20]. If CP4 is conserved by the minimum of the Higgs potential, then the model produces two mass-degenerate scalar dark matter (DM) candidates h and a. The model then resembles an enhanced version of the Inert doublet model [21][22][23][24][25], with two inert doublets and with the DM candidates stabilized not by the Z 2 -symmetry but by a CP -symmetry, albeit an unusual one. 3 Although the model is truly CP -conserving, one cannot classify h and a as being CP -even or CP -odd, as they transform under CP4 as h CP − − → −a and a CP − − → h. However, one can combine them into a single complex field ϕ, which then transforms under CP4 as The presence of the i factor and the absence of complex conjugation usually associated with a CPtransformation are highly peculiar and were discussed at length in [27].
With conserved CP4, one can quantify CP -properties of a field not by its CP -parity but by a global quantum number q defined modulo 4. One then assigns q = +1 to ϕ and q = −1 to its conjugate. In any transition between initial and final states with definite q, this quantum number is additively conserved modulo 4. When rewriting the inert self-interaction potential in terms of fields ϕ, terms such as ϕ 4 + (ϕ * ) 4 are allowed, since they also conserve the CP -charge q.
Can one go beyond CP4, while still keeping the interactions renormalizable? Can one build a multi-Higgs model invariant under a CP -symmetry of order 2k, where k > 2? On the one hand, one can certainly define generalized CP -symmetries of an arbitrary even order. However, if k contains any prime factor other than 2, one can split the group Z 2k into a pure family symmetry group and a group generated by a smaller-order CP -transformation. For example, since Z 6 ≃ Z 2 × Z 3 , imposing a CP -symmetry of order 6 would produce a model with a usual CP and a Z 3 family-symmetry group. The only way to prevent it is to take the order of the CP -symmetry 2k = 2 p , with integer p ≥ 1. The usual CP , which is of order two, can be denoted as CP2, the first non-trivial higher-order CP -symmetry is CP4, the next ones are CP8, CP16, and so on.
Next, although one can define CP8 or CP16 transformations in multi-Higgs models and impose them on the potential, it may easily happen that the model leads to accidental symmetries. For example, this is what happens in 3HDM [20]. Trying to impose CP8 leads to a model with an accidental continuous symmetry 2 It is worth mentioning that the scalar sector of 2HDM can accommodate a CP -transformation which cannot be transformed into the usual CP by any basis change [12][13][14][15]. This CP -transformation is defined by exactly the same matrix X as the one used in CP4 3HDM. However, within 2HDM, its effective order is not 4 but 2 due to the U (1) Y rephasing symmetry. This is best seen in the geometric picture where this unusual CP -transformation is described by a point reflection rather than plane reflection in the bilinear space, see detailed discussion in [14]. It is still a reflection, that is, a transformation of order 2, but it is different from the usual CP . The 2HDM based on this symmetry, which was dubbed in [14] the maximally CP -symmetric model, has a very peculiar phenomenology, especially when this symmetry is extended to the fermionic sector, [17]. It turns out, however, that the imposition of this maximal CP -symmetry entails other symmetries in 2HDM including the usual CP . The CP4 3HDM example is free from these accidental symmetries. 3 An example of models in which P -symmetry stabilizes a fermionic DM was presented in [26] and appeals to the known fact that a Majorana fermion picks up an i-factor upon P -transformation. In our case, CP4 stabilizes scalar DM candidates, and this phenomenon has a different origin.
U (1) and the usual CP , so that CP8 plays no special role in it. Since the classification performed in [20] was exhaustive, it means that one needs to move beyond three Higgs doublets, at least as long as one keeps the renormalizability. Thus, even though there seems to be no obstacles a priori, one should demonstrate explicitly how such models based on yet higher-order CP -symmetries can be built, and what novel features they involve. In this paper we perform this task. We build two examples of five-Higgs-doublet models based on CPsymmetries CP8 and CP16. Assuming that these symmetries are respected by the minimum, we derive scalar mass spectrum and discuss the properties of the DM candidates. The five doublets are grouped in a natural way: one Higgs doublet acquires the vacuum expectation value (vev) and produces the SM-like Higgs particle, while the inert sector includes two pairs of two doublets, with the CP -transformation mixing the doublets within each pair. When constructing these examples, we will explain the strategy of building models with even higher-order CP -symmetries, should an interest in such models appear.

II. NHDMS WITH HIGHER ORDER CP
A. The freedom of defining CP -symmetries A self-consistent local quantum field theory does not uniquely specify how discrete symmetries, such as C and P , act on field operators [1,11,28,29]. There is freedom in defining these transformations, which becomes especially large in the case of several fields with equal quantum numbers. These fields are not physical by themselves; any linear combination of those fields which preserves the kinetic terms will be equally acceptable as a basis choice for the theory. Therefore, any symmetry of the Lagrangian which is supposed to incorporate a physically measurable property, is defined up to an unconstrained basis change.
Focusing on several scalar fields φ i , i = 1, . . . , N with equal quantum numbers, one can define the CPtransformation as in (1). If there exists a unitary matrix X such that the Lagrangian and the vacuum of a model are invariant under this transformation, then the model is CP -conserving in the very traditional sense that all CP -odd observables are zero, and the transformation (1) plays the role of "the CP -symmetry" of the model [1]. It is only when none of transformations (1) is a symmetry of the model that we say that CP -violation takes place.
Using the basis change freedom, it is possible to bring the matrix X to a block-diagonal form [9,11], which has on its diagonal either unit entries or 2 × 2 matrices of the following type: as in Ref. [9], or 0 e iα e −iα 0 as in Ref. [11].
This is the simplest form of X one can achieve with basis transformations in the scalar space C N . Applying the transformation (1) twice, one obtains a pure family transformation a = XX * . If X contains at least one 2 × 2 block with α = 0 or π, then a = δ ij , which means that the CP -transformation (1) is not an order-2 transformation. If k is the smallest integer such that a k = δ ij , then we get the CP -transformation of order 2k, which we denote CP2k, and the resulting family symmetry group Z k , which is generated by a, the square of the CP -transformation. As we explained in the introduction, in order to avoid accidental symmetries, one needs to consider only 2k = 2 p .

B. The strategy
Before moving to specific examples, let us first outline the strategy of building N -Higgs-doublet models whose only symmetries in the scalar sector are CP -symmetries of orders 2k and their powers.
One starts by writing the Higgs potential as a sum of rephasing-invariant and rephasing-sensitive parts, V = V 0 + V 1 . The rephasing-invariant part can be generically written as with all the coefficients being real. The rephasing-sensitive part V 1 contains only those quadratic and quartic combinations which are invariant under the rephasing transformation a.
Although for small values of k and N , the phase-sensitive part of the potential can be quickly constructed by trial-and-error, there exists an algorithmic procedure described in [20] which allows one to build V 1 for a chosen rephasing symmetry group Z k with a given number of doublets N . Of course, not all discrete groups can be implemented. In the same work [20], it is proven that, staying with N doublets and renormalizable potentials, one can implement cyclic groups Z k of order k ≤ 2 N −1 . Trying to impose any symmetry whose order is larger than this bound unavoidably leads to accidental continuous symmetries. Thus, the order of generalized CP -symmetry in NHDM cannot exceed 2k = 2 N .
Even if the potential V 1 is constructed with guess, one can always find its full rephasing symmetry group via the systematic procedure based on Smith normal forms, which was developed in [20] and explained in less technical fashion in [30]. This computation can be done by hand or implemented in a computer-algebra code. It is in this way that one verifies the absence of accidental rephasing symmetries. The absence of other symmetries beyond rephasing ones is guaranteed by the fact that all free parameters in (4) are independent.
Next, having the potential invariant under Z k generated by a, one needs to check what additional conditions on its parameters one must impose to make it invariant under the desired CP -symmetry of order 2k. Since higher-order CP -transformations mix pairs of doublets, there arise obvious conditions on the parameters of V 0 such as m 2 22 = m 2 33 , etc. In addition, the parameters of V 1 are also constrained. These constraints can be analyzed term by term.
However, instead of such analysis, we will proceed in a more efficient way. We will first construct all bilinear combinations φ † i φ j and classify them according to their CP -charge q defined modulo 2k. Within each sector with definite q, there may exist several bilinears r a , all of them transforming in the same way under CP : It is sufficient to list only bilinears with 0 ≤ q ≤ k; the complex conjugated bilinears r † a with CP -charges −q will fill all other charge assignments from k to 2k. In terms of these bilinears, the total potential can be schematically written as where the non-zero coefficients M a span only those r a with q = 0, and the hermitian matrix Λ ab is blockdiagonal, with unconstrained blocks within each q sector. Once again, it is important to check that the resulting CP2k-invariant potential does not acquire any accidental symmetries. The rephasing symmetry group can again be unambiguously found with the Smith normal form technique [20,30]. The absence of the usual CP -symmetry is guaranteed by the fact that the hermitian matrix Λ ab in (6) cannot be made real by any basis change. Absence of other accidental symmetries beyond rephasing is assured by the fact that the matrix Λ ab has sufficiently many independent free parameters.
Since in this work we do not aim at producing minimal models but rather look for examples of CP -protected scalar dark matter candidates, we will make sure that it is possible to conserve this symmetry upon minimization of the Higgs potential. This will lead us, both for CP8 and CP16, to models with five Higgs doublets: one SM-like φ 1 and four inert ones φ i , i = 2, 3, 4, 5. These inert doublets form two pairs, (φ 2 , φ 3 ) and (φ 4 , φ 5 ), which get mixed by the CP8 or CP16-transformation. In each case, we will take the CP -conserving vev alignment v 1 = v, v 2,3,4,5 = 0, expand the potential around the minimum, and calculate the neutral and charged scalar mass matrices. We will confirm the general observation that the physical scalar fields in the inert sector are pairwise mass-degenerate, just as in CP4 3HDM. For neutral scalars, we will combine pairs of real mass eigenstates into complex neutral fields with definite CP -charge q and briefly discuss the emerging self-interaction pattern.

C. 3HDM is not enough
It is instructive to begin the study by demonstrating why 3HDM fails to accommodate the CP8-symmetry [20]. According to the general strategy, one first needs to write a model with rephasing symmetry Z 4 and then extend the symmetry to CP8. For three Higgs doublets φ i , i = 1, 2, 3, the Z 4 group of symmetries is generated by the transformation a 4 which, after an appropriate basis change, can be represented as Here, dots stand for the zero entries. The Higgs potential is written as V = V 0 + V 1 , where V 0 given in (4) and while the phase-sensitive part V 1 Here, both coefficients can be complex and must be non-zero. If at least one of them is zero, then the number of independent phase-sensitive terms drops below N − 1, and the potential acquires a continuous rephasing symmetry [20]. We now want to require that this potential be invariant under CP8, which is generated by φ i CP − − → X ij φ * j of order 8. The matrix X can be brought by a basis change to the form One immediately checks that applying CP8 twice produces XX * = a 4 from Eq. (7). Since CP8 mixes the doublets φ 2 and φ 3 , one must equate their respective coefficients in V 0 . In addition, one requires that V 1 stays invariant under CP8. Straightforward algebra shows that under CP8 Therefore, one must set λ = 0 to assure CP8-invariance of V 1 . Since we are left with only one rephasing-sensitive term, the potential acquires a continuous U (1) rephasing symmetry. Therefore, the true symmetry content of the resulting model is not the discrete group generated by CP8 but the continuous group of arbitrary phase rotations and the usual CP-transformation. Colloquially speaking, 3HDM potential does not offer enough room to incorporate CP8 without producing accidental symmetries.
In the Appendix we show that this observation generalizes to NHDM with any N . If one takes the largest cyclic group possible for NHDM, Z k with k = 2 N −1 , and calculates for all Higgs doublets φ i their q i charges associated with the rephasing group Z k , then one finds a very characteristic pattern of these charges, which involves successive powers of 2. However, if one starts with a CP -symmetry of order 2k, then one arrives at the same symmetry group Z k with a very distinct pattern of charges: for any doublet with charge q i there exists a doublet with charge −q i . These two patterns do not match. It means that trying to impose CP -symmetry of order 2k = 2 N on NHDM leads to a continuous symmetry, which ruins the construction.

III. BUILDING 5HDMS WITH CP8
A. 5HDM with CP8 The five-Higgs-doublet model 5HDM can incorporate cyclic groups with order up to 16. By the arguments exposed in the Appendix, the maximal cyclic symmetry Z 16 cannot be extended to CP32. However, 5HDMs with CP8 and CP16 are well possible, and in this and the next sections, we construct such models.
The 5HDM uses N = 5 Higgs doublets φ i , all with the same gauge quantum numbers. Similarly to the previously considered 3HDM case, we define, in the appropriate basis, the generator a 4 of the group Z 4 and the matrix X which defines CP8: with the same η ≡ e iπ/4 . The relation a 4 = XX * still holds. The scalar potential can again be written as a sum of phase-invariant and phase-sensitive parts V = V 0 + V 1 , with V 0 as in (4), where one implicitly assumes that the coefficients m 2 ii , λ ij , and λ ′ ij respect the symmetry under the simultaneous exchange φ 2 ↔ φ 3 and φ 4 ↔ φ 5 . The phase-sensitive part is now much richer than in 3HDM due to the fact that we have two 2 × 2 blocks in the definition of X.
Following the strategy outlined in the previous section, we write down all N 2 = 25 gauge-invariant bilinears φ † i φ j and build out of them combinations r a which are CP8-eigenstates, that is, which transform under CP8 as in (5). Using the shorthand notation i ≡ φ i , we list these CP8-eigenstates according to the value of q: CP8-even, q = 0 : 1 † 1, 2 † 2 + 3 † 3, 4 † 4 + 5 † 5, 2 † 4 + 5 † 3, 4 † 2 + 3 † 5, q = 1 : Bilinears with CP -charges from 4 to 8 are obtained by complex conjugating the states listed here. Notice that the combinations with q = 0 and q = k = 4 fall in the traditional classification of CP -even/odd states. Thus, they can be coupled with other CP -even or odd operators of the model in a CP -conserving way. In terms of these bilinears, the total potential is schematically written as in (6). Using the methods outlined above, one can verify that this potential indeed does not possess any other symmetry.

B. Charged Higgs masses
When minimizing the potential, we focus on the case of CP8-conserving vacuum, which implies that only the first doublet acquires a vev v 1 = v. All inert Higgs doublets are expanded as The terms of the Higgs potential which generate the scalar masses are The SM-like Higgs boson acquires mass m 2 h = −2m 2 11 = 2λ 1 v 2 . In the inert sector, we begin with the charged Higgs masses, for which only the first two lines are relevant, and obtain the following mass terms: The charged mass matrix splits into two blocks 2 × 2 within subspaces (H ± 2 , H ± 4 ) and (H ± 3 , H ± 5 ), with exactly the same eigenvalues in each block. Thus, the charged Higgs spectrum becomes pairwise mass-degenerate.
Instead of explicitly diagonalizing each block, one can take one step back and simplify the starting potential without loss of generality. Indeed, the pairs of doublets (φ 2 , φ 3 ) and (φ 4 , φ 5 ) transform in exactly the same way. Therefore, one can perform basis transformations that mix φ 2 and φ 4 by unitary matrix U and, simultaneously, φ 3 and φ 5 by unitary matrix U * , and this basis change keeps the symmetry transformations a 4 and X unchanged. This freedom of basis change is always there, and it allows us to find such U which removes the cross term φ † 2 φ 4 + φ † 5 φ 3 altogether. In that basis, we still have the same potential as before, but with reparameterized coefficients. In particular, the charged Higgs masses will now be given only by the first line of (15). The four charged Higgses then have the following masses:

C. Neutral Higgs masses and CP8-eigenstates
For neutral Higgses, instead of explicitly expanding all the doublets into real components, it is convenient to define neutral complex fields which are already CP8-eigenstates, in similarity to the states ϕ and Φ in CP4 3HDM. These fields can be read off the table (12); they correspond to the bilinears with q = 1 and q = 3 in which φ 0 1 set to its vev φ 0 1 = v/ √ 2: The two sectors corresponding to q = 1 and q = 3 do not mix in the mass matrix. Staying in the charged Higgs eigenstate basis defined above, we can represent the mass terms as where the two mass matrices are By diagonalizing them, we get the four different values for the neutral Higgs masses for the fields ϕ, Φ in the q = 1 sector, with m ϕ < m Φ , and ψ, Ψ in the q = 3 sector, with m ψ < m Ψ . If needed, these fields can also be written in terms of real components. In that case we have eight real fields which are pairwise mass degenerate. The scalars from the q = 1 and q = 3 sectors can interact with the Z-boson via Zϕ i ψ j vertices. Indeed, each inert Higgs doublet φ i produces, via its kinetic term, the term (ḡ/2)Z µ (h i ∂ ↔ a i ), whereḡ = g 2 + g ′2 is the combined gauge coupling and h ∂ ↔ a = h(∂ µ a) − a(∂ µ h). When expressed in terms of CP8-eigenstates (17), these interaction terms become These vertices represent the CP8-counterpart of the vertices ZHA in the 2HDM and ZϕΦ in the CP4 3HDM. They conserve the CP8 quantum number: the sum of the internal CP8-charges of ϕ i and ψ i is 1 + 3 = 4, which, for CP8-symmetry, is equivalent to being CP-odd. After diagonalization of the mass matrices in the ϕ i and ψ i sectors, these interactions render the next-to-lightest state metastable: To avoid confusion, we stress that notation ϕ * denotes the state conjugated to ϕ (which, contrary to the usual expectation, is not the antiparticle to ϕ, see detailed discussion in [27]), while Z ( * ) denotes a real or virtual Z-boson. The only exception is when the lightest states in these two sectors are orthogonal, which would forbid Zϕψ-vertex and render both scalars stable. The self-interaction in the inert sector leads to several interaction terms involving ϕ and ψ: All coefficients here are independent; λ 13 , λ 31 , and λ 22 can be complex, but it does not imply CP -violation, because the scalars here are themselves CP -eigenstates, and the CP8-charge is conserved by each term separately. These interactions switch on new channels for two identical DM candidates: although the direct annihilation ψψ → SM is forbidden by the CP8 conservation, the semi-annihilation processes ψψ → ψ * ϕ is allowed for m ψ > m ϕ . Finally, for sufficiently large mass splitting, the direct triple decays are also allowed:

IV. BUILDING 5HDM WITH CP16
In this section, we build yet another version of 5HDM, the one with CP16. We use the same strategy as before, but with the new parameter ξ ≡ exp(iπ/8) instead of η = exp(iπ/4). The first attempt is to use the same structure for X as before: and classify the bilinears according to their CP-charges q defined modulo 2k = 16. This classification is similar to Eq. (12) but with an important change: The key difference with respect to CP8 case of Eq. (12) is that now, with k = 8, the charges q = k − 2 and q = 2 are distinct, and the matrix Λ ab does not mix them. Thus, Λ ab stays block diagonal, with blocks corresponding to subspaces of distinct values of q. But then the structure of Λ ab does not depend on the value of k provided k > 4. It means that the same potential is invariant not only under CP16 but also under any higher-order CP2k, as well as under the continuous U (1) transformations generated by a 8 in Eq. (24) with η replaced by any phase rotation. In short, the CP16 leads to accidental symmetries including U (1) and the usual CP. However, we can try another quantum number assignment: Within CP8, this assignment could be reduces to the previously considered one by rephasing within the last block, while here it leads to an essentially different model. We also remark that it is absolutely inessential which block is the third power of which. One can equally well denote η ′ ≡ η 3 and then observe that η = η 9 = (η ′ ) 3 . Group-theoretically, this reflects the fact that the four inert doublets are transformed under a 8 by the four distinct generators of the rephasing group Z 8 : η, η 3 , η 5 , and η 7 . Again, classifying the bilinear transformations for CP2k, we get: CP-even, q = 0 : 1 † 1, 2 † 2 + 3 † 3, 4 † 4 + 5 † 5 q = 1 : In the case of CP16, the value of k = 8. Then, the CP-charges q = k − 4 and q = 4 are identical, and so are the charges q = k − 2 and q = 6. The corresponding bilinears can be coupled via Λ ab , and the resulting potential will not have the continuous symmetry.
Using the Smith normal form technique, one also verifies that the rephasing symmetry of this model is indeed Z 8 . In order to check that the model does not accidentally acquire the usual CP -symmetry, let us notice that the hermitian matrix Λ ab has three complex off-diagonal entries, coming from the 2 × 2 blocks with charges q = 2, 4, 6. They generate six different rephasing-sensitive terms in the potential. Using the rephasing freedom, one can make two of these entries real, but not all three of them. Thus, the coupling matrix Λ ab cannot be made real in any basis. This fact forbids the usual CP -symmetry as well as the symmetry under φ 2 ↔ φ 3 , φ 4 ↔ φ 5 . Thus, we have a viable 5HDM with CP16 and no accidental symmetry.
We proceed to the mass spectrum calculation for the case of unbroken CP16. The terms in the Higgs potential that generate the scalar masses are very similar to Eq. (14) but with a few terms omitted: Therefore, we can reuse exactly the same formulas for scalar masses as before, Eqs. (16) and (19), but just set λ 7 = λ ′ 7 = 0 in the latter. The matrices in (19) become diagonal, which is to be expected because the four complex neutral fields carry now all distinct CP-charges: The interaction vertices Zϕ i ψ j remain as in Eq. (20), and they are already written in terms of mass states. All these vertices still conserve the CP16-charge q. These vertices lead to decays (21) within the 23 and 45 subsectors, but they do not couple the two sectors, which renders the lightest states in the two inert subsectors stable.
In any of these cases, there exists an interaction term involving asymmetric combinations of the two fields, which can lead to semi-annihilation and decays, in similarity to what we found in the CP8 case.

V. DISCUSSION AND CONCLUSIONS
The key message of this study is that it is well possible to construct renormalizable multi-Higgs models whose symmetry content is given only by a higher-order CP and its powers. The CP4-symmetric 3HDM proposed initially in [8] is the simplest example of this kind, but it is not the only possibility. We have constructed here two versions of 5HDM with CP8 and CP16, and the methods we have used can be generalized to CPsymmetries of any order 2k = 2 p , should the need arise.
If the vacuum respects the higher-order CP -symmetry, then the real scalars emerging from the inert sector are pairwise mass-degenerate and can be grouped into complex neutral fields ϕ i which are CP -eigenstates. Just as in the CP4 3HDM example, their CP -properties are described by CP2k-charges q i defined modulo 2k. They generalize the notion of CP -parity (that is, CP -charge defined modulo 2) for the usual CP -symmetry of order 2.
The lightest scalar in the inert sector serves as the DM candidate, and its stability is insured by the exotic CP -symmetry rather than internal symmetry. Models with elaborate CP sectors, such as CP16 5HDM, can contain two or more DM candidates with different CP -charges.
One may ask whether there is any difference between models based on CP2k, considered here, and the more traditional multi-Higgs models based on rephasing symmetries of order 2k, see examples in [31][32][33]. Despite the symmetry group in both cases is the same, Z 2k , there are several distinctions.
First, the CP2k-based models possess vertices of the type Z µ ϕ i ∂ ↔ ϕ j , where q i + q j = k = 0. This is possible because the Lorentz structure of this interaction term is by itself CP -odd and requires the internal CP -properties of the two fields ϕ i ϕ j to organize themselves into a CP -odd combination. In the usual case, for example in the CP -conserving 2HDM, the corresponding vertex is ZHA, where H is CP -even and A is CP -odd, so that their product is CP -odd. In the CP2k-symmetric models, one just arranges q i + q j = k, which exactly corresponds to being CP -odd. In the traditional Z 2k -symmetric model, such vertices are impossible because the symmetry is internal, and therefore the Z 2k -charges of fields in any vertex must be a multiple of 2k.
Second, the mere fact that a complex field ϕ is a CP -eigenstate means that it is its own antiparticle. The one-particle state ϕ * |0 is not an antiparticle to ϕ|0 but is rather a different particle with the same mass. This doubling of spectrum is only possible for zero-charge fields and is discussed at length in [27].
This feature allows one to consider an asymmetric DM evolution regime, in which ϕ dominates over ϕ * . However, unlike the traditional asymmetric DM models [34,35], this imbalance does not imply particleantiparticle asymmetry. Constructing a model which exhibits such an imbalance and studying its late-time observational signatures is a task for future investigation.
Having demonstrated that it is possible to build CP -conserving models based on various CP2k-symmetries, one can ask whether this distinction is observable in any imaginable experiment. If it is, we arrive at the exciting possibility of determining experimentally the order of the CP -symmetry which the real world is closest to. This question was already posed in [8], but it remains unanswered.
On the theoretical side, it is interesting to see if a CP2k-symmetry can arise as a low-energy residual symmetry from a more symmetric model at high energy scale, whose high symmetry spontaneously breaks down at lower energies. All existing models of this kind generate at lower energies only usual family symmetries, not an exotic CP . On the other hand, as established in [18,19], certain symmetry groups G not only allow but even require the CP -symmetry to be of higher-order. Thus, equipping a symmetric model with higher-order symmetries is not a problem; one just needs to make sure it is the only symmetry to survive at low energies. If such a construction leading to a residual CP2k-symmetry is found, it may provide a natural explanation how the CP8 or CP16 5HDMs could emerge from a highly symmetric construction with one scalar singlet φ 1 and one quadruplet (φ 2 , φ 3 , φ 4 , φ 5 ). It will then serve as an additional motivation to look deeper at this exotic form of CP -conservation and its observable consequences. modulo k. Different q i spectra will produce nonequivalent models with the same Z k symmetry (in fact, we already encountered this situation in section IV when we were building CP16 5HDM). However, we prove below that there exists, up to permutation, a unique assignment of charges, for which the potential does not acquire accidental continuous symmetries.
Next, we briefly recap the technique based on the Smith normal form (SNF), which was developed in [20] to establish the rephasing symmetry group of any potential and the exact form of the NHDM potential with the maximal cyclic symmetry group Z k , k = 2 N −1 .
For any scalar potential one first checks how each individual term changes under a generic global rephasing transformation φ j → e iαj φ j . The i-th term picks up the phase factor d ij α j , where the integer coefficients d ij are immediately read off the expression for the i-term. For example, if the first term is (φ † 2 φ 1 )(φ † 3 φ 1 ), its coefficients are d 1j = (2, −1, −1, 0, . . . , 0). Going through all m rephasing-sensitive terms, one builds in this way the coefficient matrix d ij , which is an integer-valued rectangular matrix m × N . Then one can apply a sequence of certain elementary steps and reach its Smith normal form (SNF). The Smith normal form exists and is unique for any rectangular matrix with integer coefficients. The diagonal entries of the SNF immediately give the rephasing symmetry group of the potential.
The explicit form of the Z k -symmetric NHDM was given in [20]: Notice that each row has the following properties: j d ij = 0, which reflects the global overall rephasing symmetry, a subgroup of U (1) Y , and j |d ij | = 4, which reflects the fact that all interaction terms used here are quartic. The SNF of this matrix has on its diagonal the sequence (1, 1, . . . , 1, 2 N −1 ), which indicates the rephasing symmetry group Z 2 N −1 , in addition to the overall U (1) Y rephasings of all doublets. Since the SNF is unique and since its construction is invertible, any NHDM potential with the same symmetry group Z 2 N −1 without any accidental symmetry can be brought to this form by an appropriate basis transformation.
The Higgs potential encoded in this matrix is where all coefficients can be complex. From this expressions one can immediately obtain the spectrum of Z k charges: q i = (1, 2, 4, · · · , 2 N −2 , 0).
All charges are defined modulo k = 2 N −1 . There is a freedom in simultaneous shift of all charges by the same value, but their differences remain as they are. Up to permutation, this is the only charge assignment which is compatible with Z k and is capable of generating N − 1 different terms, thus avoiding accidental continuous symmetries. The potential (A2) has N − 1 terms constructed of N Higgs doublets. There exist no other renormalizable terms invariant under the same symmetry. This can be seen from the matrix d itself. Suppose there exists yet another term, which would appear in this matrix as N -th row d N j . Since the matrix is of rank N − 1, the new row can be written as a linear combination of the existing rows with integer coefficients. It is immediately seen by direct inspection that any such combination would produce a row with j |d N j | > 4. Therefore, any such term can only be of higher order.
Suppose now we want to impose the CP -symmetry of order 2k. Denoting the generator of the cyclic group Z k by a and the generator of the CP2k-symmetry by J, we are looking for such a construction which satisfies J 2 = a. As mentioned in the introduction, a higher-order CP -transformation J acts on doublets as φ i → X ij φ * j , where the matrix X can be brought to the block-diagonal basis, with the diagonal containing either entries 1 or 2 × 2 blocks of the form [11] 0 e iα e −iα 0 .