CP4 miracle: shaping Yukawa sector with CP symmetry of order four

We explore the phenomenology of a unique three-Higgs-doublet model based on the single CP symmetry of order 4 (CP4) without any accidental symmetries. The CP4 symmetry is imposed on the scalar potential and Yukawa interactions, strongly shaping both sectors of the model and leading to a very characteristic phenomenology. The scalar sector is analyzed in detail, and in the Yukawa sector we list all possible CP4-symmetric structures which do not run into immediate conflict with experiment, namely, do not lead to massless or mass-degenerate quarks nor to insufficient mixing or CP-violation in the CKM matrix. We show that the parameter space of the model, although very constrained by CP4, is large enough to comply with the electroweak precision data and the LHC results for the 125 GeV Higgs boson phenomenology, as well as to perfectly reproduce all fermion masses, mixing, and CP violation. Despite the presence of flavor changing neutral currents mediated by heavy Higgs scalars, we find through a parameter space scan many points which accurately reproduce the kaon CP-violating parameter $\epsilon_K$ as well as oscillation parameters in K and $B_{(s)}$ mesons. Thus, CP4 offers a novel minimalistic framework for building models with very few assumptions, sufficient predictive power, and rich phenomenology yet to be explored.


Introduction
The discovery of the Higgs boson [1,2] was a major milestone in particle physics: the last remaining piece of the Standard Model (SM) was brought to light. From this point onwards, any deviations from the observed SM-like character of particle physics could well be a sign of new physics hitherto undiscovered. And although further LHC measurements have shown that the discovered 125 GeV scalar is conforming to the SM behaviour (see Ref. [3] for the combined LHC Run I results; and the presentations of the latest Run II results at the recent EPS-HEP 2017 conference in Venice), the SM leaves too many phenomena unexplained to be regarded as a satisfactory theory. These include: the absence of dark matter (DM) candidates; complete ignorance of the origin of neutrino masses [4] and explanation of their observed mixings; inability to explain the origin of the observed CP violation (CPV) [5]; and incapacity to cast any light on the quark and lepton mass and mixing hierarchies.
Addressing these unsolved questions requires going beyond the Standard Model (BSM), and many proposed solutions involve expanding the minimal scalar sector of the SM by adding extra scalar fields; for a recent overview of models with extended Higgs sectors, see Ref. [6]. A direction, which is particularly attractive due to its conceptual simplicity, is to stay with Higgs SU (2) doublets and to extend the notion of generations to the scalar sector. One arrives in this way to N -Higgs-doublet models (NHDM). The Two-Higgs Doublet Model (2HDM) (see Ref. [7] for a review) is the most popular example but models with more Higgs doublets are also being actively explored. Such extended scalar sectors have a rich phenomenology but they come with a price: as the number of extra fields grow, so does the number of free parameters, and thus the predictive power of the theory is reduced. When building such models, it would be desirable to achieve a balance between the following two requirements: making as few extra assumptions as possible, on top of those in the SM, and obtaining a model which satisfies all experimental constraints, while also being able to make testable predictions for the ongoing or future measurements. One would like to avoid describing all available data at the cost of an excessive number of new fields and parameters, but also to avoid obtaining a neat BSM model, so tightly constrained by theoretical constraints that it fails a comparison with the experimental bounds.
A successful way of reducing the number of free parameters of BSM models, thus increasing their predictivity and even resolving some of their theoretical problems, is by imposing additional global symmetries, either continuous or discrete [4,8]. For instance, the most general 2HDM has Higgs-mediated flavour-changing neutral currents (FCNC), but by enforcing its Lagrangian to be invariant under a discrete Z 2 symmetry, those undesirable interactions are made to vanish [9][10][11]. Generically, a typical N -Higgs-Doublet Model (NHDM) will contain hundreds of free parameters in its scalar and Yukawa sectors. By imposing large non-Abelian discrete symmetry groups this number of parameters may be reduced to about a dozen, making such a model rather predictive. However, even though it is reasonably easy to fashion models with an acceptable scalar sector -i.e., models which include a scalar state of mass 125 GeV with a SM-like behaviour and extra scalars that are not yet excluded by the LHC searches -such models usually render fermionic sectors which are unphysical [12]. Indeed, for sufficiently large discrete symmetry groups, there is always some residual symmetry preserved by the vacuum, which will imply either massless or mass-degenerate fermions, or alternatively, lead to an inadequate quark mixing, or an insufficient CPV. Smaller symmetry groups may lead to good experimental fits, but they usually leave the model with an excess of free parameters, making such theories cumbersome to analyse and less attractive as an alternative to the SM.
Recently, in Refs. [13,14] a new type of a multi-Higgs model was proposed -based on a single symmetry requirement which, rather surprisingly, leads to well-shaped scalar and fermion sectors. The symmetry assumption behind this model is very simple: The minimal multi-Higgs-doublet model implementing a CP -symmetry of higher order without producing any accidental symmetry.
Two aspects here require clarification. First, the order k of a symmetry transformation is the minimal number of times one needs to apply it to arrive at the identity transformation. The traditionally defined CP -transformation is of order 2, but higher-order (generalized) CP -symmetries can also be defined and used for model building, see Section 2.1. Second, the concept of accidental symmetries, which is familiar to model builders in BSM physics, refers to the situation when requiring that a given model is invariant under a certain type of symmetry produces a Lagrangian which is, in fact, invariant under larger symmetry groups. Appearance of accidental symmetries signals a certain lack of control over the symmetry content of a model or its phenomenological manifestations.
The above theoretical request leads to a Three-Higgs Doublet Model (3HDM) with an order-4 CP -symmetry, which we label CP4. This model has a scalar potential which was initially written in Ref. [15] and explored, in the case of unbroken CP4, in Ref. [13]. If CP4 stays intact, the model incorporates a novel feature in BSM models: a complex scalar field which, although being a CP -eigenstate, is neither CP -even nor CP -odd but rather, in this case, a CP -half-odd state. In Ref. [14], it was shown that the CP4 symmetry can be extended to the Yukawa sector in a satisfactory manner. However, in order to avoid unobserved mass-degenerate fermions, the CP4 symmetry must be spontaneously broken. The question of whether the CP4 3HDM Lagrangian can fit the current experimental data was left unanswered up until now, and that is the main objective of this work.
In short, the CP4 3HDM emerges as an interesting candidate of a model which keeps a fair balance between the minimality of its theoretical assumptions, and its phenomenological richness and predictivity. We firmly believe that a deeper study of this model will reveal hidden features which are not present in other models.
In this paper, we will undertake the first detailed phenomenological exploration of the CP4 3HDM with the spontaneously broken CP4 symmetry. Our main objectives will be: • To prove that the model has a scalar sector which complies with known experimental results, at least for the SM-like scalar state of mass 125 GeV.
• To show that the Yukawa sector of the model, after spontaneous breaking of CP4, can fit the fermion masses, mixing, and observed CP -violation quantities.
• To analyse the inevitable tree-level FCNCs of the model and to show that the model has regions of the parameter space for which they are under control and conform to the experimental measurements.
In Section 2, we will discuss the scalar sector of the CP4 3HDM, including the extremization conditions, the scalar mass matrices, and the conditions for an alignment limit in the scalar sector, which will produce a SM-like scalar state keeping its FCNCs under control. In Section 3, we will construct all CP4-symmetric Yukawa sectors which do not run into an immediate conflict with experiment. We then present in Section 4 the results of a dedicated and efficient scan of the scalar and Yukawa parameter space of the model, showing that there exist parameter space regions for which the most stringent experimental constraints on FCNCs are satisfied. We close this study with a discussion and conclusions in Section 5. Several Appendices provide supplementary details.

