Revisiting $A_4$ model for leptons in light of NuFIT 3.2

We revisit the $A_4$ model for leptons in light of new result of NuFIT 3.2. We introduce a new flavon $\eta$ transforming as $A_4$ singlet $1'$ or $1"$ which couples to both charged leptons and neutrinos in next-leading order operators. The model consists of the five parameters: the lightest neutrino mass $m_1$, the vacuum expectation value of $\eta$ and three CP violating phases after inputting the experimental values of $\Delta m_{\rm atm}^2$ and $\Delta m_{\rm sol}^2$. The model with the $1"$ singlet flavon gives the prediction of $\sin^2 \theta_{12}$ around the best fit of NuFIT 3.2 while keeping near the maximal mixing of $\theta_{23}$. Inputting the experimental mixing angles with the $1\,\sigma$ error-bar, the Dirac CP violating phase is clearly predicted to be $|\delta_\text{CP}|=50^\circ- 120^\circ$, which will be tested by the precise observed value in the future. In order to get the best fit value $\sin^2\theta_{23}=0.538$, the sum of three neutrino masses is predicted to be larger than $90\,$meV. The cosmological observation for the sum of neutrino masses will also provide a crucial test of our predictions. It is remarked that the model is consistent with the experimental data only for the normal hierarchy of neutrino masses.


Introduction
The origin of the quark/lepton flavor is still unknown in spite of the remarkable success of the standard model (SM). To reveal the underlying physics of flavors is the challenging work. The recent development of the neutrino oscillation experiments provides us with important clues to investigate the flavor physics. Indeed, the neutrino oscillation experiments have determined precisely two neutrino mass squared differences and three neutrino mixing angles. Especially, the recent data of both T2K [1,2] and NOνA [3,4] give us that the atmospheric neutrino mixing angle θ 23 favors near the maximal angle 45 • . The global analysis by NuFIT 3.2 presents the best fit θ 23 = 47.2 • for the normal hierarchy (NH) of neutrino masses [5]. The closer the observed θ 23 is to the maximal mixing, the more likely we are to believe in some flavor symmetry behind it. In addition to the precise measurements of the mixing angles, the T2K and NOνA strongly indicate the CP violation in the neutrino oscillation [2,4]. Thus, we are in the era to develop the flavor structure of the lepton mass matrices with focus on the leptonic flavor mixing angles and CP violating phase.
Before the reactor experiments measured non-zero value of θ 13 in 2012 [6,7], the paradigm of the tri-bimaximal (TBM) mixing [8,9], a highly symmetric mixing pattern for leptons, has attracted much attention. It is well known that this mixing pattern is derived in the framework of the A 4 flavor symmetry [10]- [13]. Therefore, non-Abelian discrete groups have become the center of attention in the flavor symmetry [14]- [17]. In order to obtain non-vanishing θ 13 , two of the authors improved the A 4 model by the minimal modification through introducing another flavon which transforms as 1 ( ) of A 4 and couples only to the neutrino sector [18]. Then, the predicted values of θ 13 are consistent with the experimental data. This pattern is essentially the trimaximal mixing TM 2 [19,20,21] which leads to sin 2 θ 12 ≥ 1/3. However, the predicted sin 2 θ 12 is outside of 2 σ interval of the experimental data in the NuFIT 3.2 result [5]. Therefore, the A 4 model should be reconsidered in light of the new data of T2K and NOνA as the implication of the new data has been changed.
In this work, we introduce a new flavon transforming as A 4 singlet, η (1 or 1 ) which couples to both charged leptons and neutrinos in next-leading order operators. The model consists of the five parameters: the lightest neutrino mass m 1 , the vacuum expectation value (VEV) of η and three CP violating phases after inputting the observed values of ∆m 2 atm and ∆m 2 sol . The model with a 1 singlet flavon gives the prediction of sin 2 θ 12 around the best fit of NuFIT 3.2 with keeping near the maximal mixing of θ 23 . The non-vanishing θ 13 is derived from both charged lepton and neutrino sectors. Inputting the observed mixing angles with the 1 σ error-bar, the CP violating Dirac phase is clearly predicted to be |δ CP | = 50 • − 120 • . Therefore, the observation of the CP violating phase is essential to test the model in the future.
It is remarked that the model is consistent with the experimental data only for NH of neutrino masses. The inverted hierarchy (IH) of neutrino masses is not allowed in the recent experimental data. This situation comes from that the singlet 1 or 1 flavon couples to leptons in the nextleading order. It is contrast with the model in Ref. [18] where both NH and IH are allowed.
We present our framework of the A 4 model in section 2 where lepton mass matrices and VEVs of scalars are discussed. The numerical results are shown in section 3. The section 4 is devoted to the summary. Appendix A shows the lepton mixing matrix and CP violating measures which are used in this work. The relevant multiplication rules of A 4 are represented in Appendix B. The derivation of the lepton mixing matrix is given in Appendix C. Appendix D presents the distributions of our parameter which are used in our numerical calculations.

Our framework of A 4 model
We discuss our A 4 model in the framework of the supersymmetry (SUSY). In the non-Abelian finite group A 4 , there are four irreducible representations: 1, 1 , 1 and 3. The left-handed leptons l and right-handed charged leptons e c , µ c , τ c are assigned to the triplet and singlets, respectively, as seen in Table 1. The two Higgs doublets (h u , h d ) are assigned to the A 4 singlets, and their VEVs are denoted as (v u , v d ) as usual. We introduce several flavons as listed in Table  1. The flavons φ T and φ S are A 4 triplets while ξ andξ are the same A 4 singlet 1. In addition, η andη are the same non-trivial singlet 1 or 1 . The A 4 flavor symmetry is spontaneously broken by VEVs of gauge singlet flavons, φ T , φ S , ξ and η, whereasξ (1) andη (1 , 1 ) are defined to have vanishing VEVs through the linear combinations of ξ andξ and η andη, respectively, as discussed in Ref. [13]. In the original model proposed by Altarelli and Feruglio [12,13], φ T , φ S and ξ were introduced, and then the specific vacuum alignments of the triplet flavons lead to the tri-bimaximal mixing where the lepton mixing angle θ 13 vanishes. In 2011, two of the authors minimally modified the model by introducing an extra flavon η (1 ) on top of those flavons to generate non-vanishing θ 13 [18]. This modification of the model leads to the trimaximal mixing of neutrino flavors, so called TM 2 which predicts sin 2 θ 12 ≥ 1/3 [19,20,21]. Unfortunately, this prediction for θ 12 is inconsistent with the data at 2 σ C.L. given in the NuFIT 3.2 result [5]. In this work, we force the flavon η (1 or 1 ) to couple to both charged lepton and neutrino sectors in next-leading operators by assigning a Z 3 charge to η appropriately.
We impose the Z 3 symmetry to control Yukawa couplings in both neutrino sector and charged lepton sector. The third row of Table 1 shows how each chiral multiplet transforms under Z 3 with its charge ω = exp(2πi/3).
In order to obtain the natural hierarchy among lepton masses m e , m µ and m τ , we resort to the Froggatt-Nielsen mechanism [22] with an additional U (1) FN symmetry under which only the right-handed lepton sector is charged. The field Θ denotes the Froggatt-Nielsen flavon in Table  1. The U (1) FN charges are taken as (4, 2, 0) for (e c , µ c , τ c ), respectively. By assuming that Θ, carrying a negative unit charge of U (1) FN , acquires a VEV, the relevant mass ratio is reproduced through the Froggatt-Nielsen charges.
We also introduce a U (1) R symmetry in Table 1 to distinguish the flavons and driving fields φ T 0 , φ S 0 , ξ 0 and η 0 , which are required to build a non-trivial scalar potential so as to realize the relevant symmetry breaking. In these setup, the superpotential for respecting A 4 × Z 3 × U (1) FN × U (1) R symmetry is written by introducing the cutoff scale Λ as where the subscripts 1 (1 ) in (φ T l) 1 (1 ) , etc. correspond to the case of η for 1 (1 ). The Yukawa couplings y's and y 's are complex number of order one, M is a complex mass parameter while g's and λ's are trilinear couplings which are also complex number of order one. Both leading operators and next-leading ones are included in w Y , which leads to the flavor structure of lepton mass matrices including next-leading corrections.
On the other hand, w d only contains leading operators, where we can forceξ(η) to couple with φ S 0 φ S (φ T 0 φ T ), but not ξ(η) with it sinceξ and ξ(η and η) have the same quantum numbers [13]. We can study the vacuum structure and lepton mass matrices with these superpotential.

Vacuum alignments of flavons
Let us investigate the vacuum alignments of flavons. The superpotentials w T d and w S d in Eq. (1) are written in terms of the components of triplet flavons: where w S d is the same superpotential given in Ref. [13]. Note that new terms including η andη are added in w T d .
Then, the scalar potential of the F -term is given as The vacuum alignments of φ T , φ S and VEVs of η,η, ξ andξ are derived from the condition of the potential minimum, that is V T = 0 and V S = 0 in Eq.(3) as where the VEVs ofξ andη are taken to be zero by the linear transformation of ξ andξ (η and η) without loss of generality. The coefficients λ i and g i are of order one since these flavons have no FN charges. Therefore, the VEVs of η and ξ are of same order as v T and v S , respectively. In our numerical analyses, q/Λ is scanned around v T /Λ which is fixed by the tau-lepton mass.
On the other hand, the FN flavon Θ is not contained in w d due to the U (1) FN invariance. The VEV of Θ can be derived from the scalar potential of D-term by assuming gauged U (1) FN . The Fayet-Iliopolos term leads to the non-vanishing VEV of Θ as discussed in Ref. [23]. Thus, its VEV is determined independently of v T , v S , u and q.

Lepton Mass Matrices
The explicit lepton mass matrices are derived from the superpotentials w l and w ν in Eq. (1) by use of the multiplication rule of A 4 given in Appendix B. Let us begin with writing down the charged lepton mass matrices by imposing the vacuum alignments in Eq.(4) as: where α , α η and λ are defined in terms of the VEVs of φ T , η and Θ, respectively: We note that the off-diagonal elements arise from the next-leading operators.
The left-handed mixing matrix of the charged lepton is derived by diagonalizing M M † . We obtain the mixing matrix U † approximately for the cases of η being 1 or 1 of A 4 as (more explicitly presented in Appendix C): The mass eigenvalues m 2 e , m 2 µ and m 2 τ are obtained by U M M † U † as shown in Appendix C. In the leading order approximation, U depends on one real parameter α τ η and one phase ϕ for the case of η(1 ), whereas it depends on α τ η , α µ η , ϕ and ϕ for the case of η(1 ). The parameter α η is expected to be much less than 1 as discussed in the next section. As seen in Eq. (7), the off-diagonal (1,3) and (3,1) entries in U † are dominant for the case of η(1 ) while the off-diagonal (1,2) and (2,3) (also (2.1) and (3,2)) entries in U † are dominant for the case of η(1 ). Thus, it is expected that the assignments of η(1 ) and η(1 ) give rise to different predictions of the mixing and the CP violation. It is found that the effects of the next-leading terms of O(α 2 η ), O(α 3 η ) and O(α η λ 4 ) in the mixing matrix U † are negligibly small by our numerical estimation.
The neutrino mass matrix is derived from the superpotential w ν in Eq. (1) by imposing the vacuum alignments given in Eq.(4). The next-leading operator y 5 llh u h u φ S η can be absorbed in the leading one y S llh u h u φ S due to the alignment of φ S ∝ (1, 1, 1). Although the next-leading operators llh u h u φ S φ T and llh u h u φ T ξ cannot be absorbed in the leading one, their effects are expected to be suppressed because φ T /Λ is fixed to be small. We have confirmed that the effect of those next-leading operators is negligibly small in our numerical calculations.
On the other hand, the operator y 7 llh u h u ξη leads to the significant contribution to the neutrino mass matrix because η /Λ could be significantly larger than φ T /Λ as discussed in Appendix D. For η(1 ), we have where the coefficients a, b, c and d are given in terms of the Yukawa couplings and VEVs of flavons as follows: with Since the parameter d is induced from the next-leading operator llξηh u h u , the magnitude of d is expected to be much smaller than a, b and c.
For η(1 ), we get where the last matrix of the right-hand side is a different one compared with the case of η(1 ).
There are three complex parameters in the model since the coefficient b is given in terms of a. We take a to be real without loss of generality and reparametrize them as follows: where a, c and d are real parameters and φ c , φ d are CP violating phases. For the lepton mixing matrix, Harrison-Perkins-Scott proposed a simple form of the mixing matrix, so-called TBM mixing [8,9], by which M ν is diagonalized in the case of d = 0. We obtain the neutrino mass matrix in the TBM basis by rotating it with V TBM as: where upper (lower) sign in front of (1,3) and (3,1) components correspond to η transforming as 1 (1 ). The neutrino mass eigenvalues are explicitly given in Appendix C.
The mixing matrix U ν is derived from the diagonalization ofM νM † ν apart from the Majorana phases such as As shown in Appendix C, we get where θ and σ are given in terms of parameters in the neutrino mass matrix. As seen in Eq.(10), the parameter d is related with c as where y 7 and y ξ are coefficients of order one. On the other hand, a and c are given in terms of m 1 , α ν η and the experimental data ∆m 2 sol and ∆m 2 atm as shown in Appendix C. Therefore, m 1 and α ν η are free parameters as well as φ c and φ d in our model. It is remarkable that neutrino mass eigenvalues do not satisfy ∆m 2 sol > 0 for the case of IH of neutrino masses as discussed in Appendix C because of the relation, a ∼ c and c d, in our model. It is understandable by considering the case of d = 0 limit which corresponds to the exact TBM mixing. It is allowed only for NH of neutrino mass spectrum.
Finally, the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix [24,25] is given as where P is the diagonal matrix responsible for the Majorana phases obtained from where m 1 , m 2 and m 3 are real positive neutrino masses. The effective mass for the neutrinoless double beta (0νββ) decay is given as follows: where U ei denotes each component of the PMNS matrix U PMNS , which includes the Majorana phases. From Eq. (19), we can write down the three neutrino mixing angles of Appendix A in terms of our model parameters for the case of 1 singlet η, which shows how experimental results can be accommodated in our model: where the next-leading terms are omitted. It is remarkable that sin θ 13 is composed of contributions from both the charged leptons and neutrinos. On the other hand, the deviation from the trimaximal mixing of θ 12 comes from the charged lepton sector, whereas the deviation from the maximal mixing of θ 23 comes from the neutrino sector. Since these are given in terms of four independent parameters, we cannot obtain the sum rules in the PMNS matrix elements. However, the tau-lepton mass helps us to predict the allowed region of the CP violating Dirac phase δ CP and Majorana phases α 21 and α 31 as discussed in the next section.  Table 2: The best fit, 1 σ and 3 σ ranges of neutrino oscillation parameters from NuFIT 3.2 for NH [5].
At first, we present the framework of our calculations to predict the CP violating Dirac phase δ CP and Majorana phases α 21 and α 31 . We explain how to get our predictions in terms of three real parameters α τ η , α ν η and m 1 on top of three phases ϕ, φ c and φ d for NH of neutrino masses. We can put for simplicity that is |y 7 /y ξ | = |y τ /y τ | = 1 since all Yukawa couplings are of order one. The result of NuFIT 3.2 [5] is used as input data to constrain the unknown parameters. By taking m 2 3 − m 2 1 = ∆m 2 atm and m 2 2 − m 2 1 = ∆m 2 sol with 3 σ and 1 σ data given in Table 2, a, c and d are fixed in terms of m 1 , α η , φ c and φ d . There is also the CP violating phase ϕ in the charged lepton mixing matrix. In our numerical analysis, we perform a parameter scan over those three phases and m 1 by generating random numbers. The scan ranges of the parameters are −π (ϕ, φ c , φ d ) π and 0 m 1 50 meV. Note that the range of m 1 is restricted by the lower bound of cosmological data for the sum of neutrino masses, 160 meV [26]. The parameter α η is constrained by the tau-lepton mass: which gives α = 0.0316 and 0.010 for the minimal supersymmetric standard model (MSSM) with tan β = 3 and SM, respectively. Here we put |y τ | = 1. Since α η is of same order as α as seen in Eq.(4), we vary the parameter α η around α = 0.0316 (0.010) by using the Γ distribution (χ 2 distribution), which is presented in Appendix D.
We calculate three neutrino mixing angles in terms of the model parameters while keeping the parameter sets leading to values allowed by the experimental data at 1 σ and 3 σ C.L. as given in Table 2. Then, we calculate the CP violating phases and |m ee | with those selected parameter sets. Accumulating enough parameter sets survived the above procedure, we make various scatter plots to show how observables depend on the model parameters.
In subsection 3.1, we show our numerical results for η(1 ). The numerical results for η(1 ) are briefly shown in subsection 3.2. Let us show numerical results for the case of a 1 singlet η. We analyze only the case of NH of neutrino masses since the case of IH of neutrino masses is inconsistent with experimental data as discussed in Appendix C.
At first, we show the prediction of δ CP versus sin 2 θ 23 in Fig. 1 where the blue and green dots correspond to the input of 3 σ and 1 σ data in Table 2, respectively. This result is similar to the prediction of TM 2 since the deviation from the maximal mixing of θ 23 is due to the extra (1-3) family rotation of the neutrino mass matrix in Eq. (15). In order to compare our prediction with the TM 2 result [27,28], we draw its prediction by a red curve which is obtained by taking the best fit data in Table 2. We see that our predicted region is inside of the TM 2 boundary. For the maximal mixing θ 23 = π/4, the absolute value of δ CP is expected to be 60 • -90 • . It is also predicted to be 90 • |δ CP | 110 • at the best fit of sin 2 θ 23 = 0.538. All values between −180 • and 180 • are allowed for δ CP in the case of the input data at 3 σ as seen in Fig. 1. However, for the input data at 1 σ, |δ CP | is restricted to be 50 • -120 • , which is completely consistent with the present data at 1 σ, −157 • δ CP −83 • apart from its sign. Thus, the precise data of θ 23 and δ CP would provide us with a crucial test of our prediction.
Next, we show the prediction of δ CP versus sin 2 θ 12 in Fig. 2. The deviation from the trimaximal mixing of θ 12 is due to the (1-3) family rotation of the charged lepton sector as seen in Eq. (22). The model without the additional rotation to the neutrino mass matrix in the TBM basis presented a clear correlation between sin 2 θ 12 and δ CP [27,28]. We also draw its prediction by a red curve which is obtained by taking the best fit data in Table 2. Predicted points are scattered around the red curve. Our predicted region is broad for the 3 σ data of mixing angles. However, 1 σ data forces the predicted region to be rather narrow. Then, |δ CP | = 60 • -120 • is predicted at the best fit of sin 2 θ 12 = 0.307, where the maximal CP violation |δ CP | = 90 • is still allowed.
On the other hand, we cannot find any correlation between δ CP and sin 2 θ 13 since both phases σ in the neutrino mass matrix and ϕ in the charged lepton mass matrix contribute to sin 2 θ 13 as seen in Eq. (22). We omit to present the result in a figure.  Figure 6: The allowed region on α η -δ CP plane. The meaning of colors is the same as in Fig. 1.
In order to understand the role of the key parameter α η , we show how the three neutrino mixing angles and the CP violating Dirac phase depend on α η in Figs. 3-6. At first, in Fig. 3, we show the prediction of sin θ 13 versus α η where the 3 σ data is taken as input except for sin θ 13 . The red lines denote the upper and lower bounds of the 3 σ experimental data for sin θ 13 . Note that sin θ 13 depends on α η crucially as seen in Eq. (22). As shown in Fig. 3, the observed value sin θ 13 is not reproduced unless α η is larger than 0.07. The clear dependence between α η and the predicted sin 2 θ 23 can be seen in Fig.4. In order to reproduce the maximal mixing of θ 23 , α η should be larger than 0.12. The highly probable prediction of sin 2 θ 23 is near 0.47 − 0.5 for 0.1 ≤ α η ≤ 0.2.
The deviation from the trimaximal mixing of sin 2 θ 12 explicitly depends on α η as seen in Eq. (22). We show the prediction of sin 2 θ 12 versus α η in Fig. 5. The predicted sin 2 θ 12 is almost independent of α η as far as α η ≥ 0.1.
The α η dependence on δ CP gives the characteristic prediction as shown in Fig. 6. The CP conservation δ CP = 0 is excluded in the smaller region α η ≤ 0.12 for the experimental data with 3 σ. By inputting the 1 σ data in Table 2, we obtain the prediction of δ CP to be ±(50 • − 120 • ), which is almost independent of α η for α η = 0.1-0.2.
We show the prediction of the Majorana phases α 21 and α 31 in Fig. 7. While both Majorana phases are allowed in all region of −180 • ∼ 180 • , there is a clear correlation between both phases.
In Fig. 8, we present the predicted |m ee |, the effective mass for the 0νββ decay, versus m 1 which is another key parameter in our model. The parameter m 1 should be larger than 12 meV  Figure 12: The allowed region on sin 2 θ 12 -δ CP plane for η(1 ). The meaning of colors is the same as in Fig. 1.
in order to reproduce the observed mass squared differences, and it is smaller than 46 meV due to the cosmological constraint on the sum of neutrino masses [26]. In the hierarchical case of neutrino masses m 1 < m 2 m 3 , the predicted value |m ee | is at most 10 meV but close to 45 meV for the degenerated neutrino masses.
Next, we discuss the sum of three neutrino masses Σm i because the cosmological observation gives us a upper bound for it. We show the predicted region of Σm i -sin 2 θ 23 plane in Fig. 9. The minimum of the sum of three neutrino masses Σm i is 75 meV in our model. In order to get sin 2 θ 23 ≥ 0.5, Σm i should be larger than 85 meV. For the best fit of sin 2 θ 23 = 0.538, Σm i is expected to be larger than 90 meV. We show the predicted region of Σm i -δ CP plane in Fig. 10. The predicted |δ CP | is smaller than 90 • if Σm i is smaller than 85 meV. Thus, the cosmological observation for the sum of neutrino masses will be a crucial test of these predictions.
We have neglected the next-leading terms llφ S φ T h u h u and llφ T ηh u h u in the neutrino mass matrix of Eq.(9) because α = 0.0316 (0.010) is small compared with α η ≥ 0.1. We have confirmed that those effects are small with our numerical calculation by inputting 1 σ data. Indeed, the prediction of sin 2 θ 23 -δ CP almost remains inside of the red curve in Fig. 1.
It is also deserved to comment on the α η distribution in our numerical results. In order to remove the predictions for α η > 0.3 smoothly, which is about ten times larger than α = 0.0316, we have used the Gamma distribution for α η given in Eq. (56) of Appendix D. We have confirmed that our results are not changed even if we adopt another Gamma distribution presented in Eq.(57) of Appendix D although the number density of dots gets lower. We have also used α = 0.010 which corresponds to SM in our calculations. In this case, the number density of dots significantly gets lower, but the allowed region is almost unchanged. Moreover, we have found that the allowed region is also unchanged even if we use the flat-distribution of α η in the region 0 ≤ α η ≤ 0.3. Thus, our results are robust for any distribution of α η .

Case of a 1 singlet η
We show the numerical results for a 1 singlet η briefly because the correlations of the observables appear to be weak. We show the predicted δ CP versus sin 2 θ 23 in Fig. 11. The region of |δ CP | ≤ 50 • is almost excluded while the regions near ±180 • are allowed. There is no correlation between sin 2 θ 23 and δ CP .
We also show the predicted δ CP versus sin 2 θ 12 in Fig. 12. The predicted |δ CP | increases as sin 2 θ 12 decreases from the trimaximal mixing 1/3, but its correlation is rather weak.
Both results in Figs. 11 and 12 are due to both mixing of (1-2) and (2-3) families in the charged lepton sector. Thus, the model with the 1 singlet η is less attractive than that with the 1 singlet η in light of NuFIT 3.2 data.

Summary
The flavor symmetry of leptons can be examined precisely in light of the new data and the upcoming experiments [29]. We study the A 4 model with minimal parameters by using the results of NuFIT 3.2. We introduce the A 4 singlet 1 or 1 flavon η which couples to both the charged lepton and neutrino sectors in the next-leading order due to the relevant Z 3 charge for η. The model with the 1 (1 ) flavon is consistent with the experimental data of ∆m 2 sol only for NH of neutrino masses. The key parameter is α η which is derived from the VEV of the flavon η. The parameter α η is distributed around α = 0.0316 (0.010) in the Gamma distribution of the statistic. Our results are robust for different distribution of α η .
In the case of the singlet η(1 ), α η should be larger than 0.07 in order to reproduce the observed value of sin θ 13 . The numerical prediction of δ CP versus sin 2 θ 23 is similar to the prediction of TM 2 . However, our predicted region is inside of the TM 2 boundary. The absolute value of the predicted δ CP is 60 • -90 • for the maximal mixing θ 23 = π/4. For the best fit of sin 2 θ 23 = 0.538, |δ CP | is in the region of 90 • -110 • . The predicted sin 2 θ 12 is also allowed around the best fit of NuFIT 3.2 while keeping near the maximal mixing of θ 23 . Inputting the data with the 1 σ errorbar, we obtain the clear prediction of the CP violating Dirac phase to be |δ CP | = 50 • − 120 • . The lightest neutrino mass m 1 is expected to be 12 meV-46 meV, which leads to the |m ee | < 45 meV. In order to get the best fit of sin 2 θ 23 = 0.538, the sum of the three neutrino masses is expected to be larger than 90 meV. The cosmological observation for the sum of the neutrino masses will also provide a crucial test of these predictions.
The model with η(1 ) is not attractive in light of NuFIT 3.2 result because the input data given in Table 2 does not give a severe constraint for the predicted region of δ CP .
We expect the precise measurement of the CP violating phase to test the model in the future. and two Majorana phases α 21 , α 31 as follows: where c ij and s ij denote cos θ ij and sin θ ij , respectively. The rephasing invariant CP violating measure, Jarlskog invariant [30], is defined by the PMNS matrix elements U αi . It is written in terms of the mixing angles and the CP violating phase as: where U αi denotes the each component of the PMNS matrix.

B Multiplication rule of A 4 group
We use the multiplication rule of the A 4 triplet as follow: More details are shown in the review [15,16].

C Charged lepton and neutrino mass matrices
The charged lepton mass matrix is given by the multiplication rule of A 4 in Appendix B as follows: where α , α η and λ are written in terms of the VEVs of φ T , η and Θ: The coefficients y i and y i are order one parameters. The left-handed mixing matrix of the charged lepton is derived from the diagonalization of U M M † U † . The diagonalizing matrix U † l for the charged lepton is given as follows: The mass eigenvalues of the charged leptons are given in a good approximation: where the Yukawa couplings are of order one. For η (1 ), the neutrino mass matrix is given as: where a + 3b = 0 is satisfied. The coefficients a, b, c and d are given in terms of the Yukawa couplings and the VEVs of flavons: with Taking a to be real without loss of generality, we reparametrize them as follows: where a, c and d are real parameters and φ c , φ d are CP violating phases.
We obtain the neutrino mass matrix in the TBM basis by rotating it with V TBM , as follows:M where upper (lower) sign in front of (1,3) and (3,1) components correspond to the assignment of 1 and 1 for η, respectively. Next, we consider where We obtain the neutrino mass eigenvalues for NH as follows: TheM νM † ν is diagonalized by the (1-3) family rotation as: where The θ and σ are given in terms of parameters in the neutrino mass matrix: The parameters a, c and d are written in terms of m 1 and α η . As seen in Eq. (34) , the parameter d is related with c as where y 7 and y ξ are of order one coefficients. On the other hand, a and c are given in terms of m 1 , α ν η , ∆m 2 31 ≡ m 2 3 − m 2 1 and ∆m 2 21 ≡ m 2 3 − m 2 1 since we have the following relations in Eq.(42): 1 4 (∆m 2 31 ) 2 = 3c 2 d 2 sin 2 (φ c − φ d ) + 4a 2 (c 2 cos 2 φ c + d 2 cos 2 φ d − cd cos φ c cos φ d ) , ∆m 2 21 = c 2 + d 2 + 2cd cos(φ c − φ d ) − m 2 1 .
Finally, the PMNS matrix is given as where P is the diagonal phase matrix originated from the Majorana phases. The P is obtained from P U νMν U T ν P = diag{m 1 , m 2 , m 3 }, where m 1 , m 2 and m 3 are real positive. The effective mass for the 0νββ decay is given as follows:

D Distribution of α η
The magnitude of the parameter α is determined by the tau-lepton mass as seen in Eq. (24).
The key parameter α η is related with α through the vacuum structure as discussed in Eq.(4): The coefficients λ 1(2) and g 3(4) are of order one. Then, the factor in front of α in Eq.(54) could be O(10). We scan α η by using Eq.(54) after fixing α in the statistical approach. For this purpose, we use the Gamma distribution which is available to find the distribution of the order one parameter: Taking γ = 1 with α = 3/2, µ = 0 and β = 2, we obtain which is equivalent to the χ 2 distribution. When we take γ = 2 with α = 1, µ = 0 and β = √ 2, we obtain which damps as like the Gaussian distribution at the large x.
It is easy to check that f is maximal at x = 1 and f = 0 at x = 0 for both the two types of Gamma distribution. We obtain the distribution of α η by multiplying α by f , which is used in our numerical calculations. We show the distribution of α η in Figs Figure 14: α η distribution for α = 0.0316 (blue) and α = 0.010 (red) in Eq.(57) (α = 1, β = √ 2, γ = 2, µ = 0).