Neutrino mixing and leptonic CP violation from S4 flavour and generalised CP symmetries

We consider a class of models of neutrino mixing with S4 lepton flavour symmetry combined with a generalised CP symmetry, which are broken to residual Z2 and Z2 × HCPν symmetries in the charged lepton and neutrino sectors, respectively, HCPν being a remnant CP symmetry of the neutrino Majorana mass term. In this set-up the neutrino mixing angles and CP violation (CPV) phases of the neutrino mixing matrix depend on three real parameters — two angles and a phase. We classify all phenomenologically viable mixing patterns and derive predictions for the Dirac and Majorana CPV phases. Further, we use the results obtained on the neutrino mixing angles and leptonic CPV phases to derive predictions for the effective Majorana mass in neutrinoless double beta decay.


Introduction
Understanding the origin of the pattern of neutrino mixing that emerged from the neutrino oscillation data in the recent years (see, e.g., [1]) is one of the most challenging problems in neutrino physics. It is part of the more general fundamental problem in particle physics of understanding the origins of flavour in the quark and lepton sectors, i.e., of the patterns of quark masses and mixing, and of the charged lepton and neutrino masses and of neutrino mixing.
The idea of extending the Standard Model (SM) with a non-Abelian discrete flavour symmetry has been widely exploited in attempts to make progress towards the understanding the origin(s) of flavour (for reviews on the discrete symmetry approach to the flavour problem see, e.g., [2][3][4]). In this approach it is assumed that at a certain high-energy scale the theory possesses a flavour symmetry, which is broken at lower energies to residual JHEP12(2017)022 symmetries of the charged lepton and neutrino sectors, yielding certain predictions for the values of, and/or correlations between, the low-energy neutrino mixing parameters. In the reference 3-neutrino mixing scheme we are going to consider in what follows (see, e.g., [1]), i) the values of certain pairs of, or of all three, neutrino mixing angles are predicted to be correlated, and/or ii) there is a correlation between the value of the Dirac CP violation (CPV) phase δ in the neutrino mixing matrix and the values of the three neutrino mixing angles, 1 θ 12 , θ 13 and θ 23 , which includes also symmetry dependent fixed parameter(s) (see, e.g., [5][6][7][8][9][10][11][12][13] and references quoted therein). These correlations are usually referred to as "neutrino mixing sum rules". As we have already indicated, the sum rules for the Dirac phase δ, in particular, depend on the underlying symmetry form of the PMNS matrix [5][6][7][8][9] (see also, e.g., [10][11][12][13]), which in turn is determined by the assumed lepton favour symmetry that typically has to be broken, and by the residual unbroken symmetries in the charged lepton and neutrino sectors (see, e.g., [2-4, 7, 9]). They can be tested experimentally (see, e.g., [6,10,[14][15][16]). These tests can provide unique information about the possible existence of a new fundamental symmetry in the lepton sector, which determines the pattern of neutrino mixing [5]. Sufficiently precise experimental data on the neutrino mixing angles and on the Dirac CPV phase can also be used to distinguish between different possible underlying flavour symmetries leading to viable patters of neutrino mixing.
While in the discrete flavour symmetry approach at least some of the neutrino mixing angles and/or the Dirac phase are determined (directly or indirectly via a sum rule) by the flavour symmetry, the Majorana CPV phases α 21 and α 31 [17] remain unconstrained. The values of the Majorana CPV phases are instead constrained to lie in certain narrow intervals, or are predicted, in theories which in addition to a flavour symmetry possess at a certain high-energy scale a generalised CP (GCP) symmetry [18]. The GCP symmetry should be implemented in a theory based on a discrete flavour symmetry in a way that is consistent with the flavour symmetry [19,20]. At low energies the GCP symmetry is broken, in general, to residual CP symmetries of the charged lepton and neutrino sectors.
In the scenarios involving a GCP symmetry, which were most widely explored so far (see, e.g., [19,[21][22][23][24][25]), a non-Abelian flavour symmetry G f consistently combined with a GCP symmetry H CP is broken to residual Abelian symmetries G e = Z n , n > 2, or Z m ×Z k , m, k ≥ 2, and G ν = Z 2 ×H ν CP of the charged lepton and neutrino mass terms, respectively. 2 The factor H ν CP in G ν stands for a remnant GCP symmetry of the neutrino mass term. In such a set-up, G e fixes completely the form of the unitary matrix U e which diagonalises the product M e M † e and enters into the expression of the PMNS matrix, M e being the charged lepton mass matrix (in the charged lepton mass term written in the left-right convention). At the same time, G ν fixes the unitary matrix U ν , diagonalising the neutrino Majorana mass matrix M ν up to a single free real parameter -a rotation angle θ ν . Given the fact that the PMNS neutrino mixing matrix U PMNS is given by the product  Table 1. The best fit values, 2σ and 3σ ranges of the neutrino oscillation parameters obtained in the global analysis of the neutrino oscillation data performed in [26].
all three neutrino mixing angles are expressed in terms of this rotation angle. In this class of models one obtains specific correlations between the values of the three neutrino mixing angles, while the leptonic CPV phases are typically predicted to be exactly 0 or π, or else π/2 or 3π/2. For example, in the set-up considered in [19] (see also [21]), based on G f H CP = S 4 H CP broken to G e = Z T 3 and G ν = Z S 2 × H ν CP with H ν CP = {U, SU }, 3 the authors find: sin 2 θ 13 = 2 3 sin 2 θ ν , sin 2 θ 12 = 1 2 + cos 2θ ν = It follows, in particular, from the results on the neutrino oscillation parameters -best fit values, 2σ and 3σ allowed ranges -obtained in the latest global fit of neutrino oscillation data [26] and summarised in table 1, to be used in our further analysis, 4 that the predictions quoted in eq. (1.2) for sin 2 θ 12 and sin 2 θ 23 lie outside of their respective currently allowed 2σ ranges. 5 Another example of one-parametric models is the extensive study performed in [28], in which the authors have considered two different residual symmetry patterns. The first pattern is the one described above, and the second pattern has G e = Z 2 × H e CP and G ν = Z 2 × Z 2 × H ν CP as residual symmetries in the charged lepton and neutrino sectors, respectively. The authors have performed an exhaustive scan over discrete groups of order less than 2000, which admit faithful 3-dimensional irreducible representations, and classified phenomenologically viable mixing patterns. 3 S, T and U are the generators of S4 in the basis for its 3-dimensional representation we employ in this work (see subsection 3.2). 4 The results on the neutrino oscillation parameters obtained in the global fit performed in [27] differ somewhat from, but are compatible at 1σ confidence level (C.L.) with, those found in [26] and given in table 1. 5 We have used the best fit value of sin 2 θ13 to obtain the prediction of sin 2 θ12 leading to the quoted conclusion. Using the 2σ allowed range for sin 2 θ13 leads to a minimal value of sin 2 θ12 = 0.340, which is above the maximal allowed value of sin 2 θ12 at 2σ C.L., but inside its 3σ range.

JHEP12(2017)022
Theoretical models based on the approach to neutrino mixing that combines discrete symmetries and GCP invariance, in which the neutrino mixing angles and the leptonic CPV phases are functions of two or three parameters have also been considered in the literature (see, e.g., [29][30][31][32]). In these models the residual symmetry G e of the charged lepton mass term is typically assumed to be a Z 2 symmetry or to be fully broken. In spite of the larger number of parameters in terms of which the neutrino mixing angles and the leptonic CPV phases are expressed, the values of the CPV phases are still predicted to be correlated with the values of the three neutrino mixing angles. A set-up with G e = Z 2 × H e CP and G ν = Z 2 × H ν CP has been considered in [32]. The resulting PMNS matrix in such a scheme depends on two free real parameters -two angles θ ν and θ e . The authors have obtained several phenomenologically viable neutrino mixing patterns from G f = S 4 combined with H CP , broken to all possible residual symmetries of the type indicated above. Models allowing for three free parameters have been investigated in [29][30][31]. In, e.g., [30], the author has considered G f = A 5 combined with H CP , which are broken to G e = Z 2 and G ν = Z 2 × H ν CP . In this case, the matrix U e depends on an angle θ e and a phase δ e , while the matrix U ν depends on an angle θ ν . In these two scenarios the leptonic CPV phases possess non-trivial values.
The specific correlations between the values of the three neutrino mixing angles, which characterise the one-parameter models based on G e = Z n , n > 2, or Z m × Z k , m, k ≥ 2, and G ν = Z 2 × H ν CP , do not hold in the two-and three-parameter models. In addition, the Dirac CPV phase in the two-and three-parameter models is predicted to have non-trivial values which are correlated with the values of the three neutrino mixing angles and differ from 0, π, π/2 and 3π/2, although the deviations from, e.g., 3π/2 can be relatively small. The indicated differences between the predictions of the models based on G e = Z n , n > 2, or Z m × Z k , m, k ≥ 2, and on G e = Z 2 symmetries make it possible to distinguish between them experimentally by improving the precision on each of the three measured neutrino mixing angles θ 12 , θ 23 and θ 13 , and by performing a sufficiently precise measurement of the Dirac phase δ.
In the present article, we investigate the possible neutrino mixing patterns generated by a G f = S 4 symmetry combined with an H CP symmetry when these symmetries are broken down to G e = Z 2 and G ν = Z 2 × H ν CP . In section 2, we describe a general framework for deriving the form of the PMNS matrix, dictated by the chosen residual symmetries. Then, in section 3, we apply this framework to G f = S 4 combined with H CP and obtain all phenomenologically viable mixing patterns. Next, in section 4, using the obtained predictions for the neutrino mixing angles and the Dirac and Majorana CPV phases, we derive predictions for the neutrinoless double beta decay effective Majorana mass. Section 5 contains the conclusions of the present study.

The framework
We start with a non-Abelian flavour symmetry group G f , which admits a faithful irreducible 3-dimensional representation ρ. The three generations of left-handed (LH) leptons are assigned to this representation. Apart from that, the high-energy theory respects also the GCP symmetry H CP , which is implemented consistently along with the flavour symmetry.

JHEP12(2017)022
At some flavour symmetry breaking scale G f H CP gets broken down to residual symmetries G e and G ν of the charged lepton and neutrino mass terms, respectively. The residual flavour symmetries are Abelian subgroups of G f . The symmetries G e and G ν significantly constrain the form of the neutrino mixing matrix U PMNS , as we demonstrate below.

The PMNS matrix from
We choose G e to be a Z 2 symmetry. We will denote it as Z ge 2 ≡ {1, g e }, g 2 e = 1 being an element of G f of order two, generating the Z ge 2 subgroup. The invariance of the charged lepton mass term under G e implies Lets Ω e be a diagonalising unitary matrix of ρ(g e ), such that This result is obtained as follows. The diagonal entries of ρ(g e ) d are constrained to be ±1, since this matrix must still furnish a representation of Z 2 and hence its square is the identity. We have assumed that the trace of ρ(g e ) is −1, for the relevant elements g e , as it is the case for the 3-dimensional representation of S 4 we will consider later on. 6 Note that we can take the order of the eigenvalues of ρ(g e ) as given in eq. (2.3) without loss of generality, as will become clear later.

JHEP12(2017)022
The diagonal phase matrix is, however, unphysical, since it can be eliminated by rephasing of the charged lepton fields, and we will not keep it in the future. Thus, we arrive to the conclusion that the matrix U e diagonalising M e M † e reads U e = Ω e U 23 (θ e , δ e ) † P T e , (2.7) with and P e being one of six permutation matrices, which need to be taken into account, since in the approach under consideration the order of the charged lepton masses is unknown. The six permutation matrices read: Note that the order of indices in P ijk stands for the order of rows, i.e., when applied from the left to a matrix, it gives the desired order, i-j-k, of the matrix rows. The same is also true for columns, when P ijk is applied from the right, except for P 231 which leads to the 3-1-2 order of columns and P 312 yielding the 2-3-1 order.
In the neutrino sector we have a G ν = Z 2 × H ν CP residual symmetry. We will denote the Z 2 symmetry of the neutrino mass matrix as Z gν 2 ≡ {1, g ν }, with g 2 ν = 1 being an element of G f , generating the Z gν 2 subgroup. H ν CP = {X ν } is the set of remnant GCP unitary transformations X ν forming a residual CP symmetry of the neutrino mass matrix. H ν CP is contained in H CP = {X} which is the GCP symmetry of the high-energy theory consistently defined along with the flavour symmetry G f . 7 The invariance under G ν of the 7 It is worth to comment here on the notation H ν CP we use. When we write in what follows H ν CP = {Xν1, Xν2}, we mean a set of GCP transformations (Xν1 and Xν2) compatible with the residual flavour Z gν 2 symmetry (see eq. (2.13)). However, when writing Gν = Z gν 2 × H ν CP , H ν CP is intended to be a group generated by Xν1. Namely, following appendix B in [19], H ν CP is isomorphic to {I, Xν1}, where I is the unit matrix and acts again on (ϕ, ϕ * ) T . Finally, it is not difficult to convince oneself that the full residual symmetry group Gν is given by a direct product Z gν 2 × H ν CP , and the second GCP transformation Xν2 = ρ(gν ) Xν1 is contained in it. The same logic applies to the notation HCP, and, as has been shown in appendix B of [19], the full symmetry group is a semi-direct product G f HCP. Note that these notations are widely used in the literature.

JHEP12(2017)022
neutrino mass matrix implies that the following two equations hold: In addition, the consistency condition between Z gν 2 and H ν CP has to be respected: To derive the form of the unitary matrix U ν diagonalising the neutrino Majorana mass matrix M ν as m j > 0 being the neutrino masses, we will follow the method presented in [32].
Lets Ω ν1 be a diagonalising unitary matrix of ρ(g ν ), such that Expressing ρ(g ν ) from this equation and substituting it in the consistency condition, eq. (2.13), we find meaning that Ω † ν1 X ν Ω * ν1 is a block-diagonal matrix, having the form of eq. (2.5). Moreover, this matrix is symmetric, since the GCP transformations X ν have to be symmetric in order for all the three neutrino masses to be different [19,21], as is required by the data. In appendix A we provide a proof of this. Being a complex (unitary) symmetric matrix, it is diagonalised by a unitary matrix Ω ν2 via the transformation: The matrix (Ω † ν1 X ν Ω * ν1 ) d is, in general, a diagonal phase matrix. However, we can choose (Ω † ν1 X ν Ω * ν1 ) d = diag(1, 1, 1) as the phases of (Ω † ν1 X ν Ω * ν1 ) d can be moved to the matrix Ω ν2 . With this choice we obtain the Takagi factorisation of the X ν (valid for unitary symmetric matrices): with Ω ν = Ω ν1 Ω ν2 . Since, as we have noticed earlier, Ω † ν1 X ν Ω * ν1 has the form of eq. (2.5), the matrix Ω ν2 can be chosen without loss of generality to have the form of eq. (2.5) with a unitary 2 × 2 matrix in the 2-3 block. This implies that the matrix Ω ν = Ω ν1 Ω ν2 also diagonalises ρ(g ν ). Indeed, where we have used eq. (2.15). We substitute next X ν from eq. (2.18) in the GCP invariance condition of the neutrino mass matrix, eq. (2.12), and find that the matrix Ω T ν M ν Ω ν is real. Furthermore, this is a symmetric matrix, since the neutrino Majorana mass matrix M ν is symmetric. A real JHEP12(2017)022 symmetric matrix can be diagonalised by a real orthogonal transformation. Employing eqs. (2.19) and (2.11), we have (2.20) implying that Ω T ν M ν Ω ν is a block-diagonal matrix as in eq. (2.5). Thus, the required orthogonal transformation is a rotation in the 2-3 plane on an angle θ ν : (2.21) Finally, the matrix U ν diagonalising M ν reads where P ν is one of the six permutation matrices, which accounts for different order of m j , and the matrix Q ν renders them positive. Without loss of generality Q ν can be parametrised as follows: Assembling together the results for U e and U ν , eqs. (2.7) and (2.22), we obtain for the form of the PMNS matrix: Thus, in the approach we are following the PMNS matrix depends on three free real parameters 8 -the two angles θ e and θ ν and the phase δ e . One of the elements of the PMNS matrix is fixed to be a constant by the employed residual symmetries. We note finally that, since R 23 (θ ν + π) = R 23 (θ ν ) diag(1, −1, −1), where the diagonal matrix can be absorbed into Q ν , and U 23 (θ e + π, δ e ) = diag(1, −1, −1) U 23 (θ e , δ e ), where the diagonal matrix contributes to the unphysical charged lepton phases, it is sufficient to consider θ e and θ ν in the interval [0, π).

Conjugate residual symmetries
In this subsection we briefly recall why the residual symmetries G e and G ν conjugate to G e and G ν , respectively, under the same element of the flavour symmetry group G f lead to the same PMNS matrix (see, e.g., [19,22] Xν Ω * ν1 ) d = diag(1, 1, 1), and thus the matrix Ων = Ων1 Ων2 in eq. (2.18), is determined up to a multiplication by an orthogonal matrix O on the right. The matrix Ων2 O must be unitary since it diagonalises a complex symmetric matrix, which implies that O must be unitary in addition of being orthogonal, and therefore must be a real matrix.
Equation (2.19) restricts further this real orthogonal matrix O to have the form of a real rotation in the 2-3 plane, which can be "absorbed" in the R23(θ ν ) matrix in eq. (2.24).

Phenomenologically non-viable cases
Here we demonstrate that at least two types of residual symmetries {G e , G ν } = {Z ge 2 , Z gν 2 × H ν CP }, characterised by certain g e and g ν , cannot lead to phenomenologically viable form of the PMNS matrix.
• Type I: g e = g ν . In this case, we can choose Ω e = Ω ν P , with P 123 or P 132 . Then, eq. (2.24) yields U PMNS = P e U 23 (θ e , δ e ) P R 23 (θ ν ) P ν Q ν . (2.32) This means that up to permutations of the rows and columns U PMNS has the form of eq. (2.5), i.e., contains four zero entries, which are ruled out by neutrino oscillation data [26,27].

JHEP12(2017)022
Taking into account that Ω ν = Ω ν1 Ω ν2 , with Ω ν2 of the block-diagonal form given in eq. (2.5), we obtain where P T Ω ν2 , depending on P , can take one of the following forms: As a consequence, U PMNS up to permutations of the rows and columns has the form containing one zero element, which is ruled out by the data.  From 24 elements of the group there are nine elements of order two, which belong to two of five conjugacy classes of S 4 (see, e.g., [21]):

Mixing patterns from
Each of these nine elements generates a corresponding Z 2 subgroup of S 4 . Each subgroup can be the residual symmetry of M e M † e , and, combined with compatible CP transformations, yield the residual symmetry of M ν . Hence, we have 81 possible pairs of only residual flavour symmetries (taking into account remnant CP symmetries increases the number of possibilities). Many of them, however, being conjugate to each other, will lead to the same form of the PMNS matrix, as explained in subsection 2.2. Thus, we first identify the pairs of elements {g e , g ν }, which are not related by the similarity transformation given in eq. (2.25). We find nine distinct cases for which {g e , g ν } can be chosen as , S, U , SU } subgroups (see, e.g., [34]). Thus, we are left with three cases in eq. (3.5).
We have chosen g ν in such a way that it is S, U or T U for all the cases in eq. (3.5). Now we need to identify the remnant CP transformations X ν compatible with each of these three elements. It is known that the GCP symmetry H CP = {X} compatible with G f = S 4 is of the same form of G f itself [20], i.e., Thus, to find X ν compatible with g ν of interest, we need to select those X = ρ(g), which i) satisfy the consistency condition in eq. (2.13) and ii) are symmetric in order to avoid partially degenerate neutrino mass spectrum, as was noted earlier. The result reads: 9 A GCP transformation in parentheses appears automatically to be a remnant CP symmetry of M ν , if X ν which precedes this in the list is a remnant CP symmetry. This is a consequence of eqs. (2.11) and (2.12), which imply that if X ν is a residual CP symmetry of M ν , then ρ(g ν )X ν is a residual CP symmetry as well. for the neutrino mass spectrum with inverted ordering (IO). The ranges in eqs. (3.10) and (3.11) differ a little from the results obtained in [27].

Explicit forms of the PMNS matrix
First, we present an explicit example of constructing the PMNS matrix in the case of g e = S, g ν = T U and H ν CP = {U, T }, which is the first case out of the seven potentially viable cases indicated above. We will work in the basis for S 4 from [36], in which the matrices for the generators S, T and U in the 3-dimensional representation read where ω = e 2πi/3 . For simplicity we use the same notation (S, T and U ) for the generators and their 3-dimensional representation matrices. We will follow the procedure described in subsection 2.1. The matrix Ω e which diagonalises ρ(g e ) = S (see eq. (2.3)) is given by The matrix Ω ν , such that Ω ν Ω T ν = U (see eq. (2.18)), reads (3.14) Using the master formula in eq. (2.24), we obtain that up to permutations of the rows and columns U PMNS has the form  where "×" entries are functions of the free parameters θ ν , θ e and δ e . Taking into account the current data, eqs. (3.10) and (3.11), the fixed element with the absolute value of 1/ √ 2 ≈ 0.707 can be (U PMNS ) µ2 , (U PMNS ) µ3 , (U PMNS ) τ 2 or (U PMNS ) τ 3 . Note that |(U PMNS ) τ 2 | = 0.707 is outside the 3σ range in the case of the NO neutrino mass spectrum, while |(U PMNS ) µ2 | = 0.707 is at the border of the 3σ allowed ranges for both the NO and IO spectra.
Let us comment now on the following issue. Once one of the elements of the PMNS matrix is fixed to be a constant, we still have four possible configurations, namely, a permutation of two remaining columns, a permutation of two remaining rows and both of them. For instance, in the case considered above, except for P e = P ν = P 213 , we can have a fixed (U PMNS ) µ2 with (P e , P ν ) = (P 213 , P 231 ), (P 312 , P 213 ) and (P 312 , P 231 ). These combinations of the permutation matrices will not lead, however, to different mixing patterns by virtue of the following relations:
We list in table 2 the matrices Ω e and Ω ν for all seven phenomenologically viable pairs of residual symmetries {G e , G ν } = {Z ge 2 , Z gν 2 ×H ν CP }. It turns out, however, that four of these seven pairs, namely, lead to the same predictions for the mixing parameters. We demonstrate this in appendix C. Table 2. The matrices Ω e and Ω ν dictated by the residual symmetries G e = Z ge 2 and G ν = Z gν 2 ×H ν CP for all seven phenomenologically viable pairs of G e and G ν . For each pair H ν CP = {X ν1 , X ν2 } of remnant GCP transformations, the given matrix Ω ν provides the Takagi factorisation of the first element, i.e., X ν1 = Ω ν Ω T ν . 10

Extracting mixing parameters and statistical analysis
In this subsection we perform a statistical analysis of the predictions for the neutrino mixing angles and CPV phases for each of the four distinctive sets of the residual flavour and CP 10 Xν2 is instead factorised as Xν2 =ΩνΩ T ν , withΩν = Ων diag(1, i, i), as follows from This allows us to derive predictions for the three neutrino mixing angles and the three leptonic CPV phases, which, in many of the cases analysed in the present study is impossible to obtain purely analytically.
Once a pair of residual symmetries and the permutation matrices P e and P ν are specified, we have the expressions for sin 2 θ ij in terms of θ ν , θ e and δ e of the type of eqs. (3.16)- (3.18). Moreover, employing a sum rule for cos δ analogous to that in eq. (3.19) and computing the rephasing invariant which determines the magnitude of CPV effects in neutrino oscillations [37] and which in the standard parametrisation of the PMNS matrix is proportional to sin δ, we know the value of δ for any θ ν , θ e and δ e . Similarly, making use of the two charged lepton rephasing invariants, 11 associated with the Majorana phases [38][39][40][41], 25) and the corresponding real parts which in the standard parametrisation of the PMNS matrix read: I 1 = sin θ 12 cos θ 12 cos 2 θ 13 sin (α 21 /2) , I 2 = cos θ 12 sin θ 13 cos θ 13 sin (α 31 /2−δ) , (3.27) R 1 = sin θ 12 cos θ 12 cos 2 θ 13 cos (α 21 /2) , R 2 = cos θ 12 sin θ 13 cos θ 13 cos (α 31 /2−δ) , (3.28) we also obtain the values of α 21 and α 31 for any θ ν , θ e and δ e . Further, we scan randomly over θ ν ∈ [0, π), θ e ∈ [0, π) and δ e ∈ [0, 2π) and calculate the values of sin 2 θ ij and the CPV phases. We require sin 2 θ ij to lie in the corresponding 3σ ranges given in table 1. The obtained values of sin 2 θ ij and δ can be characterised by a certain value of the χ 2 function constructed as follows: where x = {x i } = (sin 2 θ 12 , sin 2 θ 13 , sin 2 θ 23 , δ) and χ 2 i are one-dimensional projections for NO and IO taken from [26]. 12 Thus, we have a list of points (sin 2 θ 12 , sin 2 θ 13 , sin 2 θ 23 , δ, 11 In their general form, when one keeps explicit the unphysical phases ξj in the Majorana condition C νj T = ξj νj, j = 1, 2, 3, the rephasing invariants related to the Majorana phases involve ξj and are invariant under phase transformations of both the charged lepton and neutrino fields (see, for example, eqs. (22)- (28) in [38]). We have set ξj = 1. 12 We note that according to the latest global oscillation data, there is an overall preference for NO over IO of ∆χ 2 IO−NO ≈ 3.6. Nevertheless, we take a conservative approach and treat both orderings on an equal footing. A discussion on this issue can be found in [26].

JHEP12(2017)022
. To see the restrictions on the mixing parameters imposed by flavour and CP symmetries we consider all 15 different pairs (a, b) of the mixing parameters. For each pair we divide the plane (a, b) into bins and find a minimum of the χ 2 function in each bin. We present results in terms of heat maps with colour representing a minimal value of χ 2 in each bin. The results obtained in each case are discussed in the following subsection.

Results and discussion
In this subsection we systematically go through all different potentially viable cases and summarise their particular features. All these cases can be divided in four groups corresponding to a particular pair of residual symmetries {G e , G ν }.
In each case we concentrate on results for the ordering for which a better compatibility with the global data is attained. Note that results for NO and IO differ only i) due to the fact that the 3σ ranges of sin 2 θ 13 and sin 2 θ 23 depend slightly on the ordering and ii) in the respective χ 2 landscapes. Moreover, we present numerical results for the Majorana phases obtained for k 1 = k 2 = 0, where k 1 and k 2 are defined in eq. (2.23). However, one should keep in mind that all four (k 1 , k 2 ) pairs, where k i = 0, 1, are allowed. Whenever k 1(2) = 1, the predicted range for α 21(31) shifts by π. The values of the k i are important for the predictions of the neutrinoless double beta decay effective Majorana mass (see, e.g., [38,[42][43][44]), which we obtain in section 4.
Using the corresponding matrices Ω e and Ω ν from table 2 and the master formula for the PMNS matrix in eq. (2.24), we find the following form of the PMNS matrix (up to permutations of rows and columns and the phases in the matrix Q ν ): with c e ≡ cos θ e , s e ≡ sin θ e , c ν ≡ cos θ ν , s ν ≡ sin θ ν and From eq. (3.30), we see that the absolute values of the elements of the first row are fixed. Namely, the modulus of the first element is equal to 1/ √ 2, while the moduli of the second and third elements equal 1/2. Taking into account the current knowledge of the mixing parameters, eqs. (3.10) and (3.11), this implies that there are only two potentially viable In this case we obtain This means that only a narrow interval sin 2 θ 23 ∈ [0.510, 0.512] is allowed using the 3σ region for sin 2 θ 13 . From the equality |(U PMNS ) µ1 | = 1/2, which we find to hold in this case, it follows that cos δ satisfies the following sum rule: cos δ = 1 − 4 sin 2 θ 12 cos 2 θ 23 − 4 cos 2 θ 12 sin 2 θ 23 sin 2 θ 13 2 sin 2θ 12 sin 2θ 23 sin θ 13 , where the mixing angles in addition are correlated among themselves. We find that sin 2 θ 13 is constrained to lie in the interval (0.0213, 0.0240 (2)] for NO (IO) and, hence, sin 2 θ 23 in [0.5109, 0.5123 (4)]. This range of values of sin 2 θ 23 is not compatible with its current 2σ range. Moreover, sin 2 θ 12 is found to be between approximately 0.345 and 0.354, which is outside its current 2σ range as well. What concerns the CPV phases, the predicted values of δ are distributed around 0, namely, δ ∈ [−0.11π, 0.11π], of α 21 around π, α 21 ∈ (0.93π, 1.07π), while the values of α 31 fill the whole range, i.e., α 31 ∈ [0, 2π). These numbers, presented for the NO spectrum, remain practically unchanged for the IO spectrum. However, the global minimum χ 2 min of the χ 2 function, defined in eq. (3.29), yields approximately 22 (19) for NO (IO), which implies that this case is disfavoured by the global data at more than 4σ.
It is worth noting that we should always keep in mind the correlations between the mixing angles in expressions of this type. The values of δ in this case lie around π, in the interval [0.89π, 1.11π]. As in the previous case, the global minimum of χ 2 is somewhat large, χ 2 min ≈ 18.5 (15) for NO (IO), meaning that this case is also disfavoured.

JHEP12(2017)022
For this choice of the residual symmetries, the PMNS matrix reads (up to permutations of rows and columns and the phases in the matrix Q ν ): Equation (3.41) implies that the absolute value of one element of the PMNS matrix is predicted to be 1/ √ 2. Thus, we have four potentially viable cases.
Note that from eqs. (3.10) and (3.11) it follows that this magnitude of the fixed element is inside its 3σ range for NO, but slightly outside the corresponding range for IO. Hence, we will focus on the results for NO. The characteristic feature of this case is the following sum rule for cos δ: cos δ = 2 cos 2 θ 12 cos 2 θ 23 + 2 sin 2 θ 12 sin 2 θ 23 sin 2 θ 13 − 1 sin 2θ 12 sin 2θ 23 sin θ 13 , (3. 46) which arises from the equality of |(U PMNS ) µ2 | to 1/ √ 2. The pair correlations between the mixing parameters in this case are summarised in figure 1. The colour palette corresponds to values of χ 2 for NO. As can be seen, while all values of sin 2 θ 13 in its 3σ range are allowed, the parameters sin 2 θ 12 and sin 2 θ 23 are found to lie in [0.250, 0.308] and [0.381, 0.425) intervals, respectively. The predicted values of δ span the range [0.68π, 1.32π]. Thus, CPV effects in neutrino oscillations due to the phase δ can be suppressed. The Majorana phases instead are distributed in relatively narrow regions around 0, so the magnitude of the neutrinoless double beta decay effective Majorana mass (see section 4 and, e.g., [38,[42][43][44]) is predicted (for k 1 = k 2 = 0) to have a value close to the maximal possible for the NO spectrum. Namely, α 21 ∈ [−0.16π, 0.16π] and α 31 ∈ (−0.13π, 0.13π). In addition, δ is strongly correlated with α 21 and α 31 , which in turn exhibit a strong correlation between themselves. Finally, χ 2 min ≈ 7 for both NO and IO, i.e., this case is compatible with the global data at less than 3σ. 13

JHEP12(2017)022
• Case B2: |(U PMNS ) τ 2 | = 1/ √ 2 (P e = P 321 , P ν = P 213 ). Note that this value of |(U PMNS ) τ 2 | is compatible at 3σ with the global data in the case of IO spectrum, but not in the case of NO spectrum, as can be seen from eqs. (3.10) and (3.11). Thus, below we present results for the IO spectrum only. As in case B1, the whole 3σ range for sin 2 θ 13 is allowed. The obtained ranges of values of α 21 and α 31 are the same of the preceding case. The range for sin 2 θ 12 differs somewhat from that obtained in case B1, and it reads sin 2 θ 12 ∈ [0.250, 0.328]. 14 The predictions for sin 2 θ 23 and δ are different. Now the following sum rule, derived from |(U PMNS ) τ 2 | = 1/ √ 2, holds: cos δ = 1 − 2 cos 2 θ 12 sin 2 θ 23 − 2 sin 2 θ 12 cos 2 θ 23 sin 2 θ 13 sin 2θ 12 sin 2θ 23 sin θ 13 . (3.47) The values of δ are concentrated in [−0.38π, 0.38π]. For sin 2 θ 23 we find the range (0.575, 0.636]. The correlations between the phases are of the same type as in case B1. We summarise the results in figure 2. Finally, χ 2 min ≈ 6 in the case of IO and χ 2 min ≈ 12.5 for NO, which reflects incompatibility of this case at more than 3σ for the NO spectrum. This occurs mainly due to the predicted values of sin 2 θ 23 , which are outside its current 2σ range for NO. the angles θ 13 and θ 23 are correlated as in case A1, i.e., according to eq. (3.35). For IO this leads to sin 2 θ 23 ∈ [0.5097, 0.5124] due to the fact that the whole 3σ range of sin 2 θ 13 is found to be allowed, as can be seen from figure 3. Note that this range is outside the current 2σ range of sin 2 θ 23 . In addition, we find that the whole 3σ range of the values of sin 2 θ 12 can be reproduced. In contrast to case A1, |(U PMNS ) µ1 | does not equal 1/2, but depends on θ ν in the following way: From this equation we find cos δ = 1 − 4 sin 2 θ 12 cos 2 θ 23 − 4 cos 2 θ 12 sin 2 θ 23 sin 2 θ 13 − sin 2θ ν 2 sin 2θ 12 sin 2θ 23 sin θ 13 , (3.49) i.e., cos δ depends on θ ν explicitly (not only via θ 12 , θ 23 and θ 13 ). With this relation, any value of δ between 0 and 2π is allowed (see figure 3). The Majorana phases, however, are constrained to lie around 0 in the following intervals: α 21 ∈ [−0.23π, 0.23π] and α 31 ∈ (−0.18π, 0.18π). Moreover, both phases α 21 and α 31 are correlated in one and the same peculiar way with the phase δ. The correlation between α 21 and α 31 is similar to those in cases B1 and B2 (cf. figures 1 and 2). Due to the predicted values of sin 2 θ 23 , which belong to the upper octant, IO is preferred over NO, the corresponding χ 2 min being approximately 5 and 8.5.

JHEP12(2017)022
• Case B4: |(U PMNS ) τ 3 | = 1/ √ 2 (P e = P ν = P 321 ). The predicted ranges of all the mixing parameters are the same of case B3, except for sin 2 θ 23 , which respects the relation in eq. (3.38), and thus belongs to [0.4876, 0.4903] in the case of IO spectrum. As in the previous case, this interval falls outside the 2σ range of sin 2 θ 23 . The results obtained in this case for the IO spectrum are presented in figure 4. Similarly to the preceding case, we find which leads to cos δ = sin 2θ ν + 4 sin 2 θ 12 sin 2 θ 23 + 4 cos 2 θ 12 cos 2 θ 23 sin 2 θ 13 − 1 2 sin 2θ 12 sin 2θ 23 sin θ 13 . (3.51) The correlation between the Majorana phases is similar to that in the previous case. Also in this case, χ 2 min ≈ 4.5 for IO is lower than that of approximately 6.5 for NO, the reason being again the predicted range of sin 2 θ 23 .

JHEP12(2017)022
i.e., cos δ explicitly depends on θ e , and eventually this relation does not constrain δ. Instead the Majorana phase α 21 is predicted to be exactly π (exactly 0) for k 1 = 0 (k 1 = 1). While the second Majorana phase α 31 itself remains unconstrained, the difference α 31 −2δ = 0 (π) for k 2 = 0 (k 2 = 1), i.e., we have a strong linear correlation between δ and α 31 (see figure 5). The reason for these trivial values of α 21 and α 31 −2δ is the following. In the standard parametrisation of the PMNS matrix, α 21 and the combination (α 31 − 2δ) may be extracted from the phases of the first row of the PMNS matrix, as can be seen from eqs. (3.25)-(3.28). In case C1, none of the phases of the first row elements of the PMNS matrix depend (mod π) on the free parameters θ ν , θ e and δ e . Namely, the phases of (U PMNS ) e1 , (U PMNS ) e2 and (U PMNS ) e3 are fixed (mod π and up to a global phase) to be −π/6, π/3 and −π/6, respectively. Notice that only in groups B and C the relative phases of the first row can be predicted (mod π) to be independent of θ ν , θ e and δ e . Furthermore, case C1 stands out since it is, out of these relevant cases, the only one which survives the constraints on the magnitudes of the PMNS matrix elements given in eqs. (3.10) and (3.11). Finally, χ 2 min ≈ 7 for both mass orderings.
• Case C2: |(U PMNS ) µ1 | = 1/2 (P e = P 213 , P ν = P 123 ). The correlations between the mixing parameters obtained in this case for NO are summarised in figure 6 (the results for IO are very similar). This case accounts for the whole 3σ range of sin 2 θ 13 , but constrains the values of the two other angles. Namely, we find sin 2 θ 12 ∈ [0.285, 0.354] and sin 2 θ 23 ∈ [0.381, 0.524]. This case enjoys the sum rule for cos δ given in eq. (3.37), since |(U PMNS ) µ1 | = 1/2 as it was also in case A1. As a consequence, we find δ to be constrained: δ ∈ (−0.38π, 0.38π). Both Majorana phases are distributed in relatively narrow intervals around π: α 21 ∈ (0.85π, 1.15π) and α 31 ∈ [0.91π, 1.09π]. The phase δ is correlated with each of the two Majorana phases in a similar way. The latter in turn are correlated linearly between themselves. Overall, NO is slightly preferred over IO in this case. The corresponding values of χ 2 min read 4.5 and 5.5, respectively.
• Case C3: |(U PMNS ) τ 1 | = 1/2 (P e = P 321 , P ν = P 123 ). This case shares some of the predictions of case C2. Namely, the whole 3σ range of sin 2 θ 13 is allowed, and the ranges of α 21 and α 31 are the same as in the preceding case, as can be seen from figure 7, in which we present the results for the IO neutrino mass spectrum.  6 and 7). Due to the predicted range of sin 2 θ 23 , this case is favoured by the data for IO, for which χ 2 min ≈ 1.5, while for NO we find χ 2 min ≈ 8.5.

JHEP12(2017)022
• Case C4: |(U PMNS ) µ2 | = 1/2 (P e = P ν = P 213 ). From eqs. (3.10) and (3.11) it follows that the value of |(U PMNS ) µ2 | = 1/2 is allowed at 3σ only for IO. Thus, below we present results obtained in the IO case. In the case under consideration there are no constraints on the ranges of sin 2 θ 12 and sin 2 θ 13 . The atmospheric angle is, in turn, found to lie in the upper octant, sin 2 θ 23 ∈ (0.505, 0.636]. As can be seen in figure 8, δ ∈ [−0.54π, 0.54π], which is a consequence of the following correlation between cos δ and the mixing angles: cos δ = 4 cos 2 θ 12 cos 2 θ 23 + 4 sin 2 θ 12 sin 2 θ 23 sin 2 θ 13 − 1 2 sin 2θ 12 sin 2θ 23 sin θ 13 , The allowed values of δ span the range [0.51π, 1.49π]. The correlations between the mixing parameters in this case are summarised in figure 9. Finally, we have χ 2 min ≈ 0.5 for both NO and IO. For this last group of cases, we find that the PMNS matrix takes the following form (up to permutations of rows and columns and the phases in the matrix Q ν ):

JHEP12(2017)022
Therefore, the absolute value of the fixed element of the neutrino mixing matrix yields 1/2. Thus, we have again five potentially viable cases.
• Case D1: |(U PMNS ) e2 | = 1/2 (P e = P 123 , P ν = P 213 ). In this case we find  (31) and all the other mixing parameters in this case, as can be seen in figure 10, in which we summarise the results for NO. Finally, χ 2 min ≈ 4.5 for NO, and it is slightly higher, χ 2 min ≈ 5.5, for IO.
• Case D3: |(U PMNS ) τ 1 | = 1/2 (P e = P 321 , P ν = P 123 ). As in the previous case, the whole 3σ range of sin 2 θ 13 gets reproduced. The allowed ranges of sin 2 θ 12 , α 21 and , which holds in this case since |(U PMNS ) τ 1 | = 1/2, leads to the values of δ distributed around π in a rather broad range of (0.59π, 1.41π). The correlations between the Majorana phases and δ are as in the previous case, but again with an approximate shift of δ by π (see figure 11). The minimal value χ 2 min ≈ 1.5 in the IO case, while for the NO spectrum we get approximately 8.5. This difference is due to the allowed values of sin 2 θ 23 .
• Case D4: |(U PMNS ) µ2 | = 1/2 (P e = P ν = P 213 ). This case can account only for a part of the 3σ range of sin 2 θ 13 , namely, sin 2 θ 13 ∈ [0.0214, 0.0240(2)] for NO (IO) spectrum. The constraints on two other angles are more severe. We find that only a narrow region of the values of sin 2 θ 23 , which falls outside its present 2σ range, is allowed, namely, sin 2 θ 23 ∈ [0.505, 0.512]. For the solar mixing angle we have sin 2 θ 12 ∈ [0.345, 0.354], which is also outside the current 2σ range of this parameter. The sum rule in eq. (3.57), which is also valid in this case, constrains δ to lie in a narrow interval around 0: δ ∈ [−0.11π, 0.11π]. The Majorana phases instead are distributed in narrow intervals around π. Namely, α 21 ∈ (0.83π, 1.17π) and α 31 ∈ [0.92π, 1.08π]. However, the global minimum of χ 2 is somewhat large in this case for both NO and IO orderings. Namely, we find χ 2 min ≈ 22 (19) for NO (IO), i.e., this case is disfavoured at more than 4σ by the current global data.

JHEP12(2017)022
• Case D5: |(U PMNS ) τ 2 | = 1/2 (P e = P 321 , P ν = P 213 ). This last case shares the predicted ranges for sin 2 θ 12 , sin 2 θ 13 , α 21 and α 31 with case D4. Therefore, this case is also not compatible with the 2σ range of the values of sin 2 θ 12 . For sin 2 θ 23 instead we find the narrow interval in the lower octant, sin 2 θ 23 ∈ [0.488, 0.495], which lies outside the 2σ range of sin 2 θ 23 . We find cos δ to satisfy the sum rule in eq. (3.58), which in this case gives us the values of δ in a narrow interval around π, δ ∈ [0.89π, 1.11π]. Thus, all the three CPV phases are concentrated in narrow ranges around π. Finally, we find χ 2 min ≈ 18.5 (15) for NO (IO), which implies that this case is also disfavoured by the latest global neutrino oscillation data.
The PMNS matrix in case A2 is related with that in case A1 by the permutation matrix P 312 as U A2 PMNS = P 312 U A1 PMNS . Given that P 312 = P 132 P 321 , one can see that these matrices are related by µ − τ interchange, after an unphysical exchange of the first and third rows of U A1 PMNS has been performed (which amounts to a redefinition of the free parameter θ e , as shown in eq. (3.21)). The same also holds for the following pairs of cases: (B1, B2), (B3, B4), (C2, C3), (C4, C5), (D2, D3) and (D4, D5). As can be seen from the discussion above and figures 1-4 and 6-11, cases inside a pair share some qualitative features. Namely, i) the predicted ranges of sin 2 θ 12 , sin 2 θ 13 , α 21 and α 31 are approximately the same; ii) the predicted range of sin 2 θ 23 gets approximately reflected around 1/2, i.e., sin 2 θ 23 → 1 − sin 2 θ 23 ; iii) the predicted range of the CPV phase δ experiences an approximate shift by π, i.e., δ → δ + π.
In tables 3 and 4 we summarise the predicted ranges of the mixing parameters obtained in all the phenomenologically viable cases discussed above. The corresponding best fit values together with χ 2 min are presented in tables 5 and 6. Finally, in table 7 we show whether the cases compatible with the 3σ ranges of the three mixing angles are also compatible with their corresponding 2σ ranges.
The results shown in tables 3-6 allow to assess the possibilities to critically test the predictions of the viable cases of the model and to distinguish between them. We recall that the current 1σ uncertainties on the measured values of sin 2 θ 12 , sin 2 θ 13 and sin 2 θ 23 are [26] 5.8%, 4.0% and 9.6%, respectively. These uncertainties are foreseen to be further reduced by the currently active and/or future planned experiments. The Daya Bay collaboration plans to determine sin 2 θ 13 with 1σ uncertainty of 3% [45]. The uncertainties on sin 2 θ 12 and sin 2 θ 23 are planned to be reduced significantly. The parameter sin 2 θ 12 is foreseen to be measured with 1σ relative error of 0.7% in the JUNO experiment [46,47]. In the proposed upgrading of the currently taking data T2K experiment [48,49], for example, θ 23 is estimated to be determined with a 1σ error of 1.7 • , 0.5 • and 0.7 • if the best fit value of sin 2 θ 23 = 0.50, 0.43 and 0.60, respectively. This implies that for these three values of sin 2 θ 23 the absolute 1σ error would be 0.0297, 0.0086 and 0.0120. This error on sin 2 θ 23 will be further reduced in the future planned T2HK [50] and DUNE [51][52][53] experiments. If δ = 3π/2, the CP-conserving case of sin δ = 0 would be disfavoured for the NO mass spectrum in the same experiment at least at 3σ C.L. Higher precision measurements of δ are planned to be performed in the T2HK and DUNE experiments.  Figure 11. Correlations between the neutrino mixing parameters in case D3. The values of all the three mixing angles are required to lie in their respective 3σ ranges. Colour represents values of χ 2 for the IO neutrino mass spectrum.

JHEP12(2017)022
We turn now to the possibilities to discriminate experimentally between the different cases listed in tables 3-6 using the prospective data on sin 2 θ 12 , sin 2 θ 13 , sin 2 θ 23 and δ. The first thing to notice is that the predicted ranges for sin 2 θ 12 , sin 2 θ 13 , sin 2 θ 23 and δ in cases A1 and A2 practically coincide with the predictions respectively in cases D4 and D5. However, cases A1, D4 and cases A2, D5 are strongly disfavoured by the current data: for the NO (IO) neutrino mass spectrum A1 and D4 are disfavoured at 4.7σ (4.4σ), while A2 and D5 are disfavoured at 4.3σ (3.9σ). In all these cases sin 2 θ 12 , in particular, is predicted to lie in the interval (0.345,0.354) compatible with the current 3σ range and, given the current best fit value of sin 2 θ 12 and prospective JUNO precision on sin 2 θ 12 , it is very probable that future more precise data on sin 2 θ 12 will rule out completely these scenarios. We will not discuss them further in this subsection.
It follows also from tables 5 and 6 that the combined results on the best fit values of sin 2 θ 12 , sin 2 θ 23 and δ we have obtained in the different viable cases (excluding A1, A2, D4 and D5) differ significantly. Assuming, for example, that the experimentally determined best fit values of sin 2 θ 12 and sin 2 θ 23 will coincide with those found by us for a given viable case, it is not difficult to convince oneself inspecting tables 5 and 6 that the cited prospective 1σ errors on sin 2 θ 12 and sin 2 θ 23 will allow to discriminate between the different viable cases identified in our study. More specifically, considering as an example only the case of NO neutrino mass spectrum, the prospective high precision measurement of sin 2 θ 12 will allow to discriminate between case C1 and all other cases B1-B4, C2-C5, D2 and D3. The same measurement will make it possible to distinguish i) between case B1 and all the other cases except B2, ii) between case B2 and all the other cases except B1, B3 and B4, and similarly iii) between case B3 and all the other cases except B2, B4, C4 and C5. However, the differences between the best fit values of sin 2 θ 23 in cases B1, B2 and B3 (or B4) are sufficiently large, which would permit to distinguish between these three cases if sin 2 θ 23 were measured with the prospective precision. It follows from table 5, however, that it would be very challenging to discriminate between cases B3 and B4: it will require extremely high precision measurement of sin 2 θ 23 . These two cases would be ruled out, however, if the experimentally determined best fit value of sin 2 θ 23 differs significantly from the results for sin 2 θ 23 , namely, 0.511 and 0.489, we have obtained for sin 2 θ 23 in the B3 and B4 cases.
In the remaining cases C2-C5 and D2-D3, the results we have obtained for sin 2 θ 12 , as table 6 shows, are very similar. However, the predictions for the pair sin 2 θ 23 and δ differ significantly in cases C2 or D2, and C3 or D3. The cases within each pair would be ruled out if the experimentally determined values of sin 2 θ 23 and δ differ significantly from the predicted best fit values.
Thus, the planned future high precision measurements of sin 2 θ 12 and sin 2 θ 23 , together with more precise data on the Dirac phase δ, will make it possible to critically test the predictions of the cases listed in tables 3-6. A comprehensive analysis of the possibilities to distinguish between the different viable cases found in our work in the considered S 4 model can only be done when more precise data first of all on sin 2 θ 12 and sin 2 θ 23 , and then on δ, will be available. We schematically summarise in figure 12 the predicted 3σ allowed regions in the plane (sin 2 θ 23 , sin 2 θ 12 ) for all currently viable cases from figures 1-11. In this figure we also present the best fit point in each case used in the preceding discussion. When future more precise data on sin 2 θ 23 and sin 2 θ 12 become available, the experimentally allowed region in the (sin 2 θ 23 , sin 2 θ 12 ) plane will shrink, and only a limited number of cases, if any, will remain viable. It will be possible to distinguish further between some or all of the remaining viable cases with a high precision measurement of δ.
Finally, we note that the sum rules for sin 2 θ 23 (sin 2 θ 12 in case C1) and/or cos δ obtained in the present study follow from those derived in [7] for certain values of the parameters sin 2 θ • ij , fixed by G f = S 4 and the residual Z ge 2 and Z gν 2 flavour symmetries, and the additional constraints provided by the GCP symmetry H ν CP . Note that in [7] only flavour symmetry, without imposing a GCP symmetry, has been considered. As we have seen in subsection 2.1, a GCP symmetry does not allow for a free phase δ ν coming from the neutrino sector, which is present otherwise. This, in turn, leads to the fact that in certain cases the parameter sinθ ν ij (see eq. (213) in [7]), which is free in [7], gets fixed by the GCP symmetry. Thus, we find additional correlations between θ ij and between θ ij and cos δ in these cases. We provide the correspondence between the phenomenologically viable cases of the present study and the cases considered in [7] in appendix D.  Table 3. Ranges of the mixing parameters for the viable cases, i.e., those cases for which the predicted values of all the three mixing angles lie inside their respective 3σ allowed ranges. The cases presented here correspond to G e = Z ge 2 and G ν = Z gν 2 × H ν CP with {g e , g ν } = {T U, S}, for which the magnitude of the fixed element is 1/ √ 2 (p.f.e. denotes its position in U PMNS ). For each case, the upper and lower rows refer to NO and IO, respectively.   Table 5. Best fit values of the mixing parameters and the corresponding value of the χ 2 function, χ 2 min , for the viable cases, i.e., those cases for which the predicted values of all the three mixing angles lie inside their respective 3σ allowed ranges. The cases presented here correspond to G e = Z ge 2 and G ν = Z gν 2 × H ν CP with {g e , g ν } = {T U, S}, for which the magnitude of the fixed element is 1/ √ 2 (p.f.e. denotes its position in U PMNS ). For each case, the upper and lower rows refer to NO and IO, respectively.

Neutrinoless double beta decay
As we have seen, in the class of models investigated in the present article the Dirac and Majorana CPV phases, δ and α 21 , α 31 , are (statistically) predicted to lie in specific, in most cases relatively narrow, intervals and their values are strongly correlated. The only exception is case C1, in which the exact predictions α 21 = 0 or π and (α 31 − 2δ) = 0 or π hold.
These results make it possible to derive predictions for the absolute value of the neutrinoless double beta ((ββ) 0ν -) decay effective Majorana mass, m (see, e.g., refs. [1,[42][43][44]), as a function of the lightest neutrino mass. As is well known, information about | m | is provided by the experiments on (ββ) 0ν -decay of even-even nuclei 48 Ca, 76  can trigger the process of (ββ) 0ν -decay. In this case the (ββ) 0ν -decay amplitude has the following general form (see, e.g., refs. [42][43][44] Table 7. Compatibility of the cases under consideration with the 3σ and 2σ experimentally allowed ranges of the three neutrino mixing angles for both types of the neutrino mass spectrum. limits on | m | have been obtained by the KamLAND-Zen [54] and GERDA Phase II [55] experiments searching for (ββ) 0ν -decay of 136 Xe and 76 Ge, respectively: | m | < (0.061 − 0.165) eV [54] and | m | < (0.15 − 0.33) eV [55] , on the 136 Xe and 76 Ge (ββ) 0ν -decay half-lives (for a review of the limits on | m | obtained in other (ββ) 0ν -decay experiments and a detailed discussion of the NME calculations for (ββ) 0ν -decay and their uncertainties see, e.g., [56]). It is important to note that a large number of experiments of a new generation aims at a sensitivity to | m | ∼ (0.01−0.05) eV, which will allow to probe the whole range of the predictions for | m | in the case of IO neutrino mass spectrum [57,58] (see, e.g., [56,59] for reviews of the currently running and future planned (ββ) 0ν -decay experiments and their prospective sensitivities).

JHEP12(2017)022
In what follows, we obtain predictions for | m | using the phenomenologically viable neutrino mixing patterns found in subsection 3.4. In figures 13-16 we present | m | as a function of the lightest neutrino mass m min (m min = m 1 for the NO spectrum and m min = m 3 for the IO spectrum) in cases B1-B4, C1-C3, C4 and C5, and D2 and D3. The solid and dashed lines limit the found allowed regions of | m | calculated using the predicted ranges for θ 12 , θ 13 , α 21 , (α 31 − 2δ). In the left panels we require the predicted values of sin 2 θ 12 , sin 2 θ 13 and sin 2 θ 23 to lie in their corresponding experimentally allowed 3σ intervals, while in the right panels we require them to be inside the corresponding 2σ ranges. The mass squared differences ∆m 2 21 and ∆m 2 31(23) in the case of NO (IO) spectrum are varied in their appropriate ranges given in table 1. The light-blue (light-red) areas in the left and right panels are obtained varying the neutrino oscillation parameters θ 12 , θ 13 , ∆m 2 21 and ∆m 2 31(23) in their full 3σ and 2σ NO (IO) ranges, respectively, and varying the phases α 21 and (α 31 − 2δ) in the interval [0, 2π). The horizontal brown and grey bands indicate the current most stringent upper limits on | m |, given in eq. (4.1), set by KamLAND-Zen and GERDA Phase II, respectively. The vertical grey line represents the prospective upper limit on m min ∼ < 0.2 eV from the KATRIN experiment [61].
Several comments are in order. Firstly, for given values of (k 1 , k 2 ) and a given ordering we find | m | to be inside of a band, which occupies a certain part of the allowed parameter space. Secondly, we note that most cases are compatible with both 3σ and 2σ ranges of all the mixing angles for both neutrino mass orderings (see table 7). There are several exceptions. Namely, cases B2, C3, C4 and D3, in which, due to the correlations imposed by the employed symmetry, the predictions for sin 2 θ 23 for the NO spectrum are not compatible with its 2σ allowed range (see tables 3 and 4). Moreover, there is incompatibility for both orderings of cases B3 and B4 with the allowed 2σ ranges of sin 2 θ 23 (see table 3), and of case C1 with the 2σ range of sin 2 θ 12 (see table 4). Thirdly, the predictions for | m | compatible with the 3σ ranges of all the mixing angles are almost the same for the following pairs of cases: (B1, B2), (B3, B4), (C2, C3), (C4, C5) and (D2, D3). As discussed at the end of subsection 3.4, the cases in each pair share some qualitative features, in particular, the allowed ranges of θ 12 , θ 13 , α 21 and (α 31 − 2δ) are approximately equal. We note also that case C1 stands out by having relatively narrow bands for | m | due to the predicted values of α 21 = k 1 π and (α 31 − 2δ) = k 2 π. Finally, the results shown in figures 13-16 and derived using the predictions for the CPV phases and the mixing angles θ 12 and θ 13 in the case when the predicted values of all the three mixing angles θ 12 , θ 13 and θ 23 are compatible with their respective 3σ experimentally allowed ranges, can be obtained analytically in the limiting cases of NH, IH and QD spectra using eqs.  (1, 1) Figure 13. The magnitude of the effective Majorana mass | m | versus the lightest neutrino mass m min . The lines limit the allowed regions of | m | calculated using the predictions for the relevant mixing angles and the CPV phases obtained in cases B1-B4 and compatible with the 3σ (left panels) and 2σ (right panels) ranges of all the three mixing angles. The light-blue (light-red ) areas are obtained varying the neutrino oscillation parameters θ 12 , θ 13 , ∆m 2 21 and ∆m 2 31(23) for NO (IO) in their allowed 3σ and 2σ ranges in the left and right panels, respectively, and the phases α 21 and (α 31 − 2δ) in the interval [0, 2π). The horizontal brown and grey bands indicate the current upper bounds on | m | quoted in eq. (4.1) set by KamLAND-Zen [54] and GERDA Phase II [55], respectively. The vertical grey line represents the prospective upper limit on m min ∼ < 0.2 eV from KATRIN [61]. Cases B3 and B4 are compatible with the 3σ ranges of the mixing angles, but not with their 2σ ranges.   residual Z ge 2 and Z gν 2 × H ν CP symmetries in the charged lepton and neutrino sectors, respectively, where Z ge 2 = {1, g e }, Z gν 2 = {1, g ν } and H ν CP = {X ν }, 1 being the unit element of S 4 . The massive neutrinos are assumed to be Majorana particles with their masses generated by the neutrino Majorana mass term of the left-handed (LH) flavour neutrino fields ν lL (x), l = e, µ, τ . We show that in this class of models the three neutrino mixing angles, θ 12 , θ 23 and θ 13 , the Dirac and the two Majorana CP violation (CPV) phases, δ and α 21 , α 31 , are functions of altogether three parameters -two mixing angles and a phase, θ e , θ ν and δ e .
The S 4 group has 9 different Z 2 subgroups. Assuming that the LH flavour neutrino and charged lepton fields, ν lL (x) and l L (x), l = e, µ, τ , transform under a triplet irreducible unitary representation of S 4 , we prove that there are only 3 pairs of subgroups Z ge 2 and Z gν 2 which can lead to different viable (i.e., compatible with the current data) predictions for the lepton mixing. For these three pairs, |〈m 〉|

[eV]
Case D3; 2σ Figure 16. The same as in figure 13, but for cases D2 and D3.
sentation of S 4 (eq. (3.12)). In what concerns the residual GCP symmetry in the neutrino sector, H ν CP = {X ν }, we show that the constraints on X ν (following from the conditions of consistency between Z gν 2 and H ν CP and of having non-degenerate neutrino mass spectrum, X ν = X T ν ) are satisfied in the following cases: them case by case and have classified all phenomenologically viable mixing patterns they lead to. In all four groups of cases the PMNS neutrino mixing matrix is predicted to contain one constant element which does not depend on the three basic parameters, θ e , θ ν and δ e . The magnitude of this element is equal to 1/ √ 2 in the "Group A" cases , S}, and in the "Group B" cases of In the approach to the neutrino mixing based on S 4 flavour and GCP symmetries employed by us, the PMNS matrix is determined up to permutations of columns and rows. This implies that theoretically any of the elements of the PMNS matrix can be equal by absolute value to 1/ √ 2 in the Group A and Group B cases, and to 1/2 in the Group C and Group D cases. However, the data on the neutrino mixing angles and the Dirac phase δ imply that, taking into account the currently allowed 3σ ranges of the PMNS matrix elements (see eqs. (3.10) and (3.11)), only 4 elements, namely, (U PMNS ) µ2 , (U PMNS ) µ3 , (U PMNS ) τ 2 or (U PMNS ) τ 3 , can have an absolute value equal to 1/ √ 2 ≈ 0.707, and only 5 elements, namely, (U PMNS ) e2 , (U PMNS ) µ1 , (U PMNS ) τ 1 , (U PMNS ) µ2 or (U PMNS ) τ 2 , can have an absolute value equal to 1/2. It should be added that i) |(U PMNS ) τ 2 | = 0.707 lies outside the respective currently allowed 3σ range in the case of NO neutrino mass spectrum, ii) |(U PMNS ) µ2 | = 0.707 is slightly outside the 3σ allowed range for the IO spectrum, and that iii) the value of |(U PMNS ) µ2 | = 1/2 is allowed at 3σ only for the IO spectrum.
We have derived predictions for the six parameters of the PMNS matrix, θ 12 , θ 23 and θ 13 , δ, α 21 and α 31 , in the potentially viable cases of Groups A-D. This was done for both NO and IO neutrino mass spectra in the cases compatible at 3σ with the existing data. We have performed also a statistical analysis of the predictions for the neutrino mixing angles and CPV phases for each of these cases. We have found that in certain cases the predicted values of the neutrino mixing angles are ruled out, or are strongly disfavoured, by the existing data (see subsection 3.4 for details). These are: i) in Group A, the cases of |(U PMNS ) µ3 | = 1/ √ 2 (strongly disfavoured), and |(U PMNS ) τ 3 | = 1/ √ 2 (strongly disfavoured); ii) in Group D, the cases of |(U PMNS ) e2 | = 1/2 (ruled out), |(U PMNS ) µ2 | = 1/2 (strongly disfavoured), and |(U PMNS ) τ 2 | = 1/2 (strongly disfavoured).
The results of the statistical analysis in the viable cases are presented graphically in figures 1-11. The predicted ranges of the neutrino mixing parameters and the their corresponding best fit values are summarised in tables 3-6.
Given the difference in the currently allowed 2σ ranges of sin 2 θ 23 (see table 1), the prediction for the allowed values of sin 2 θ 23 in certain phenomenologically viable cases makes the IO (NO) spectrum statistically somewhat more favourable than the NO (IO) spectrum. At the same time, we have found that in a large number of viable cases the results we have obtained for the NO and IO spectra are very similar.
As a consequence of the fact that, in the class of models we consider, the six PMNS matrix parameters, θ 12 , θ 23 , θ 13 , δ, α 21 and α 31 , are fitted with the three basic parameters, JHEP12(2017)022 θ e , θ ν and δ e , it is not surprising that we have found that there are strong correlations i) between the values of the Dirac phase δ and the values of the two Majorana phases α 21 and α 31 , which in turn are correlated between themselves (figures 1, 2, 6-9), and depending on the case ii) either between the values of θ 12 and θ 13 (figure 5), or between the values of θ 23 and θ 13 (figures 3 and 4) or else between the values of θ 12 and θ 23 (figures 1, 2, 6-11). In certain cases our results showed strong correlations between the predicted values of θ 23 and the Dirac phase δ and/or the Majorana phases α 21,31 (figures 8-11).
In the cases of i) Group B with |(U PMNS ) µ2 | = 1/ √ 2, or |(U PMNS ) τ 2 | = 1/ √ 2, ii) Group C with |(U PMNS ) µ1 | = 1/2, or |(U PMNS ) τ 1 | = 1/2, or |(U PMNS ) µ2 | = 1/2, or |(U PMNS ) τ 2 | = 1/2, and iii) Group D with |(U PMNS ) µ1 | = 1/2, or |(U PMNS ) τ 1 | = 1/2, the cosine of the Dirac phase δ satisfies a sum rule by which it is expressed in terms of the three neutrino mixing angles θ 12 , θ 23 and θ 13 . Taking into account the ranges and correlations of the predicted values of the three neutrino mixing angles, δ is predicted to lie in certain, in most of the discussed cases rather narrow, intervals (subsection 3.4). In the remaining viable cases of Groups B and C, cos δ was shown to satisfy sum rules which depend explicitly, in addition to θ 12 , θ 23 and θ 13 , on one of the three basic parameters of the class of models considered, θ e or θ ν . In these cases, as we have shown, cos δ can take any value.
We have derived also predictions for the Majorana CPV phases α 21 and α 31 in all viable cases of Groups B, C and D (subsection 3.4). With one exception -the case of |(U PMNS ) e2 | = 1/2 of Group C -the values of α 21 and α 31 , as we have indicated earlier, are strongly correlated between themselves. In case C1 there is a strong linear correlation between α 31 and δ.
Using the predictions for the Dirac and Majorana CPV phases allowed us to derive predictions for the magnitude of the neutrinoless double beta decay effective Majorana mass, | m |, as a function of the lightest neutrino mass for all the viable cases belonging to Groups B, C and D. They are presented graphically in figures 13-16. All viable cases in the class of S 4 models investigated in the present article have distinct predictions for the set of observables sin 2 θ 12 , sin 2 θ 23 , sin 2 θ 13 , the Dirac phase δ and the absolute value of one element of the PMNS neutrino mixing matrix. Using future more precise data on sin 2 θ 12 , sin 2 θ 23 , sin 2 θ 13 and the Dirac phase δ, which will allow also to determine the absolute values of the elements of the PMNS matrix with a better precision, will make it possible to test and discriminate between the predictions of all the cases found by us to be compatible with the current data on the neutrino mixing parameters.
Future data will show whether Nature followed the S 4 H CP flavour + GCP symmetry "three-parameter path" for fixing the values of the three neutrino mixing angles and of the Dirac (and Majorana) CP violation phases of the PMNS neutrino mixing matrix. We are looking forward to these data. A Symmetry of X ν If the neutrino sector respects a residual GCP symmetry H ν CP = {X ν }, the neutrino mass matrix satisfies eq. (2.12), namely, The GCP transformation matrices X ν must be unitary due to the GCP invariance of the neutrino kinetic term. In what follows we show that these matrices are additionally constrained to be symmetric if the neutrino mass spectrum is non-degenerate, as is known to be the case. Expressing M ν from eq. (2.14) and substituting it in eq. (A.1) yields where d ν ≡ diag(m 1 , m 2 , m 3 ) andX ≡ U † ν X ν U * ν is unitary. Being 3 × 3 unitary,X can be parametrised as the product of three complex rotations U ij and a diagonal matrix of phases Ψ as follows: From the non-degeneracy of the neutrino mass spectrum it follows that sin ϑ 13 = sin ϑ 23 = sin ϑ 12 = 0. Thus,X is constrained to be diagonal and hence symmetric,X T =X. This finally implies that also X T ν = X ν , i.e., a phenomenologically relevant X ν must be symmetric.

JHEP12(2017)022
B Conjugate pairs of S 4 elements As detailed in subsection 2.2, residual flavour symmetries Z ge 2 and Z gν 2 which are conjugate to each other lead to the same form of the PMNS matrix. For G f = S 4 , there are nine group elements of order two, given in eqs.
Taking {G e , G ν } = {Z T U 2 , Z S 2 × H ν CP } with H ν CP = {U, SU } as a reference case and denoting the corresponding diagonalising matrices as Ω e and Ω ν , we find the values of φ i , δ e • , θ e • , θ ν • and k for which eq. (C.1) holds, if Ω e and Ω ν are the diagonalising matrices in one of the three remaining cases under consideration. We summarise these values in table 8.

D Correspondence with earlier results
The sum rules for cos δ or sin 2 θ 23 (sin 2 θ 12 in case C1) can formally be obtained from the corresponding sum rules derived in [7]. In certain cases, this requires an additional input which is provided by the residual GCP symmetry H ν CP considered in the present article.

JHEP12(2017)022
is found for G f = S 4 and the specific choice of the residual symmetries (see table 10 in [7]). Moreover, eq. (3.56) for cos δ can formally be obtained from the corresponding formula in case C5 of table 4 in [7] setting sin 2θe 23 = sin 2 θ e . vi) Cases C2 and D2 correspond to case C4 of [7]. The sum rule for cos δ in eq. (3.37), valid in cases C2 and D2, follows from that of case C4 in [7] (see table 4 therein) for sin 2 θ • 12 = 1/4, which is in agreement with table 10 in [7].
vii) Cases C3 and D3 correspond to case C3 in [7]. Equation (3.40) for cos δ, which holds in these cases, can be obtained from the corresponding sum rule for case C3 from table 4 in [7] with sin 2 θ • 13 = 1/4. As it should be, we find this value in table 10 of [7].
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.