The scalar sector of CP4 3HDM
We will be working with three SU (2) doublets (i.e. 3HDM), each with hypercharge Y = 1 and denoted as φ 1 , φ 2 and φ 3 . A priori, the scalar potential of such a model may have (before one uses the liberty to redefine the fields of the model) 54 independent real parameters, as opposed to the two parameters of the SM scalar potential or 14 free parameters of the general 2HDM scalar sector. The imposition of global-discrete or continuoussymmetries on this model is therefore a good idea. The proposal by Weinberg of a 3HDM equipped with natural flavor conservation in the Yukawa sector [16], for instance, included two discrete symmetries generated by φ 2 → −φ 2 and, separately, φ 3 → −φ 3 , making the symmetry group of the model Z 2 × Z 2 . This symmetry is then spontaneously broken by the vacuum expectation values (vevs) in the Higgs doublets. That model had a total of 18 real parameters. Our approach is based not on family symmetries like the Z 2 ones described in the above example, but rather on generalized CP symmetries (or GCP), which relate the fields with their complex conjugates. For the reader's convenience and also to set up the notation, we begin with a general reminder on the freedom in choosing CP transformations one has when building a model.

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 [5,[17][18][19]. There is a great amount of freedom in defining these transformations, which becomes especially large in the case of several fields with equal quantum numbers. This is due to the fact that such fields are not physical by themselves; only the mass eigenstates obtained after spontaneous symmetry breaking will correspond to physical particles. Any linear combination of those fields which preserves the kinetic terms of the model 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 choice shift. Focusing now on a CP transformation acting on several scalar fields φ i , i = 1, . . . , N , one often considers the following GCP transformations [20,21]: If there exists a unitary matrix X such that the Lagrangian of a model is invariant under this GCP transformation, then the model is explicitly CP -conserving and J X can play the role of "the CP -symmetry" of the model [5]. Notice that the "conventional" definition of , is only one of many possible definitions and is, in fact, a basis-dependent choice.
We find this terminological issue so important that, in abuse of the reader's patience, we spell it out once again. When we say that the model is CP -conserving, we may refer to any form of GCP symmetry (2.1), with whatever fancy X. In particular, of the "conventional" definition of CP transformations fails to leave the lagrangian invariant, but a more complicated GCP transformation does, then the model is still CP -conserving in the very traditional sense that all CP -odd observables are zero. It is only when none of the transformations (2.1) is a symmetry of the model that we say that CP violation takes place [5].
The shape of the X matrix may have important consequences for the phenomenological behaviour of the models. In the 2HDM, for instance, three different (and relevant) choices for X are possible, each leading to different accidental symmetries, and three different CPinvariant models with vastly different phenomenologies emerge from those choices [22,23].
Notice now that applying J X twice generates a pure family transformation: Using the redefinition freedom one has in the choice of the basis of scalar fields, it is possible to bring the matrix X to a block-diagonal form [19,20], with the blocks being either 1 × 1 phases or 2 × 2 matrices of the following type: [20], or 0 e iα e −iα 0 as in Ref. [19]. (2.3) This is the simplest form of X one can achieve with basis transformations in the scalar space C N . If X contains at least one 2 × 2 block with α = 0 or π, then (J X ) 2 = XX * = I. This then means that the CP transformation (2.1) is not an order-2 transformation. If k is the smallest integer such (J X ) k = I, the GCP transformation J X is said to be of order k.
One immediately sees that k is necessarily an even number: one needs to perform conjugation an even number of times to obtain the identity transformation. However, imposing the GCP of a generic even order k immediately leads to accidental symmetries, including the GCP of a smaller order. Indeed, if k has prime factors other than two, one can factor them out and obtain a smaller-order GCP. The only way to prevent this possibility is to take k = 2 p , with p ≥ 1 an integer number. Which means that the usual CP is of order two (CP2), the first non-trivial higher-order CP symmetry is CP4, the next one is CP8, and so on.
Here we would like to reiterate again the point made three paragraphs earlier. When we label a model as CP -conserving or CP -violating, we do not need to specify whether it conserves or violates CP2, CP4, or a higher-order GCP. CP -odd observables do not distinguish them. If there exists at least one GCP transformation that leaves the model invariant, then it is CP -conserving. Conversely, when we speak of CP -violation, we mean that all possible GCP transformations fail to leave the model invariant.
Although the CP -odd observables do not distinguish different classes of GCP transformations (CP2, CP4, etc), the parameters of the lagrangian definitely do. Since higher-order GCPs involve transformation between a pair of fields, imposing it will certainly constrain the parameters stronger than the conventional CP or, in general, any expression for CP2. As a result, it may happen that the model itself does not offer enough freedom to implement a higher-order CP symmetry. For example, imposing a higher-order CP symmetry on the scalar potential of the 2HDM produces accidental symmetries, which include the usual CP [23]. Thus, imposing higher-order CP symmetry has always been viewed as a compact way of defining a model, but not as a path towards new models that could not be achieved through the usual "order-2 CP + family symmetry" combination. A rare exception is discussed in Ref. [24], where the higher-order CP symmetries were classified as distinct opportunities for model building. Further, extending these GCP symmetries to the Yukawa sector within 2HDM was problematic [25][26][27], as they ran into trouble when confronted with the experimental data (predicting some massless fermions and an insufficient CPV, for instance). In a sense, 2HDM does not offer the model builder enough room to fully incorporate such a strongly constraining symmetry as CP4, and one needs to extend the number of doublets to at least three.

The scalar potential
How many different global symmetries can one impose upon a given BSM model? Given the basis redefinition freedom present in many such models, apparently different symmetries may, in fact, be related by basis choices. For instance, in the 2HDM a symmetry of the form φ 1 ↔ φ 2 is equivalent, in a different basis, to the usual Z 2 symmetry, φ 1 → φ 1 and φ 2 → −φ 2 . In each basis, however, the Lagrangian of the model looks completely different, with seemingly diverse relations between parameters. In the 2HDM, the work of Ref. [22] proved that there are only six different symmetry classes, which are not related among themselves by basis choices.
In the 3HDM the situation is more complicated. In Ref. [28] a first attempt at finding different classes of the 3HDM symmetries was undertaken, but that study has been restricted to simple Abelian groups. A systematic and constructive search for all discrete symmetry groups in the scalar sector of the 3HDM was performed in Ref. [15] for Abelian and in Ref. [29,30] for discrete non-Abelian groups. In each case, it was checked whether the family symmetry group can be further enlarged to include a general CP symmetry without producing any further accidental group.
This construction showed that, up to a basis choice, there exists only one 3HDM with a CP symmetry of higher order, to be specific, CP4, which does not lead to accidental symmetries. In the suitable basis, CP4 acts on the three Higgs doublets in the following way: The most general potential respecting this CP4 symmetry can be written as with all parameters being necessarily real, and with real λ 5 and complex λ 8 , λ 9 1 .
Applying the transformation (2.4) twice leads to J 2 = XX * = diag(1, −1, −1) = I. It is trivial to see that one recovers the identity transformation only after applying J four times: J 4 = I. Thus, the transformation J is indeed a GCP of order 4. For generic values of the coefficients, this potential has no other Higgs-family or GCP symmetries, apart from powers of J [15]. In particular, this potential is not invariant under the "conventional" CP -symmetry or, in general, under any other CP2. Nevertheless, the model is still CPconserving because there exists at least one GCP (namely, CP4) which is a symmetry of the model. The fact that the potential has no CP2 symmetry is just irrelevant.
At this point, notice that there is no basis transformation that would make all the coefficients of the scalar potential V real [13]. Indeed, if it were possible to find such a real basis, then the potential would have an order-2 GCP. But such a symmetry is absent in the CP4 3HDM; therefore, the real basis does not exist. The absence of the real basis does not contradict explicit CP -conservation (and therefore, explicit T -conservation), because all basis-invariant combinations of the scalar couplings are CP -even. This model completely settles the issue of whether explicit CP conservation is equivalent to the existence of a real basis [31]: they are equivalent only for CP symmetries of order 2 and not for higher-order GCP.
As is conventional practice in building the models with extended Higgs sectors, it is necessary to require that the quartic parameters λ i are such that the potential is bounded from below (BFB). In other words, for quasiclassically large values of the Higgs fields along any direction in the scalar space, the potential must rise to plus infinity. One usually assumes the strong version of the BFB condition, which requires the quartic potential to strictly grow in any direction 2 . Necessary and sufficient conditions for BFB were obtained for the 2HDM earlier in Ref. [32], but no such deduction has hitherto been possible for the general 3HDM, or even for the CP4 3HDM we are dealing with. Nonetheless, we established in Appendix B a set of sufficient BFB conditions, which, although somewhat overly restrictive, will guarantee that the potential is indeed bounded from below. We will apply these conditions in our numerical analysis later on, which will therefore be a conservative one.

Extrema: generic solutions
When a scalar potential is explicitly CP -conserving, the vacuum of the theory can preserve CP , or spontaneously break it. The 2HDM, for instance, was first conceived as a model wherein spontaneous CPV has occurred [33]. The CP4 3HDM potential introduced in the previous section is explicitly CP conserving. Since we aim at extending the CP4 symmetry 1 In fact, in V1 one can write additional terms invariant under the same GCP transformation (2.4), which was indeed done in the previous publications on this model [13,14]. However, using the residual freedom of basis transformations which leave J invariant, one can simplify V1 to the form of Eq. (2.6). We provide an explanation of this procedure in Appendix A. 2 In principle, potentials stable in a weak sense, in which flat directions of their quartic potential are protected by the growing quadratic terms, are also acceptable. However, they correspond to measure zero regions in the parameter space, and we can avoid them in the phenomenological analysis.
to the Yukawa sector, it must be spontaneously broken; otherwise, the model would feature mass-degenerate fermions [14]. Without lack of generality, one can write the most generic charge-preserving vevs as where v 1 > 0, u ≡ v 2 2 + v 2 3 and we used the standard notation c ψ ≡ cos ψ, s ψ ≡ sin ψ. Later on, we will also use t ψ ≡ tan ψ. We then expand the doublets around the extremum using the following conventions: (2.8) By substituting this expansion in the Higgs potential V and setting the coefficients of the linear terms to zero we obtain the minimisation equations, also known as the tadpole conditions. The coefficient of the linear term in a 1 gives us the following relation: Let us for the moment consider the generic situation with sin(2ψ) = 0. It leads to γ 3 = −γ 2 ≡ −γ. The tadpole conditions for a 2 and a 3 then produce an additional relation, For given λ 8 and λ 9 , this equation relates the phase γ with the angle ψ. In order to simplify the analysis, we find it convenient to switch now to the real vev basis by rephasing φ 2 → φ 2 e −iγ and φ 3 → φ 3 e iγ . With this basis transformation, all real parameters in the potential stay the same, while λ 8 and λ 9 are rephased so that the quantity 2 arg(λ 9 )−arg(λ 8 ) remains unchanged. In the real vev basis, the tadpole condition above is written as Note that the latter can not be considered as an expression for ψ in terms of the parameters of the original potential since the phases of λ 8 and λ 9 depend now on γ. However, if Im (λ 8 ) is a free parameter and if ψ is known, one can deduce Im (λ 9 ). The tadpoles for h i lead to the following three relations: where we used the following shorthand notations: In total, therefore, we have five minimisation conditions to solve. One of them produces a relationship between the phases of the vevs of φ 2 and φ 3 , the remaining four must still be solved. As is usually the case in multi-Higgs models, obtaining the analytical expressions for the vevs (and their phases) in terms of the potential's parameters is exceedingly difficult. It is much simpler, in a numerical study, to take the values of the vevs and their phases as inputs. Therefore, we consider ψ as a free parameter and use Eq. (2.14) and v 2 1 + u 2 ≡ v = 246 GeV to determine v 2 1 and u 2 in terms of ψ and the scalar quartic couplings. Then, Eqs. (2.12) and (2.13) may be used to extract m 2 11 and m 2 22 . We will follow a similar procedure bearing also in mind our desire to find an acceptable scalar (and fermionic) mass spectrum.
If ψ is considered as input in a numerical scan over parameter space of the model, one needs to specify its range. Here, we argue that 0 ≤ ψ ≤ π/2, which corresponds to positive v 2 and v 3 , faithfully covers the entire set of relevant cases. The arguments go as follows. When the discrete symmetry CP4 is spontaneously broken, the potential has four degenerate minima, all related by the broken symmetry transformations but all corresponding to the same physics. These four minima are obtained by consecutively applying the transformation v 2 ↔ v 3 , γ → γ + π/2, while the real v 1 stays unchanged. Now, suppose we allow the free parameter ψ to take any value, thus allowing v 2 and v 3 to be either positive or negative. By the previous argument, we immediately see that , which, in turn, corresponds to the same model as (v 1 , |v 3 |, |v 2 |). Therefore, whatever value ψ takes, the model it leads to can also be found in the first quadrant of ψ.

Extrema: special points
Let us now consider two special values of ψ. The first case is when s 2ψ = 0, and without loss of generality, we can set ψ = 0 meaning v 3 = 0. In this case, the tadpole condition (2.9) is also satisfied. The other tadpole conditions get simplified, and after some algebra we arrive at the following relations: These relations are exactly those which we would get from the previous subsection in the limit s 2ψ → 0. Thus, we do not need to include this point as a separate case; we simply allow ψ to start from zero. The second singular point is c 2ψ = 0, implying v 2 = v 3 . Then, the tadpole condition (2.10) can be satisfied, in the real vev basis, either when u = 0 or when λ 8 is real, while λ 9 is purely imaginary. The former option would lead to an unphysical fermion spectrum, as we have already mentioned, while in the latter case the symmetry content of the model increases and involves now several accidental symmetries including the CP symmetry of order 2. This model, not being CP4-driven, falls beyond the scope of the present study. However the points in the vicinity of this limit, c 2ψ 1, are acceptable and, when accompanied with correspondingly small Re (λ 9 ), can potentially lead to realistic models.

Scalar mass matrices
In the real vev basis, we can expand the potential around the chosen extremum up to quadratic terms and construct the charged and neutral scalar mass matrices. For charged scalar fields, we get terms of the form Here, for simplicity, the dots indicate the duplicated entries of this symmetric matrix. One can explicitly verify that the charged Goldstone boson, which corresponds to the combination of fields given by is an eigenvector of this matrix with zero eigenvalue, as expected. The masses of the physical charged Higgs bosons can be explicitly calculated from the traces of powers of this matrix as usual, For the neutral scalars, the fact that CP is spontaneously broken by the vacuum implies the presence of a mixing between the real and imaginary neutral components of the Higgs doublets, h i and a j . Thus, the resulting 6 × 6 mass matrix can be written via symmetric 3-by-3 blocks, Finally, the h/a mixing terms depend on the imaginary parts of λ 8 and λ 9 , as could be expected. Using Eq. (2.11), we can parametrize these coefficients as Then, the mixing terms can be grouped as −λ 89 u 2 (s ψ h 2 −c ψ h 3 )(s ψ a 2 −c ψ a 3 ), which implies that the off-diagonal block in Eq. (2.21) is given by (2.25) One can again explicitly check that the neutral Goldstone boson, given by the following combination of fields is an eigenvector of the neutral mass matrix M with a zero eigenvalue. The resulting five neutral physical Higgs bosons mix and do not possess definite CP -properties, as should be expected considering that the CP symmetry has been broken. In particular, none of the neutral scalar states possess the parities or CP -charges defined in Refs. [13,14] for the CP4 unbroken case. As shown in Refs. [13,14], in the case of a vacuum with unbroken CP4, i.e. with u = 0, the physical Higgs spectrum organizes itself in mass-degenerate pairs, thus resembling that of the 2HDM. In such a vacuum, the model has a pair of mass-degenerate charged scalars, one SM-like neutral Higgs and two pairs of mass-degenerate neutral scalars, with masses m and M (one is an analogue of the heavier CP -even scalar in the 2HDM, H, the other is an analogue of the 2HDM pseudoscalar state, A). The spontaneous breaking of CP4 will induce a splitting in the mass spectrum, which vanishes in the u → 0 limit, assuming that the quartic couplings remain fixed.

Scalar alignment limit
An important property of viable multi-Higgs models in the post LHC era is that they should have regions in their parameter space for which one of their neutral scalars possesses a mass of about 125 GeV and closely resembles the SM Higgs boson. By this statement we mean that the couplings of this scalar to the SM gauge bosons and fermions must be very similar in magnitude to the corresponding SM values for those couplings. This may be achieved either in a "natural" way via symmetries (like in the case of the inert 2HDM [34,35]), or by a fine-tuning of the multi-Higgs model considered. In many scenarios, the desired region of parameter space can arise more or less "naturally" in a decoupling regime, wherein the extra scalar states are much heavier than the lightest SM-like one (see Ref. [36] for its introduction to the 2HDM). An intermediate case is the alignment limit, where some extra scalar masses may be low if certain relations by the couplings are satisfied. Again, within the 2HDM, the alignment without decoupling, and possible symmetry-based pathways to it, were considered in [37][38][39][40] and for maximally symmetric models beyond two doublets [41].
Let us examine how exact scalar alignment may arise in our model. In the original basis, the vevs of the doublets are given by Eq. (2.7), which we now rewrite as follows The Higgs basis is defined as a basis in which only the first doublet gets a vev. In that basis, we write the neutral complex fields (lower components of the Higgs doublets) as In general, the fields ρ i and η i are not mass eigenstates, and the neutral mass matrix includes mixing terms between them. The exact scalar alignment refers to the situation where one of the states, e.g. ρ 1 , is a mass eigenstate, identified with the 125 GeV Higgs boson H 125 . If this is the case, then its tree-level couplings to the W and Z bosons and to fermions are exactly the same as in the SM. Therefore, in such alignment limit the SM-like state will not mediate any tree-level FCNCs. The other (heavier) neutral scalars can, and generally do, have tree-level FCNCs, as will be discussed in Section 3.
To find the condition for scalar alignment, it is necessary to study the neutral scalar mass matrix in the Higgs basis. To do so, we first remark that the Higgs basis is not uniquely defined. In the 3HDM, any unitary transformation between Φ 2 and Φ 3 in Eq. (2.28) preserves the definition of the Higgs basis 3 . We make the following traditional (and convenient) choice for the Higgs basis: If one starts in the real vev basis, one should just set the phase γ to zero. For completeness, we give in Appendix C all expressions in the original basis with a non-zero γ.
To proceed to the physical scalars, we switch to the 6-dimensional real scalar field space and consider the 6 × 6 neutral mass matrix M, see Eq. (2.21). We order the fields in the real vev basis as ϕ a = (h 1 , h 2 , h 3 , a 1 , a 2 , a 3 ) and in the chosen Higgs basis as Φ a = (G 0 , ρ 1 , ρ 2 , ρ 3 , η 2 , η 3 ). The rotation to the Higgs basis is done by Φ a = P ab ϕ b , where the 6 × 6 matrix P ab is given by Eq. (C.6). To obtain the form of the neutral mass matrix in the Higgs basis, we start from M from Eq. (2.21) and use the rotation matrix P so that M H = P MP T . With this rotation, one immediately finds that the first column and the first row of M H only contain zeros, by virtue of the neutral Goldstone decoupling.
The exact Higgs alignment happens when the second row and column also decouple from the rest, i.e.
Therefore, in order to establish the scalar alignment conditions, one needs to calculate the second row of M H and to set its off-diagonal elements to zero. We did that and observed that M H 23 = M H 25 = M H 26 = 0 are automatically fulfilled for generic vevs v 1 , u and angle ψ.
The only non-zero entries of the second row can be written, after some algebra, in a very compact form: Exact scalar alignment is achieved when M H 24 = 0, which implies with arbitrary Λ abcd constrained only by BFB conditions. Then this potential automatically incorporates the exact scalar alignment, which can be verified by direct differentiation. One scalar mass eigenstate is always aligned with the direction of vevs, and its mass squared is 2m 2 . The non-trivial result of the above exercise is that, within CP4 3HDM with fully spontaneously broken CP4 symmetry, setting m 2 11 = m 2 22 is the only way to impose scalar alignment.

Yukawa models with CP4
The CP4 symmetry can be extended to the Yukawa sector, provided the CP4 transformation also mixes the fermion generations, as follows Such an extension has been performed in the framework of 2HDM in Refs. [25][26][27] and ran into difficulties with excessive accidental symmetries. In this work, we implement such an extension for the CP4 3HDM. The CP4 symmetry strongly constrains the Yukawa interaction matrices. In Ref. [14], some examples of such interactions were given under the simplifying assumption that the right-handed up and down fermions, as well as the left-handed doublets, transform in the same way, i.e. Y L = Y d = Y u . In this work, we lift this assumption and derive all possible forms of the Yukawa interaction matrices compatible with CP4 and not running into immediate conflict with experiment, that is, not leading to massless fermions or an insufficient mixing.
In this work, we focus on the quark sector only while the lepton sector can be incorporated in a similar way. The quark Yukawa Lagrangian in which we explicitly indicated the Higgs family index and omitted for clarity the quark flavor indices, is invariant under CP4 if and only if With the explicit expression for X given by Eq. (2.4), we get As usual, an appropriate change of the basis in the q L , u R , d R spaces can bring all the matrices Y to the block-diagonal form and we allow the parameters α L , α d , and α u to be all different. Here, we selected the third fermion to be a CP4 singlet but any other choice would be equivalent. As we have seen above, in order to properly define inequivalent generalized CP transformations on the fermion sector, the possible values of the phases in Eq. (3.5) must be such that the matrices Y * Y are of finite order 2 p−1 with a distinct p defining an inequivalent CP2 p symmetry. When solving Eqns. (3.4), we do not assume that α corresponds to the CP4 case. In fact, those equations only imply that the fermion bilinears coupled to φ 2 and φ 3 must faithfully transform under CP4, but the transformation law of the fermion fields individually may be of even higher order. We leave open this possibility and just solve the coupled equations.
We found that there are only four classes of Yukawa matrices Γ a and ∆ a satisfying Eq. (3.3) for some α's and not running into an immediate conflict with the data. We label them as cases A, B 1 , B 2 , B 3 . We describe them below for Γ's in terms of their independent complex parameters g ij , and then briefly comment on matrices ∆ a . Note that although for notational simplicity we use the same names for the parameters g ij for the scenarios A, B 1 , B 2 , B 3 , they should be regarded as different ones.
Before listing these results, one comment is in order. When solving all equations for α's and Γ's, we often obtain seemingly different solutions with extra minus factors in some rows or columns. We checked that in all cases they represent the same model and differ just by a sign flip in one or several fermion fields. The four cases shown below are the ones which cannot be reduced to one another.
• Case A. e iα L = 1 and e iα d = 1, giving It may seem that Γ 1 is not as generic as it would be in the CP -conserving version of the SM because the off-diagonal elements are related to each other. However, when α = 0, the CP-symmetry within the Yukawa sector is effectively of order 2. As a consequence, there exists a fermion basis in which the CP-transformation is the canonical one in the fermion sector and, in this basis, Γ 1 is simply an arbitrary real matrix.
• Case B 1 . e iα L = i and e iα d = 1, giving • Case B 2 . e iα L = 1 and e iα d = i, giving All parameters apart from g 33 can be complex in all cases. Notice also that in all cases the matrices Γ 2,3 are expressed in terms of the same complex parameters and have the same textures.
For the up-quark sector, we get the same structures for ∆'s. Indeed, the equations for ∆'s are the same as for Γ's with the exchange 2 ↔ 3. However, when constructing a viable model, we are not forced to select the same case for Γ's and ∆'s. We only must ensure that the transformation properties of the left-handed doublets (i.e. the values of α L ) are the same in both sectors. Therefore, we have two series of four different combinations each for the down and up quark sectors: Note, in the present analysis, we do not consider the case (A,A) which corresponds to a situation when the CP4 symmetry does not affect the fermion sector. The CKM matrix in this case is real at the tree level. However, the original CP symmetry is broken by the Higgs sector giving a plausible chance to generate CPV effects in the fermion sector by means of radiative corrections. Whether this mechanism is capable of producing the experimentally observed amount of CPV deserves a closer study that is left for future work.
By inspecting the textures of the Yukawa matrices in cases B 1 , B 2 , B 3 , one immediately sees that CP4 must be spontaneously broken. The unbroken case corresponds to the vev alignment (v, 0, 0), and in this case M d = Γ 1 v/ √ 2 will produce pathological quark masses in all cases B, by forcing some quarks to be either massless or mass-degenerate [14]. Also, the strength of CP4 breaking cannot be too weak, because it is the value of u multiplied by some elements of Γ 2,3 that drives the mass splitting across the two fermion generations. Fortunately, as was mentioned in the previous section, there is no reason to expect u v 1 in a generic situation.
Let us also remark that when solving matrix equations (3.4) we found other non-trivial solutions. We do not list them because they immediately enter in conflict with the quark masses/mixing properties. For example, there exists a solution with α L = α d = α u = π/4 with the following Yukawa textures (3.12) As mentioned above, the fact that the phases in Eq. (3.5) are π/4 implies that the CP symmetry acts on fermionic fields as CP8, and it only becomes CP4 for the fermion bilinears. However, when multiplied by vevs, these Yukawa matrices lead to block-diagonal M d and M u , both of which are diagonalized by rotations between the first and second fermion generations. This means that V CKM will have a block-diagonal form with only one mixing angle, which contradicts experimental results. Therefore, this solution and other similar cases are disregarded.

Reconstructing Yukawa matrices
Having picked up one of the above cases and having fixed the vevs v i , one can obtain the quark mass matrix. For example, for the down quarks we have Even if the Γ's have a given structure, they all collapse to a single matrix M d where this structure may not be present anymore. Of course, in texture-constrained situations in which the Yukawa matrices of different Higgs fields do not share common non-zero entries, one can reconstruct Γ a from M d . But in a generic situation with overlapping Yukawa matrices, this is impossible even knowing v i . The remarkable property of the CP4 3HDM is that, despite overlapping entries, it is possible to reconstruct all Γ's from M d once we know the vevs and which of the above cases A, B 1 , B 2 , B 3 we deal with. This is so because in all the cases there is a pair of entries of M d which is determined by two elements g ij only. Consider, for example, the case B 3 and the off-diagonal elements in the third column: (3.14) By inverting this matrix relation, one reconstructs the elements g 13 and g 23 in terms of (M d ) 13 and (M d ) 23 . The same holds for other elements. The resulting mass matrices for the models B 1 , B 2 and B 3 are constrained only by the following relations: However, in any of these cases, the matrices In CP4 3HDM, we see two remarkable properties: • all hermitean H d and H u are reachable for an appropriate choice of (Γ a , ∆ a ) and vevs; • the first step (Γ a , ∆ a ) → (M d , M u ) is invertible once vevs are knows.
If one wishes, one can perform the scan of the Yukawa parameters space starting with the physical quark masses, mixing, and CPV and then generating g ij elements which, by construction, will exactly reproduce the measured values. These features serve as an "existence proof" of good parameter space points, in the sense that we will be able to fit, within CP4 3HDM, all the fermion masses, mixing, and the amount of CP violation, to arbitrary precision. However, in the numerical scan to be described in the next section we will not employ this reverse engineering algorithm but rather stick to the standard procedure. While the reverse engineering algorithm is useful to prove the existence of a solution fitting the masses and mixing, it can very easily fall in a small region of the allowed parameters space. Since we will be interested in fitting additional observables, we must try to explore the largest region of the parameters space possible. We observed that the standard χ 2 minimization procedure, starting with a seed point, quickly converges to a nearby parameter space point which reproduces quark masses, mixing, and CPV arbitrarily well. The largest amount of the computer time in the scan was spent to find points which fit sufficiently well the meson oscillation observables.
The Yukawa matrices in cases B 1 (3.7) and B 2 (3.8) have 7 complex and one real free parameters, making in total 15 real free parameters in each sector. In case B 3 , we have in total 13 real free parameters per sector. We remind that in the SM, when we start from arbitrary complex matrices Γ and ∆, we have 18 initial free parameters in each sector. Certainly many of these free parameters are superficial and can be removed by basis changes in the left-handed and right-handed quark sectors. In CP4 3HDM, we have fewer superficial free parameters. For example, rotations of the right-handed quark fields do modify some observables, namely, the off-diagonal FCNC elements linked to extra neutral Higgs fields. However, these elements are, themselves, correlated, which renders counting the independent physical parameters not a straightforward exercise. We postpone this calculation to a future work.

Flavour-changing neutral currents
The generic presence of Higgs-exchange-induced large FCNCs at tree level is a notorious problem of all multi-Higgs-doublet models. For example, in the down quark sector, we get after electroweak symmetry breakinḡ When diagonalizing the quark mass matrix by a bi- The Yukawa interactions can be written as In the last line, we switched to physical quarks d ph. and introduced the following matrices The key observation is that the CP4-symmetric Yukawa structures Γ i generically produce matrices Γ (H,ph.) 2,3 with unsuppressed off-diagonal terms, typically of the same order as the off-diagonal elements of Γ 2,3 . Even if some of these off-diagonal elements are small due to accidental cancellations assisted by vevs, this cannot happen with all off-diagonal elements, again due to the specific structure of the CP4 3HDM Yukawa sector. Notice also that one cannot assume that entries g ij are very small in Γ 2,3 since those, when multiplied by v 2 and v 3 , must provide sufficient mass splitting for quarks. Thus, the only solution in sight is to assume that a strong or exact alignment takes place in the scalar sector, as discussed in Section 2.6. In this case, the field ρ 1 in Eq. (3.19) is identified with the 125 GeV Higgs boson H 125 and it is protected from tree-level FCNCs. However, dangerous FCNCs arise from the other neutral Higgs bosons' exchanges, unless they are sufficiently heavy or other cancellations take place. One then faces the necessity of a numerical study in order to check if the resulting models can pass the tight experimental constraints.

Numerical scan
In this Section, we shall describe the procedure of scanning of the parameter space of the model in the Higgs sector (scalar potential) and in the Yukawa sector, with different combinations of Yukawa matrices Γ a and ∆ a discussed above. Our intent is to show that the CP4-3HDM scenarios can realistically fit the current data avoiding large departures from the SM predictions in well-measured observables.
The Yukawa sector scan will take as input the vevs and the diagonalizing rotation matrices of the scalar fields. Therefore, we first analyze the scalar sector which, at least at tree level, requires no information from the fermionic sector. Then, we perform an efficient numerical scan over the Yukawa sector parameters. At each run, we start from a seed point, converge to a point in the Yukawa parameter space that gives a good description of quark masses and the CKM matrix, and then check K-meson parameters K and ∆m K . Keeping only those points which are sufficiently close to their experimental values, we proceed to the B-meson oscillation parameters. At the end, we have a fair selection of good points, which pass all the criteria we have imposed. Let us describe the numerical procedure as well as the main results in detail.

Parameter space of the scalar sector
As we have have seen in Section 2.6, we can guarantee the presence of a SM-like Higgs field H 125 by working in the scalar alignment limit. In the present analysis, we shall take this limit and explore the corresponding predictions.
The basic procedure adopted in our exploration of the parameter space can be summarized as follows: • The alignment limit, i.e. m 2 11 = m 2 22 , has been imposed in the mass matrices (2.17) and (2.21). It can be formulated as In this limit, as shown in Section 2.6, the SM-like Higgs mass has a simple analytic expression, Eq. (2.33). We can then use this expression to eliminate one extra parameter, i.e.
• At this stage we have 9 free parameters: λ 2,3,4,5 , λ 3,4 Im(λ 8 ), and two vevs ratios, which can be expressed in terms of angles ψ and β. We choose an initial random point satisfying the BFB conditions, presented in Eq. (B.4). In order to stay in the perturbative regime we took 4 |λ ( ) For the vevs we allowed the ratios u/v 1 = t β and v 3 /v 2 = t ψ to be between 0 and 100.
• We now use this point to compute the mass eigenstates and the respective field rotations. We use iminuit 1. 2 [43] in order to find the values of the free parameters that minimize the function where m lightest is the numerical value of the mass for the lightest massive neutral scalar, and σ Higgs is the allowed standard deviation. Since one of the mass eigenstates is H 125 , whose mass by construction is set to 125 GeV, this χ 2 scalar -minimization promotes the lightest massive neutral state to be the observed SM-like Higgs boson. We then take an extremely small deviation σ Higgs = 10 −10 GeV 2 , such that the method can quickly converge to the scenario where the massive state that is aligned with the Higgs basis is the lightest massive neutral scalar with mass of 125 GeV. Guaranteeing the existence of this state automatically ensures that the neutral scalar Hessian matrix is non-negative. We also make sure that the resulting points of the scalar parameter space produce positive mass squared for the charged Higgs bosons.
The assumption that the 125 GeV Higgs is the lightest one is not obligatory. Whether CP4 3HDM can accommodate lighter Higgs states and not run into an immediate conflict with the experimental data deserves a dedicated study, which we delegate to a future work.
• Once the minimization of χ 2 scalar is achieved, we check if the final values of the scalar potential couplings satisfy the BFB conditions. Also, as a safety check, we verify if all the fields' masses squared are positive and if the lightest neutral massive one is actually the SM-like Higgs boson.
The presence of additional neutral and charged scalars can lead to large deviations in the well measured electroweak observables. The most constraining electroweak precision observables can be summarized into the three oblique parameters S, T, U [44][45][46]. We have computed these parameters, as well as three additional ones V , W and X (with no relevant impact on the parameter space), using the corresponding expressions found in Ref. [47]. The current electroweak precision measurements lead to [48] S = 0.05 ± 0.10 , T = 0.08 ± 0.12 , U = 0.02 ± 0.10 . (4.5) These parameters are strongly correlated: there is a 91% correlation between S and T parameters, while the U parameter is −61% and −82% anti-correlated with S and T , respectively. In Fig. 1 we present the T S-and U T -plane where the third oblique parameter is set to its best fit value. We can see that both S and U place no constraints on the parameter space, while T restricts the allowed parameters space significantly. The points surviving the 3σ constraints from the ST U oblique parameters are represented by black dots while the remaining points are shown by gray dots. This representation will be used throughout this section. S with the ellipses drawn for 90%, 95% and 99% confidence level with U fixed at its the best fit value. Right: U vs. T with the ellipses drawn for 90%, 95% and 99% confidence level with S fixed at its the best fit value.
We can now explore some trends for the vevs and the scalar spectrum which emerge in the CP4 3HDM scan. In Fig. 2 we show the region of the vev ratios for the points which pass the scalar sector constraints discussed above. Both ratios tend to be O(1), even though we initially allowed for much wider ranges. The shape of this region arises from an interplay of several factors: the BFB and perturbativity constraints, the scalar alignment assumption Eq. (2.32), the requirement that H 125 be the lightest massive neutral scalar, and finally the relation among the vev ratios and λ's encoded in Eq. (2.14). Qualitatively, choosing the vev ratios very far from unity would push some of the quartic couplings beyond the allowed range. If the scalar alignment requirement or the constraint on the lightest massive neutral scalar were dropped, larger regions on this plot would be accessible. In Fig. 3 (left) we present the masses of the two charged scalars H + 1 and H + 2 . One sees that the mass of lightest charged scalar can go down to as low as 90 GeV and still be compatible with the electroweak precision data. We label the five neutral physical scalars  These states are mass ordered. In Fig. 3, right, we plot the next-to-lightest neutral scalar mass m H 3 vs. the ratio of two vevs. Again, the same constraints that were shaping Fig. 2 are at work here. In addition, one can directly compute the sum of all neutral masses squared as i m 2 H i = TrM n . The trace of the neutral scalar mass matrix (2.21) is given in terms of v 2 , vev ratios, and quartic couplings, and therefore it is bounded from above. Satisfying all scalar sector constraints maximizes i m 2 H i when all vevs are of the same order. Since the scalars (4.6) are mass-ordered, the largest value of m H 3 is attained when all additional neutral Higgses are equally heavy, which explains the rather sharp upper boundary of Fig. 3, right.
The mass spectra of yet heavier Higgses H 4 , H 5 , H 6 are shown in Fig. 4. Here, the upper boundaries pass higher and are less pronounced. For example, maximizing M H 5 implies keeping H 3 and H 4 close to H 125 and setting M H 5 = M H 6 .

Parameter space of the Yukawa sector
Once the scalar sector analysis is done, we have a numerical set of vevs and scalar rotation matrices, which can be used as input for the Yukawa sector scan. In this sector, we are interested in finding the parameter space points where the experimentally known masses of the quarks as well as the CKM mixing matrix can be fully accommodated. We used the following experimental data [48]: for the quark masses, sin θ 12 = 0.22496 ± 0.001 , sin θ 23 = 0.0416496 ± 0.001 , sin θ 13 = 0.00361726 ± 0.0001 (4.8) for the CKM mixing angles, and sin δ = 0.949 ± 0.01 (4.9) for the CP violating phase.
On top of fitting the masses and mixings in the quark sector, we also need to guarantee that we do not have too large FCNCs. The meson oscillation observables set strong constraints on these models. For the K-meson sector, we have the CPV observable K and the mass difference ∆m K between K L and K S . The observables read [49] with where V is the CKM mixing matrix. The loop functions S 0 (x) and S 0 (x, y) are given in Ref. [50] (see also Ref. [49]). Finally, ϕ = 43.51 • and κ = 0.94 were taken from Ref. [48,51,52].
In the above expression, the function ∆S(K) is the only model dependent contribution. For our model, the contribution of the neutral scalar fields is given by where the functions f XY (µ H ) encode the information on the Wilson coefficients and hadronic matrix elements. Their explicit form, as well as numerical values, can be found in Ref. [49]. The relevant part here is the model dependent ∆ sd X (H i ), which in our model reads We have included in the χ 2 -fit the K-meson sector, with the experimental values for the observables [48] | K | exp = 2.228 × 10 −3 , ∆m exp K = 3.5 × 10 −15 GeV , (4.14) and we allowed a 50% deviation from the experimental measurement. The scan procedure is very similar to the one adopted in the scalar sector, with the χ 2 function to be now minimized defined as (4.15) Here, m i are the quark masses, s θ i correspond to the sine of mixing angles and CPV phase. In the numerical scan, the absolute values of the Yukawa parameters were taken in the range [0, 5] with arbitrary phases. Once the numerical quark mass matrices are found, we extract the mixing angles and the CPV phase through the reshaping invariants (4.16) For the B-meson sector, the mass difference can also set strong constraints on the pa-  rameter space. From Ref. [49] we have  For the three models surviving the meson oscillations constraints despite non-trivial FCNCs in the down-quark sector, we looked at the magnitude of the off-diagonal Yukawa entries in the mass basis. These are shown in Fig. 8. Namely, we computed the Yukawa matrices in the Higgs basis and for quark mass eigenstates Γ (H,ph.) 2,3 , which were defined in Eq. 3.20, and extracted the off-diagonal entry in either of these two matrices with largest absolute value, denoted as Max(|Γ mass ij |). We repeated the same for the up-quark sector extracting largest off-diagonal element from ∆ (H,ph.) 2,3 , denoted as Max(|∆ mass ij |). Then, for each point of the scan, we plotted these two values to get the scattering plots Fig. 8. These plots exhibit certain clustering around √ 2m b /v and √ 2m t /v for down and up sectors which is natural when FCNCs are unsuppressed. Most of these points, however, do not pass the tight meson oscillation constraints. Those which do are shown in red. For the down-quark sector, the largest absolute value for the off-diagonal entry typically stay below 10 −3 , although a few outliers exist. In the up-quark sector, the off-diagonal elements can indeed be large.
What emerges from this analysis is that the magnitudes of the off-diagonal elements do not necessarily need to be extremely small to provide a good description of the data. Several neutral Higgs bosons, which mediate FCNCs, can interfere destructively due to the intrinsic structural properties of the model, and lead to acceptable results. This is also reminiscent of what was observed in certain classes of the 2HDM, where the real and imaginary parts of the complex neutral Higgs field destructively interfere and bring FCNCs under control, see e.g. Ref. [56].
In order to illustrate the typical textures which arise in this scan, we give in Table 1  In both cases, the value of χ 2 yuk. given by (4.15) is very small, about 0.01, which indicates that the model is capable of describing the quark masses, mixing, and CPV to arbitrary precision, as discussed in Section 3.2.

Discussion and conclusions
In this paper, we undertook the first step in the phenomenological investigation of a unique variant of 3HDM: the simplest model which realizes the CP symmetry of order 4 (CP4) without any accidental symmetry. The single symmetry CP4 strongly constrains both the scalar and Yukawa sectors, but it is nevertheless capable of accurately fitting all quark masses and mixing parameters, including CP violation, as well as the SM-like Higgs boson properties.
When developing the model, we obtained convenient sufficient conditions for the stability and boundedness from below for the potential, analytic expressions for the minima and scalar mass matrices, and the condition for alignment in the scalar sector. We also established all forms of the CP4-invariant Yukawa Lagrangian with four distinct Yukawa textures in the up and down quark sectors.
The spontaneous breaking of the CP4 symmetry induces proper splittings in the fermion mass spectra consistent with observations. In the scalar alignment limit, the SMlike Higgs boson H 125 , which we assumed to be the lightest scalar, was shown to conform to the experimental constraints on the FCNCs and the electroweak precision observables. The heavier Higgs partners do induce tree-level FCNC effects, but, thanks to a destructive interference among them, the resulting deviations from the SM values of kaon and B-meson mixing and CPV are kept under control.
To make the analysis quantitative, we have set up a conservative scan over the parameter space of the model adopting the most relevant theoretical and phenomenological constraints from the scalar and Yukawa sectors. Our χ 2 fit has revealed viable parameter space regions for several distinct scenarios for the combinations of up and down Yukawa coupling matrices and showed a promising potential to rule out some of them. The benchmark points found in this analysis could be very useful for follow-up in-depth phenomenological explorations of the CP4 3HDM in the future.
In this exploratory study we did not address many specific issues which merit close examination. These directions of research include: • checking the compatibility of the points which pass meson oscillation constraints (red dots in the above plots) with direct searches for new neutral and heavy Higgs scalars  Table 1. Two benchmark points for our model implementation in cases (B 1 , B 1 ) and (B 1 , B 3 ) and their resulting mass spectra and predictions for the flavor observables. at the LHC. Given that the masses of extra Higgses are rather low, one may suspect that these points are already ruled out by the negative results of the heavy Higgs boson searches. This is not the case since, in the scalar alignment regime, all extra Higgses decouple from the gauge bosons. In addition, ttH i coupling-and therefore the gg-fusion production cross section-may well be suppressed with respect to the SM, but even if it is not, the extra Higgses can easily avoid the main LHC searches which focus on the ZZ/W W or γγ final states. One needs to look into the qq decay channels, especially into flavor-violating ones, to see their manifestation. In any case, a dedicated study is required to check what percentage of the points we have found are ruled out by the direct LHC searches.
• checking whether the light charged Higgs bosons can satisfy B-decay constraints, especially b → sγ. It is known that, within 2HDM type-II, this decay places a very stringent lower bound on the charged Higgs mass m H + ∼ > 580 GeV, [57]. In our case, there are two charged Higgs bosons which may destructively interfere due to the CP4-driven structural relations between the Yukawa matrices for Φ 2 and Φ 3 . We expect that, at least for some of the scan points, this interference will significantly reduce the combined charged Higgs contribution to b → sγ and relax the constraints.
• checking the magnitude of the electric dipole moments (EDMs) induced within CP4 3HDM and its compatibility with the negative results of the electron and neutron EMD searches, see e.g. [58][59][60]  • investigating how close the CP4 3HDM can approach the SM and which observables exhibit the strongest deviations. In multi-Higgs models such as 2HDM, it is customary to define the SM-like limit of the theory. In CP4 3HDM, it is not guaranteed that this limit exists at all, again due to the structural relations among various parameters imposed by CP4 symmetry. In our scan, the SM-limit was definitely absence due to the scalar alignment limit condition we imposed, m 2 11 = m 2 22 .
It remains an open question if such limit exists if m 2 22 is taken to be negative and large.
• checking whether the CP4 3HDM can alleviate any of the several tensions persistently observed over the last few years in various B-meson decays; • checking whether the particular variants of the CP4 3HDM such as the Yukawa combination (A, A) or the presence of new light scalars, which we briefly mentioned in the text, are ruled out or still compatible with the data; • incorporating the charged leptons and verifying that the experimental bounds on lepton flavor violating processes can be satisfied; • building neutrino mass models based on CP4 symmetry. A first step in this direction was take in [61].
• verifying if the CP4 3HDM can say something useful about the strong CP problem.
The key message of this work is that, unlike in many other multi-Higgs models, a single symmetry requirement strongly shapes the model and brings us very far without falling into a direct conflict with experiment. It is true that on the way we had to stick to the exact scalar alignment, whose condition m 2 11 = m 2 22 is not symmetry-protected. However, it may turn out that the CP4 3HDM emerges as the lower energy limit of yet another model with a higher symmetry. In this case, m 2 11 = m 2 22 would arise naturally from a higher symmetry requirement at a certain energy scale, and then, when running down to the electroweak scale, the parameters would not deviate much from the alignment condition. Again, verifying that this construction is viable requires a dedicated study.
Finally, one can imagine multi-Higgs models with CP symmetries of even higher order: CP8, CP16, etc. These will require more than three Higgs doublets since the full classification of symmetries possible within the 3HDM is complete [15,29,30] and it does not include such symmetries. Whether these can be extended to the Yukawa sector in a satisfactory manner remains to be seen.
In short, models based on a single higher-order CP symmetry arise as an attractive minimalistic framework, with many phenomena yet to be explored.
for of V 1 analyzed in [13,14]: 2) Using rephasing of φ 2 and φ 3 , both λ 5 and λ 6 can be made real, but then λ 8 , λ 9 are in general complex. It turns out that this form of the CP4 3HDM potential is not the simplest one. Indeed, there exists a residual freedom of basis changes φ i → φ i = R ij φ j , which leaves invariant the symmetry transformation J if RXR T = X. This condition enforces R to be of the block-diagonal form with R 11 = 1 and the 2 × 2 block being a generic Sp(2, C) transformation. Absorbing the rephasings into redefinition of φ's, one gets one extra transformation freedom: an orthogonal rotation between doublets φ 2 and φ 3 .
This rotation reparametrizes the general CP4 3HDM potential. The potential remains of the same form but new parameters are expressed via the old parameters in a non-trivial way. In particular, it induces an orthogonal transformation between λ 5 and λ 6 with twice the (φ 2 , φ 3 ) rotation angle. This is so because, in group-theoretic terms, there exist two bilinear combinations which pick up i under CP4: φ † 1 φ 2 + φ † 3 φ 1 and φ † 1 φ 3 − φ † 2 φ 1 , and it is the product of these two bilinears with their conjugates that is encoded in the λ 5 and λ 6 terms.
Therefore, starting from V 1 given by Eq. (A.2), one can always find a basis in which λ 5 = 0 or λ 6 = 0. In the present work, we choose the latter option and eliminate λ 6 . We could have used the former one; in fact, that choice would be preferred for investigation of 3HDM with unbroken CP4. We stress that when setting one of these parameters to zero we do not lose any generality.

B Boundedness from below
In the 2HDM, the exact necessary and sufficient boundedness-from-below (BFB) conditions are known explicitly for the most general renormalizable scalar potential [32]. Beyond two doublets, this general result is still unknown. Although in certain simplified or symmetric cases one can establish them, we could not find such necessary and sufficient conditions for the CP4 3HDM.
For the purpose of a numerical scan, we limited ourselves to a set of sufficient conditions, which constrain the parameters somewhat stronger than needed but which guarantee that the models selected have stable potentials. Here we outline how these conditions were derived.

C Three scalar bases
Here we define the three bases for the scalar fields, appropriate for the discussion of the scalar alignment: the original basis, one specific choice of the Higgs basis, and the physical scalar basis.
Finally, the physical scalar basis, which is defined as the basis in which the neutral mass matrix M n is diagonal, must be written via real fields, which we label as In 3HDM, the Higgs basis is not uniquely defined: one can always perform a unitary transformation between Φ 2 and Φ 3 in Eq. (C.2) preserving the definition of the Higgs basis. Here, we use the following convenient and traditional choice: The physical scalar basis is linked with the Higgs basis and the original basis via matrices Q and R, respectively: Notice that since the Goldstone does not mix with the physical Higgses, the matrix Q is of block-diagonal form: