Testing Lepton Flavor Models at ESSnuSB

We review and investigate lepton flavor models, stemming from discrete non-Abelian flavor symmetries, described by one or two free model parameters. First, we confront eleven one- and seven two-parameter models with current results on leptonic mixing angles from global fits to neutrino oscillation data. We find that five of the one- and five of the two-parameter models survive the confrontation test at $3\sigma$. Second, we investigate how these ten one- and two-parameter lepton flavor models may be discriminated at the proposed ESSnuSB experiment in Sweden. We show that the three one-parameter models that predict $\sin\delta_{\rm CP}=0$ can be distinguished from those two that predict $|\sin\delta_{\rm CP}|=1$ by at least $7\sigma$. Finally, we find that three of the five one-parameter models can be excluded by at least $5\sigma$ and two of the one-parameter as well as at most two of the five two-parameter models can be excluded by at least $3\sigma$ with ESSnuSB if the true values of the leptonic mixing parameters remain close to the present best-fit values.


Introduction
The flavor puzzle, which consists in understanding the number of fermion generations as well as the patterns of masses and mixing fermions obey, remains one of the most important problems in particle physics. It has no description in the Standard Model (SM) and thus calls for physics beyond the SM. In particular, this new physics is supposed to provide a mechanism for neutrino mass generation and, ideally, reveal an organizing principle behind the quark and lepton mixing patterns emerged from the data, if such a principle exists.
The tremendous experimental progress made in neutrino physics in the last two decades allowed us to establish that two leptonic mixing angles are large, while the third one is relatively small, but non-zero [1]. Such a mixing pattern appears to be very different from that in the quark sector, where all three mixing angles are small. In the attempts to quantitatively describe the peculiar structure of leptonic mixing, flavor symmetries have gained a significant interest. In particular, models based on discrete non-Abelian symmetries, which naturally allow for rotation in flavor space by fixed large angles, have been extensively studied over the past years (see Refs. [2][3][4][5][6][7][8][9] for reviews).
The ESSnuSB experiment [62,63] is a future LBL facility proposed to be built in Sweden that will significantly improve the precision on the leptonic Dirac CPV phase δ CP . The novel feature of ESSnuSB is the measurement of δ CP at the second oscillation maximum, where the measurement of δ CP is much less sensitive to systematic errors compared to the measurement at the first oscillation maximum. Therefore, it has the capability of measuring δ CP with excellent precision even with low statistics. As we will discuss, this would be of paramount importance for discriminating among various flavor models. The capability of ESSnuSB in measuring the unknown neutrino oscillation parameters within the standard three-flavor oscillation scenario has recently been studied in Refs. [64][65][66]. Other future LBL facilities, which are also capable of measuring θ 23 and δ CP with very good precision, are DUNE [67,68] and T2HK [69] that are expected to become operational in the next decade. In addition, the proposed mediumbaseline JUNO experiment [70,71] will improve the precision on the leptonic mixing angle θ 12 .
In the present work, we consider several well-motivated lepton flavor models, which lead to different mixing patterns. Our classification of these patterns is based on the number of free parameters they depend upon. First, we comment on fully-fixed mixing patterns, which do not contain any free parameter, and thus, all mixing angles (and sometimes δ CP ) are predicted to have certain values. For small (in terms of the number of elements) discrete non-Abelian groups, such as A 4 , S 4 , and A 5 , some of the angles (typically, but not universally, θ 13 ) turn out to be many sigmas away from their measured values. Next, we examine scenarios for which U PMNS depends on one real continuous parameter θ. Such mixing patterns arise from breaking of a flavor symmetry group G f combined with a generalized CP symmetry to a Z 2 × CP residual symmetry [24]. As illustrative examples, we consider the patterns derived from G f = A 4 [28], S 4 [24], and A 5 [31][32][33] combined with the CP symmetry. Further, we explore the mixing patterns obtained from breaking the same flavor symmetries, but with no CP symmetry, to a Z 2 residual symmetry in either the charged lepton or neutrino sector [20]. They are characterized by two real continuous parameters -an angle θ and a phase φ. Finally, we briefly discuss cases when the leptonic mixing matrix depends on three free parameters (either two angles and one phase or three angles).
Concentrating on the one-and two-parameter scenarios, we first confront their predictions with the current global neutrino oscillation data [72,73] (see Refs. [74][75][76] for alternative global analyses) to single out those which are well compatible with the data. For the selected models, we investigate in detail the potential of ESSnuSB to discriminate among them.
This article is organized as follows. In Section 2, we review lepton mixing patterns derived from discrete non-Abelian flavor symmetries and classify them according to the number of free parameters they depend upon, whereas in Section 3, we confront the predictions of the one-and two-parameter scenarios with global neutrino oscillation data. Then, in Section 4, we describe the ESSnuSB experimental setup, whereas in Section 5, we present the details of the simulation method used and the statistical analysis performed. Next, in Section 6, we present and discuss the results of this analysis. Finally, in Section 7, we present a summary of our work and draw our conclusions.

Lepton Mixing Patterns from Residual Symmetries
In the discrete symmetry approach to lepton flavor, it is assumed that at energies higher than a certain scale Λ, there exists a flavor symmetry described by a discrete non-Abelian group G f . The group needs to be non-Abelian, since only in this case, it has multi-dimensional (in particular, three-dimensional) irreducible representations to which three lepton SU(2) L doublets can be assigned. This in turn allows for predictions of the leptonic mixing matrix (see, e.g., Ref. [9] for more details). At energies below Λ, the flavor symmetry must be completely broken to account for three distinct charged lepton and neutrino masses. However, the charged lepton and neutrino mass matrices, M e and M ν , may be separately invariant under non-trivial Abelian subgroups G e and G ν of G f . These so-called residual symmetries constrain M e and M ν , and hence, the form of the unitary matrices U e and U ν , which diagonalize the mass matrices as follows where m i , i = 1, 2, 3, are three neutrino masses. 1 Thus, the form of the leptonic mixing matrix given by is also constrained. In the so-called standard parametrization [1], U PMNS is expressed in terms of the three leptonic mixing angles, θ 12 , θ 13 , and θ 23 , the Dirac, δ CP , and two Majorana, α 21 and α 31 , CPV phases: where c ij ≡ cos θ ij , s ij ≡ sin θ ij , P ≡ diag 1, e iα 21 /2 , e iα 31 /2 and θ ij ∈ [0, π/2] = [0, 90 • ] whilst δ, α 21 , α 31 ∈ [−π, π) = [−180 • , 180 • ). The diagonal matrix P has a physical meaning only if massive neutrinos are of Majorana nature. It can further be shown that depending on G e and G ν the leptonic mixing matrix is either completely fixed (up to permutations of rows and columns and external phases) or predicted to depend on a number of free parameters. In the following subsections, we consider lepton mixing patterns that arise from particular G e and G ν . We classify them according to the number of free parameters entering the predicted form of U PMNS .

Fully-fixed mixing patterns
If G e = Z k , k > 2 or Z m × Z n , m, n ≥ 2 and G ν = Z 2 × Z 2 , which is the maximal symmetry of the neutrino Majorana mass matrix, 2 U PMNS is completely predicted (up to permutations of rows and columns and external phases). In this case, putting it into the standard parametrization [1], one can obtain the values of the leptonic mixing parameters sin 2 θ 12 , sin 2 θ 13 , and sin 2 θ 23 , as well as the Dirac CPV phase δ CP (if θ 13 = 0).
A complete classification of all possible mixing patterns fully-fixed by the residual symmetries G e and G ν has been performed in Ref. [90]. From 17 sporadic cases and one infinite series of mixing matrices found therein, only the latter can lead to phenomenologically viable values of the leptonic mixing angles. Namely, this infinite series is given (up to permutations of rows and columns) by 3 where σ = e 2πip/n and ω = e 2πi/3 are roots of unity and p and n are coprime. Thus, for all viable mixing patterns, we have the following prediction Furthermore, sin δ CP = 0, i.e., no Dirac CP violation is predicted. A given choice of σ defines the flavor symmetry group G f . It turns out that the smallest viable group leading to n = 9, 18 is (Z 18 × Z 6 ) S 3 , which has 648 elements. The next "minimal" group implying n = 11, 22, 33, 66 is ∆(6 × 22 2 ), which has order 2904. Such groups are much more complex (less natural) than those originally proposed for description of the observed lepton mixing pattern.

Models with one free parameter
A discrete flavor symmetry can be consistently combined with a generalized CP symmetry [24][25][26]. The latter is given by a CP transformation, which can act non-trivially in flavor space. This action is represented by a unitary symmetric matrix X. The full symmetry group, in this case G CP , is given by a semi-direct product G CP = G f CP [24]. Breaking this symmetry to G e = Z k , k > 2 or Z m × Z n , m, n ≥ 2 and G ν = Z 2 × CP leads to U PMNS that is defined up to a rotation R ij (θ) in the neutrino sector. Here, ij refers to the plane of rotation, e.g., where θ is a free real continuous parameter and its fundamental interval is [0, π) = [0, 180 • ) [24]. The advantage of this approach is that the Majorana CPV phases are also predicted (contrary to the case without CP).
In what follows, we consider the lepton mixing patterns originating from S 4 CP [24] and A 5 CP [31][32][33] broken to the residual symmetries as described above. Note that the results | sin δ CP | 1 0 1 0 Table 1: Predictions for the leptonic mixing parameters in terms of the free parameter θ in Cases I, II, IV, and V originating from G CP = S 4 CP broken to G e = Z 3 and G ν = Z 2 × CP derived in Ref. [24].
obtained from A 4 CP [28] are contained in those for S 4 CP [24]. Thus, we do not need to consider the A 4 CP case separately. In Ref. [24], it has been shown that fixing G e = Z 3 , there are five inequivalent choices of the Z 2 × CP transformations leaving the neutrino sector invariant (due to different Z 2 subgroups of S 4 and CP transformations X compatible with them). Four of them, denoted Cases I, II, IV, and V, lead to the results summarized in Table 1. Case III has been found to be phenomenologically not viable and we do not present it here. Cases I and II realize the so-called trimaximal (TM) mixing pattern 2 (TM 2 ) [91], for which |(U PMNS ) 2 | 2 = 1/3, = e, µ, τ , while Cases IV and V give rise to TM 1 [92], characterized by |(U PMNS ) e1 | 2 = 2/3 and |(U PMNS ) µ1 | 2 = |(U PMNS ) τ 1 | 2 = 1/6. In particular, Cases I and II lead to Eq. (2.6), 4 i.e., they predict sin 2 θ 12 ≥ 1/3, while for Cases IV and V sin 2 θ 12 = 1 − 3 sin 2 θ 13 3 1 − sin 2 θ 13 < 1 3 . (2.8) Furthermore, the following sum rule for cos δ CP holds in Cases I and II cos δ CP = 1 − 2 sin 2 θ 13 cot 2θ 23 sin θ 13 2 − 3 sin 2 θ 13 . (2.9) In Case I, sin 2 θ 23 is predicted to be 1/2, and thus, cos δ CP = 0. To establish for which values of θ the Dirac CPV phase δ CP = −90 • and for which it is 90 • , one needs to look at the rephasing invariant J CP [93] that controls the magnitude of CPV effects in neutrino oscillations [94]. In terms of the elements of U PMNS , it can be chosen as sin δ CP 0 Table 2: Predictions for the leptonic mixing parameters in terms of the free parameter θ in the case originating from G CP = S 4 CP broken to G e = Z 4 (or Z 2 × Z 2 ) and G ν = Z 2 × CP derived in Ref. [24].
For G e = Z 4 , only one viable case has been found. It leads to the predictions shown in Table 2. The existence of two solutions, marked as Cases VI-a and VI-b, arises from the freedom to exchange the second and third rows of the leptonic mixing matrix. For both solutions, we have 19) and in addition, we find that (ϕ cos θ+sin θ) 2 4ϕ 2 −(cos θ−ϕ sin θ) 2 Table 3: Predictions for the leptonic mixing parameters in terms of the free parameter θ in Case V (VII) originating from G CP = A 5 CP broken to G e = Z 3 (Z 2 ×Z 2 ) and G ν = Z 2 ×CP derived in Ref. [31]. The quantity ϕ = (1 + √ 5)/2 is the golden ratio. | sin δ CP | 1 0 1 Table 4: Predictions for the leptonic mixing parameters in terms of the free parameter θ in Cases II, III, and IV originating from G CP = A 5 CP broken to G e = Z 5 and G ν = Z 2 × CP derived in Ref. [31]. The quantity ϕ = (1 + √ 5)/2 is the golden ratio.
Note that these equations are consistent with the fact that under exchange of the second and the third rows of U PMNS , θ 12 and θ 13 remain unchanged, while θ 23 → π/2 − θ 23 and δ CP → δ CP + π. The J CP invariant vanishes, which means that |cos δ CP | = 1. Substituting the expressions for sin 2 θ ij given in Table 2 in Eqs. (2.20) and (2.21), we find cos δ CP = ± sgn (1 + 3 cos 2θ) √ 2 sin θ − cos θ , (2.22) where "+" corresponds to Case VI-a and "−" to Case VI-b. This in turn implies that Similarly, lepton mixing patterns originating from G CP = A 5 CP broken to G ν = Z 2 ×CP in the neutrino sector have been derived in Refs. [31][32][33]. The results for G e = Z 3 (Case V) and G e = Z 2 × Z 2 (Case VII) are presented in Table 3 and those for G e = Z 5 (Cases II, III, and IV) in Table 4. The case numbering follows Ref. [31]. Several comments are in order. First, Cases I, VI, and VIII have been found to be not viable, so we do not present them here. Secondly, since θ is a free parameter and the corresponding rotation can be performed either clockwise or counterclockwise, we can replace θ → π/4 ± θ. Then, the results for Case V in Table 3 match those for Case I in Table 1. Thus, these two cases are identical (cf. Ref. [33]). Further, two solutions in Case VII arise from the exchange of the second and the third rows of the leptonic mixing matrix. Finally, the predictions in all cases except for Case V involve the golden ratio ϕ characteristic for the A 5 group.
To establish the values of θ leading to δ CP = 0 (−90 • ) and those giving rise to δ CP = 180 • (90 • ) in Cases VII and III (Cases II and IV), one needs to perform an analysis similar to that we have described for G CP = S 4 CP. In this work, we do not present the corresponding analytical expressions for the sake of brevity. However, we have established such a correspondence and used it in our main analyses in Sections 3 and 6.
It is worth noting that all cases in Tables 1-4 predict either CP conservation or maximal CP violation. In addition, the cases predicting the latter also lead to sin 2 θ 23 = 1/2. Finally, for these cases, sin 2 θ 12 and sin 2 θ 13 remain invariant under θ → π − θ, whereas δ CP changes from −90 • to 90 • . This fact will manifest itself in Section 3, where we fit the models to global neutrino oscillation data.

Models with two free parameters
Now, if we relax the assumption of generalized CP invariance and break the flavor symmetry group G f to either the leptonic mixing matrix U PMNS is defined up to a complex rotation U ij (θ l ij , δ l ij ) in either the charged lepton (l = e) or neutrino (l = ν) sector and thus contains two free real continuous parameters. Here, ij refers to the plane in which the rotation is imposed, e.g., All possible cases of this type have been considered in Ref. [20], where they have been classified according to the plane in which the free rotation is performed. Below, we summarize the expressions for sin 2 θ ij , cos δ CP , and J CP in terms of the two free parameters. We will use them further in our statistical analysis. Case A1. In this case, a fixed part of the leptonic mixing matrix parametrized in terms of the angles θ • ij defined by the choice of residual symmetries G e , G ν ⊂ G f , is corrected from the left by U 12 (θ e 12 , δ e 12 ). The leptonic mixing parameters and J CP are given by where θ ∈ [0, π) = [0, 180 • ) and φ ∈ [0, 2π) = [0, 360 • ) are free parameters related to θ e 12 and δ e 12 (see Appendix B of Ref. [20] for details). From Eqs. (2.26) and (2.28), we see that in this case there are two relations not explicitly involving the free parameters: (i) between sin 2 θ 23 and sin 2 θ 13 and (ii) among cos δ CP , θ 12 , and θ 13 . This can be expected, since four observables (sin 2 θ ij and δ CP ) have been expressed in terms of two parameters (θ and φ). To obtain the value of δ CP for given (θ, φ), the expression for J CP in Eq. (2.29) has to be compared to its expression in the standard parametrization given in Eq. (2.10). We note that sin 2 θ ij in Eqs. (2.25)-(2.27), and hence, cos δ CP in Eq. (2.28), depend on θ and φ via sin 2 θ and the product sin 2θ cos φ. There are two transformations that leave them invariant: Under both of these transformations, J CP changes sign, and as a consequence, also δ CP → −δ CP . Thus, for given (θ, φ), we have four solutions that lead to the same values of sin 2 θ ij , of which two -(θ, φ) and (π − θ, π + φ) -give δ CP and the other two -(π − θ, π − φ) and (θ, 2π − φ) -yield −δ CP . This observation applies to all two-parameter cases considered below.
Case A2. In this case, the corresponding free rotation matrix is U 13 (θ e 13 , δ e 13 ). The leptonic mixing angles, cos δ CP , and J CP are defined by The free parameters θ and φ are related to θ e 13 and δ e 13 . As in the previous case, we find sum rules for sin 2 θ 23 and cos δ CP .
Case A3. The correction to a fixed part of U PMNS due to U 23 (θ e 23 , δ e 23 ) leads to (2.40) Therefore, while the mixing angles θ 12 and θ 13 are predicted to have certain fixed values, the mixing angle θ 23 and the Dirac CPV phase δ CP remain unconstrained in this case. Case B1. For pattern B, the leptonic mixing matrix is defined up to a free complex rotation from the right. If this rotation is U 13 (θ ν 13 , δ ν 13 ), one has sin 2 θ 13 = cos 2 θ • 12 sin 2 θ , (2.41) where θ and φ are related to θ ν 13 and δ ν 13 and the angles θ • ij are fixed once the residual symmetries G e and G ν originating from breaking G f are specified. This case is characterized by the sum rules for sin 2 θ 12 and cos δ CP , i.e., Eqs. (2.43) and (2.44), respectively.
Case B3. For U 12 (θ ν 12 , δ ν 12 ), the following simple relations hold i.e., while the mixing angles θ 13 and θ 23 are predicted to have certain fixed values, the mixing angle θ 12 and the Dirac CPV phase δ CP remain unconstrained. It is worth noting that Cases A1, A2, B1, and B2 lead to non-trivial predictions for the Dirac CPV phase δ CP . The study of these cases for G f = A 4 , S 4 , and A 5 and all their Abelian subgroups, which can play the role of the residual symmetries G e and G ν , has been performed in Refs. [20,21] in light of global neutrino oscillation data. For G f = S 4 , only two viable cases have been found [21]. We display the corresponding residual symmetries and the values of the parameters sin 2 θ • ij fixed by them in Table 5. For G f = A 5 , there are six viable cases [21], which we summarize in Table 6. Irrational values of sin 2 θ • ij quoted therein can be expressed Table 5: Values of the fixed parameters sin 2 θ • ij for viable cases originating from G f = S 4 broken to the residual symmetries G e and G ν , as derived in Ref. [20] (see also Ref. [21]). The entries marked with "−" are not relevant. Table 6: Values of the fixed parameters sin 2 θ • ij for viable cases originating from G f = A 5 broken to the residual symmetries G e and G ν , as derived in Ref. [20] (see also Ref. [21]). The entries marked with "−" are not relevant.

Models with three free parameters
Here, we briefly discuss the scenarios for which U PMNS depends on three free parameters. For instance, this happens if G f is broken to G e = Z 2 and G ν = Z 2 . In this case both U e and U ν are defined up to complex rotations. Thus, naively, we have four free parameters -two angles and two phases. However, as demonstrated in Ref. [20], the number of free parameter reduces to three -two angles and one phase -after an appropriate rearrangement. Since four observables are now expressed in terms of three parameters, one relation between the observables can be expected. Indeed, depending on the planes in which the two complex rotations are imposed, either cos δ CP is determined by the mixing angles or one relation among the latter is found.
Another example of a three-parameter setup is given by G e = Z 2 and G ν = Z 2 × CP. In this case, the free rotation in the neutrino sector is real. Such a breaking pattern has been investigated in Ref. [37] for G CP = S 4 CP and in Ref. [34] for G CP = A 5 CP. Furthermore, if G e is larger than Z 2 and a single residual CP transformation is preserved in the neutrino sector, U PMNS depends on three free angles [36].
Finally, there exists a possibility that G f is fully broken in the charged lepton (neutrino) sector, while the form of the matrix U ν (U e ) is fully determined by a residual G ν (G e ) symmetry (which is larger than Z 2 ). In such a case, the unitary matrix U e (U ν ) is unconstrained, unless its form is motivated by some additional arguments. A well-studied example is given by  U † e = U 12 (θ e 12 , δ e 12 )U 23 (θ e 23 , δ e 23 ) and U ν being the BM, TBM, or GR mixing matrix. In this case, U e provides necessary corrections to U ν , in particular, generating a non-zero θ 13 and shifting θ 12 from its leading order value θ ν 12 . 6 Such scenario leads to the following sum rule [17] cos δ CP = tan θ 23 sin 2θ 12 sin θ 13 cos 2θ ν 12 + sin 2 θ 12 − cos 2 θ ν It has been studied in detail in Refs. [18,43] and the effects of renormalization group evolution of the leptonic mixing parameters on its predictions have been further investigated in Refs. [95,96]. In Ref. [47], this sum rule has been thoroughly analyzed in the context of DUNE and T2HK. Alternative sum rules arising from the charged lepton corrections have been derived in Ref. [19].
In the next section, we will confront eleven one-parameter mixing patterns (see Tables 1-4) and seven two-parameter ones (see Tables 5-6) with the current global neutrino oscillation data. Note that we will not consider the two-parameter Cases A3 and B3, since they do not lead to predictions for cos δ CP [cf. Eqs. (2.39) and (2.54)]. Finally, we will consider neither fully-fixed mixing patterns nor those depending on three free parameters. While the former to be viable requires very large discrete groups (see Subsection 2.1), the predictive power of the latter is less than that of models with one or two free parameters.

Confronting the Flavor Models with Global Neutrino Oscillation Data
Fits to global neutrino oscillation data have basically been performed by three groups (see Refs. [72,74,76]). In our analysis, we use the results of the so-called NuFIT group announced in July 2019 [73]. In Table 7, we list the current best-fit values of the leptonic mixing parameters as well as their corresponding 1σ errors and 3σ ranges. In order to confront the one-and two-parameter lepton flavor models discussed in Section 2 with global neutrino oscillation data and to determine the allowed values of the model parameters, we define the χ 2 function as follows χ 2 (θ) = sin 2 θ 12 (θ) − sin 2 θ 12 σ(sin 2 θ 12 ) 2 + sin 2 θ 13 (θ) − sin 2 θ 13 σ(sin 2 θ 13 ) 2 + sin 2 θ 23 (θ) − sin 2 θ 23 σ(sin 2 θ 23 ) 2 (3.1) 6 We abuse the notation θ ν 12 , which here denotes the value of θ12 for BM, TBM, or GR mixing, and not a free parameter as in Subsection 2.3.
In Table 8, we present the results of our fits including the minimal values of the χ 2 function (i.e., χ 2 min ) for the eleven relevant one-parameter models and the seven relevant two-parameter models and rearrange them according to how well they fit the data of NuFIT 4.1. In addition, we give the best-fit values of the corresponding model parameters. In the case of the oneparameter models, the best-fit parameter and its corresponding 3σ interval are denoted θ bf and θ 3σ , respectively, whereas in the case of the two-parameter models, the best-fit parameters and their corresponding 3σ intervals are denoted by θ bf , θ 3σ and φ bf , φ 3σ , respectively. To calculate the 3σ intervals for both one-and two-parameter models, we impose χ 2 − χ 2 min = 3 2 = 9, corresponding to one degree of freedom (d.o.f.). In fact, this holds also for the two-parameter models, since presenting the 3σ interval for one of the two parameters, we minimize χ 2 over the other parameter. Note that the model parameters θ and φ have different meaning for different models. In Table 9, as an outcome of the minimization of the χ 2 function, we display the bestfit values of the leptonic mixing parameters sin 2 θ 12 (bf), sin 2 θ 13 (bf), and sin 2 θ 23 (bf), which are evaluated at the best-fit values of the model parameters. Note that some of the models have fixed values of some of the leptonic mixing parameters, which stems from the fact that these parameters are independent from the corresponding model parameters. For example, sin 2 θ 23 (bf) = 1/2 for Models 1.3, 1.5, 1.6, and 1.9. Furthermore, we calculate the predicted values of the CPV phase δ CP for the various one-and two-parameter models at the best-fit values of the corresponding model parameters (given in Table 8), which are also presented in Table 9. We observe that all one-parameter models have fixed values of δ CP , which are either 0 and 180 • (CP conservation) or −90 • and 90 • (maximal CP violation) and these values are, by construction, given by the respective model. For the two-parameter models, it is interesting to note that the negative values of δ CP for Models 2.3-2.5 and 2.7 are within the 1σ interval of (−167 • , −100 • ) [73].
In Figs. 1 and 2, we show the predictions of the leptonic mixing parameters in different planes for the eleven one-and the seven two-parameter models, respectively. To produce Fig. 1, we vary the model parameter θ in the 3σ interval for the one-parameter models as given in Table 8 and obtain the values of the leptonic mixing angles for each value of θ using the relations given in Tables 1-4. We plot the obtained leptonic mixing parameters against each other. Since all relations between θ and the angles are rather simple, this method is sufficient to describe the allowed parameter space of the one-parameter models. On the other hand, for the two-parameter models, some of the relations between the free parameters θ and φ and the leptonic mixing parameters are simpler and some of them are more complex. Therefore, to generate Fig. 2, we generate values of θ and φ with a probability density proportional to e −χ 2 /2 for each model, where χ 2 is the function given in Eq. (3.2). We then draw the smallest possible regions in the predicted parameter space containing 95 % of the generated combinations. For the cases, where the models predict a direct relation between the two parameters plotted, e.g., if θ 23 can be written as a function of θ 13 only, we draw a curve on which 95 % of the generated parameter combinations lie.   [73]. See also Table 7. Note that there are four degenerate best-fit points (θ bf , φ bf ) for each two-parameter model (i.e., Models 2.1-2.7) that lead to the same value of the χ 2 min function. Two best-fit points are presented for each model and the other two can be found performing the replacement rule φ bf → 2π − φ bf . In addition, the two presented best-fit points for each model are related by the simultaneous replacement of θ bf → π − θ bf and φ bf → π − φ bf .   Figure 1: Prediction of leptonic mixing angles for the eleven one-parameter models. The lines show the allowed values of the leptonic mixing angles for the 3σ intervals θ 3σ of the respective model parameters θ given in Table 8. The white areas show the allowed 3σ regions of the leptonic mixing angles from fits to global neutrino oscillation data, whereas the gray-shaded areas show the corresponding excluded regions. A star (" ") indicates the present best-fit values of the leptonic mixing angles from global neutrino oscillation data. Top panel: sin 2 θ 13 -sin 2 θ 12 . Middle panel: sin 2 θ 13 -sin 2 θ 23 . Bottom panel: sin 2 θ 12 -sin 2 θ 23 .
From Fig. 1 for the eleven one-parameter models, we note that the allowed values of (i) sin 2 θ 12 (i.e., the values within the 3σ intervals of the model parameter) for Models 1.  Table 8   corresponding values are all below 3σ. Therefore, we deem Models 1.6-1.11 as excluded at 3σ or more by the current data. Let us understand why, among the allowed one-parameter models (i.e., Models 1.1-1.5), Model 1.1 gives the best fit to the data and Model 1.5 the worst. From Fig. 1, we note that these models do not give any restrictions on sin 2 θ 13 (except from Model 1.4), and therefore, the fits mainly depend on the pulls of the predictions for sin 2 θ 12 and sin 2 θ 23 . The distance between the prediction of a given model and the best-fit point in the model parameter space basically determines the value of χ 2 min of that model in terms of the pulls of the three different leptonic mixing angles. Thus, this distance increases with increasing model number, and therefore, Model 1.1 leads to the best fit and Model 1.5 to the worst.
Next, from Fig. 2 for the seven two-parameter models, we observe that all seven models predict leptonic mixing angles within the 3σ regions of the global fit of NuFIT 4.1, since none of the contours displayed in the six panels lie totally outside these 3σ regions. Furthermore, it is interesting to note that Models 2.6 and 2.7 predict the lower octant (LO) of θ 23 (see top-right, middle-left, and bottom-right panels), whereas Models 2.1-2.5 predict the higher octant (HO). In addition, the intervals of δ CP predicted by the different two-parameter models are: Here, we also note that the two-parameter models do not put any restrictions on sin 2 θ 13 , and therefore, the fits depend on their predictions for sin 2 θ 12 and sin 2 θ 23 . From the top-right panel of Fig. 2, we see that Model 2.1 predicts sin 2 θ 12 and sin 2 θ 23 closest to the current best-fit point, whereas Model 2.7 predicts sin 2 θ 12 and sin 2 θ 23 furthest from the current best-fit point. This is the reason why Model 2.1 gives the best fit to the data and Model 2.7 the worst. From this panel, it is also easy to understand why the fits of Models 2.2-2.6 lie in between the fits of Models 2.1 and 2.7. Note that among the seven two-parameter models, Models 2.6 and 2.7 have values of χ 2 min larger than 9, i.e., 3σ for 1 d.o.f. Thus, confronting the 18 one-and two-parameter lepton flavor models with global neutrino oscillation data by NuFIT 4.1, we conclude that Models 1.1-1.5 and 2.1-2.5 are allowed at 3σ, whereas Models 1.6-1.11, 2.6, and 2.7 are excluded at 3σ or more by the current data. For the analysis of the models with ESSnuSB, we consider only those models that are allowed by the current data at 3σ. Therefore, in the next three sections, we will address these ten models, i.e., Models 1.1-1.5 and 2.1-2.5, with ESSnuSB.

Experimental Setup of ESSnuSB
The ESSnuSB experiment [62,63] is a proposed long-baseline neutrino oscillation experiment in Sweden and mainly designed to measure potential leptonic CP violation with excellent precision. In our work, we use exactly the same configuration of ESSnuSB as was used to generate the results of Ref. [66]. In the original experimental setup of ESSnuSB, the second peak of the neutrino oscillation probability was identified to be optimal to maximize the CP violation discovery potential. The source of neutrinos is a beam of power 5 MW capable of delivering protons of energy 2.5 GeV corresponding to 2.7 × 10 23 protons on target per year, situated at the European Spallation Source (ESS) in Lund, Sweden. The detector is 1 Mt MEMPHYS-like water-Cherenkov detector [97] located at a distance of 540 km away from the source in the mine in Garpenberg, Sweden. We assume a total running time of 10 years with 5 years in neutrino mode and 5 years in antineutrino mode. In our analysis, we also consider an identical near detector having mass of 0.1 kt and located at a distance of 500 m away from the source. For the near detector, the fluxes are simulated at a distance of 1 km from the  target station, whereas for the far detector, they are calculated at 100 km. We use correlated systematics between the near and the far detectors. The specification for the systematic errors is adopted from Ref. [98] and listed in Table 10 for convenience.

Details of Simulation and Statistical Analysis
All of our numerical simulations are performed using the GLoBES software [99,100] with the experimental description implemented as described in Section 4. In order to judge whether or not two models can be distinguished using the data of ESSnuSB, we follow the statistical procedure laid out in this section. In order to accomplish this, we introduce the chi-square function where χ 2 ESSnuSB is the chi-square function for the ESSnuSB data alone, χ 2 0 is the prior chisquare coming from already performed experiments, θ is the set of parameter values in a particular model, and D is the data obtained in ESSnuSB. For χ 2 ESSnuSB , we use the default GLoBES chi-square function where D i is the number of events in bin i andD i (θ) is the prediction for bin i of the model being tested given the parameter set θ. We assume χ 2 0 to be the same prior as used in Section 3. In order to answer the question of whether Model A can be excluded if Model B is true, we assume that the data gathered in ESSnuSB will be the Asimov data [101] corresponding to Model B being true with particular parameter values θ B true . We then minimize χ 2 (θ A , D) for the predictions of Model A in order to find the best-fit parameters of Model A and take this χ 2 as a measure of whether or not Model A could be ruled out if Model B were true. In the cases, where the true Model B is taken to be a one-parameter model, we present our results as a function of the true parameter θ B true , and when it is taken to be a model with two parameters, our results are presented as the profiled χ 2 , i.e., as a function of one of the true parameters minimized over the other.
Note that, since we include the current priors in the χ 2 function, it may occur that even if Model B is true, Model A provides a better fit to the ESSnuSB data together with the prior. This will be particularly true when the assumed true Model B is disfavored relative to Model A by current data. The interpretation of such a result should therefore be that the data of ESSnuSB will not be able to provide sufficient discriminatory power to overcome the fact that Model B is disfavored by current data. On the other hand, if χ 2 of the test Model A is significantly larger than that of the assumed true Model B, then Model A could be ruled out by the ESSnuSB result if Model B is the model implemented in Nature.
In addition to testing models against each other, we simulate the case where the oscillation parameters take particular values and consider the level at which the different models could be rejected in such a case. In order to do this, we follow the very same procedure as outlined above with the difference that all models are compared to the global χ 2 min , i.e., we minimize χ 2 with respect to all of the leptonic mixing parameters (θ 12 , θ 13 , θ 23 , and δ CP ) and compare this to χ 2 min in each model. Note that the number of d.o.f. is taken to be 3 for the one-parameter models and 2 for the two-parameter models as the models are being compared to the case where all four mixing parameters (three angles and δ CP ) are left free. In all cases, we fix the mass-squared differences to their current best-fit values from Table 7. 6 Results: Addressing Flavor Models with ESSnuSB In Fig. 3, we show the comparison of the one-parameter models 1.1-1.5, where each panel corresponds to different assumed true models and the horizontal axes represent the assumed true parameter value in the true model. The different curves show the resulting χ 2 min in the different models, including both the prior and the Asimov data from ESSnuSB. In each panel, the curve of the assumed true model is drawn in black in order to highlight it. Note that, due to the addition of the prior, χ 2 of the true model is not equal to zero despite using the Asimov data for ESSnuSB. For each model, χ 2 min is therefore equal to the minimum of the prior for that model. This also means that it is possible for an assumed true model to give a worse fit than another model, even in the case where the ESSnuSB data are generated from that model. The interpretation of this should be that the ESSnuSB data will not be sufficient to overcome the preference for the other model that is present in the current data.
For the one-parameter models, we see from Fig. 3 that the models essentially split into two different groups, one group including Models 1.1, 1.2, and 1.4 and another including Models 1.3 and 1.5. The reason is that Models 1.1, 1.2, and 1.4 predict sin δ CP = 0, whereas Models 1.3 and 1.5 predict | sin δ CP | = 1. If the ESSnuSB data is generated assuming a model from one group, then χ 2 of all models in that group remains more or less the same, whereas the models from the other group are strongly disfavored due to the ESSnuSB precision to δ CP . For example, looking at the top-left panel, where Model 1.1 is assumed to be true, χ 2 of Models 1.1, 1.2, and 1.4 remain essentially the same regardless of the assumed true value of the parameter θ, while the other models would be excluded at a level of around 7σ, where we use ∆χ 2 , where ∆χ 2 is the difference in χ 2 min between the models, as the sensitivity estimator [102]. Correspondingly, Model 1.1 would be strongly excluded when data is generated from any of Models 1.3 and 1.5.
In Figs. 4-8, we present the results for the two-parameter models. As the two-parameter models have several parameters, each panel shows the results for particular choices of the true parameter shown on the horizontal axes profiled over the other parameter, e.g., in Fig. 4  (upper-left panel), the curves show how χ 2 min of the different models behave as a function of the true θ in Model 2.1 when minimized over the true value of φ. Unlike in the case of the one-parameter models, the two-parameter models do not split into different groups. Still, χ 2 of the different models are highly dependent on the assumed true model, but also on the true parameter values assumed for those models. A remarkable feature is that, for the   Figure 3: Comparison among the five allowed one-parameter models as test models for ESSnuSB using priors. Each panel shows the quantity χ 2 min , given a true model, as a function of the model parameter θ in its 3σ interval θ 3σ presented in Table 8. In each panel, χ 2 min for the given true model is displayed with a black solid curve (when acting as a test model), whereas it is marked for the five test models by curves corresponding to their respective model colors.   Figure 4: Comparison among the five allowed two-parameter models as test models for ESSnuSB using priors in the cases the model parameters φ (upper row) and θ (lower row) are minimized and corresponding to the first (left column) and the second (right column) best-fit points in Table 8  three models that provide the best fit to the current neutrino oscillation data, the ESSnuSB data make those models the preferred one if they are assumed as the true model for most of the parameter space. As such, the ESSnuSB data can also help in the rejection of other models, in particular those that are relatively close in χ 2 at the present time, even for the two-parameter models where the regions of the leptonic mixing parameters are more extended than for the one-parameter models, cf. Figs. 1 and 2. The addition of the ESSnuSB data is sufficient to disfavor most other models at at least 3σ for Models 2.1-2.3. As discussed earlier, the two-parameter models give degenerate fits of θ and φ to current data, corresponding to δ CP < 0 or δ CP > 0, since the values of the other leptonic mixing parameters are the same. In other words, the two-parameter models are symmetric for the degenerate values of θ and φ apart from the prediction for δ CP and this feature is reflected in Figs. 4-8, where the exact symmetry is broken by the different phenomenology for different values of δ CP at ESSnuSB. For example, the symmetry is more prominent in Fig. 4 and less prominent in Fig. 6. This can be understood, since for Model 2.1, the neutrino oscillation probabilities are not very different when δ CP → −δ CP , while they are different for Model 2.3. Another interesting feature to note is that the sensitivities for Models 2.4 and 2.5 are very similar. The reason is that these two models predict similar values of all leptonic mixing parameters except for θ 12 .
It is also possible to compare models with different number of parameters. As performing the full analysis using all ten models from the one-and two-parameter cases would result in   Figure 5: Comparison among the five allowed two-parameter models as test models for ESSnuSB using priors in the cases the model parameters φ (upper row) and θ (lower row) are minimized and corresponding to the first (left column) and the second (right column) best-fit points in Table 8  figures that would be extremely cluttered, we present the comparison of the one-and twoparameter models that are the current best fits to the available data, i.e., Models 1.1 and 2.1.
The resulting comparisons are shown in Figs. 9 and 10 assuming Model 1.1 and 2.1 to be the true model, respectively. From these figures, we can see that the two models could be separated by around 6σ or more, depending on the true model parameters. The reason for this is that Model 1.1 predicts a value of δ CP of 180 • , whereas the values predicted by Model 2.1 are spread, but closer to 0 than 180 • . Thus, the measurement of δ CP by ESSnuSB will be able to provide a good discriminator between the models. In Fig. 11, we present the capability of ESSnuSB to exclude the one-parameter models in the sin 2 θ 23 (true)-δ CP (true) plane. For the other neutrino oscillation parameters, we assume the following true values sin 2 θ 12 = 0.310, sin 2 θ 13 = 0.02237, ∆m 2 21 = 7.39 · 10 −5 eV 2 , and ∆m 2 31 = 2.528 · 10 −3 eV 2 . If the true values of the leptonic mixing parameters chosen by Nature fall in the shaded region and the Asimov data is assumed, then the model under test will be excluded at the shown confidence level. We observe that, for all five models, no model would survive at 1σ and the region of true parameters for which the models would survive at 3σ shrink in size as we progress from Model 1.1 to Model 1.5. The reason is the following: as we go from Model 1.1 to Model 1.5, the fit of the models to the current data becomes worse, and therefore, they can be excluded at higher confidence level by ESSnuSB, in particular as they are relatively bad fits at the present time with the current best-fit values for the   Figure 6: Comparison among the five allowed two-parameter models as test models for ESSnuSB using priors in the cases the model parameters φ (upper row) and θ (lower row) are minimized and corresponding to the first (left column) and the second (right column) best-fit points in Table 8  leptonic mixing parameters. Since Models 1.1, 1.2, and 1.4 predict δ CP = 180 • , the regions where these models cannot be ruled out are around δ CP (true) = 180 • , and since Models 1.3 and 1.5 predict δ CP = ±90 • , the regions remain close to δ CP (true) = ±90 • , which is due to the excellent δ CP measurement capability of ESSnuSB. Therefore, if the true value of δ CP is close to the present best-fit point (indicated by stars), then Models 1.1, 1.2, and 1.4 can be excluded at more than 5σ, whereas Models 1.3 and 1.5 will still be consistent with the current data at 5σ. From the different panels, we also note that for Models 1.3-1.5, for some of the true values of sin 2 θ 23 (true) the models could also be excluded by ESSnuSB at 3σ, but this is not the case at 5σ. This is due to the relatively poor capability of ESSnuSB to measure the leptonic mixing angle θ 23 . Furthermore, we have checked that the global χ 2 min is larger at δ CP (true) = 90 • and smaller at δ CP (true) = −90 • in the lower octant of θ 23 (true). For Models 1.3 and 1.5, this is the reason why the exclusion around δ CP (true) = −90 • is stronger than around δ CP (true) = 90 • .
Finally, in Fig. 12, we present the same as in Fig. 11, but for the two-parameter models. For the same reason as discussed for the one-parameter models, as we progress from Model 2.1 to Model 2.5, the regions where the models cannot be ruled out become smaller. We observe that for Models 2.1 and 2.2, there is a region where the models cannot be ruled out even at 1σ, whereas for Models 2.3-2.5, such regions do not exist. If the true values of θ 23 and δ CP remain close to the current best-fit values, then Models 2.1 and 2.2 would be compatible with   Figure 7: Comparison among the five allowed two-parameter models as test models for ESSnuSB using priors in the cases the model parameters φ (upper row) and θ (lower row) are minimized and corresponding to the first (left column) and the second (right column) best-fit points in Table 8

Summary and Conclusions
In this work, we have investigated the capability of the proposed long-baseline neutrino oscillation experiment ESSnuSB to discriminate between a class of lepton flavor models along with testing the viability of such models. In the framework of the discrete symmetry approach to lepton flavor, we have reviewed various lepton mixing patterns arising from breaking a flavor symmetry G f (or its CP version G CP ) to residual symmetries G e and G ν of the charged lepton and neutrino mass matrices. We have classified these patterns according to the number of free True Model-2.5 (δ CP < 0)  Figure 8: Comparison among the five allowed two-parameter models as test models for ESSnuSB using priors in the cases the model parameters φ (upper row) and θ (lower row) are minimized and corresponding to the first (left column) and the second (right column) best-fit points in Table 8 for Model 2.5 as the true model. For description of each panel, see Fig. 3.
parameters entering the predicted form of the leptonic mixing matrix U PMNS . Further, we have concentrated on one-and two-parameter patterns originating from relatively small (in terms of the number of elements) discrete groups. Namely, we have considered eleven one-parameter models with the free parameter θ, arising from G CP = S 4 CP [24] and A 5 CP [31] (see also Refs. [32,33]) broken to G e > Z 2 and G ν = Z 2 × CP, and seven two-parameter models with the free parameters θ and φ, originating from G f = A 4 , S 4 , and A 5 broken to G e (G ν ) = Z 2 and G ν (G e ) > Z 2 [20,21]. Since both residual symmetries are non-trivial, but one of them is Z 2 , these models have a relatively high predictive power, providing at the same time necessary freedom in fitting them to the data. First, to test the compatibility of the models with the current neutrino oscillation data, we have fitted the models with a simple χ 2 function using the current best-fit values of the three leptonic mixing angles as the true values, see Eqs. (3.1) and (3.2). We have listed all eleven one-and seven two-parameter models according to how well they fit the current data. We have also calculated the 3σ ranges of the free parameter θ for the one-parameter models as well as θ and φ for the two-parameter models, see Table 8. Second, using the 3σ ranges of the free parameters, we have calculated the predicted values of the leptonic mixing parameters for these models. We have found that among the eleven one-parameter models, six models are already excluded by the current data at more than 3σ, see Table 9 and as partly shown in Fig. 1. However, all seven two-parameter models are consistent with the individual 3σ ranges  Figure 9: Comparison between the best one-parameter model (Model 1.1) and the best twoparameter model (Model 2.1) for ESSnuSB using priors. The panel shows the quantity χ 2 , given Model 1.1 as the true model, as a function of the model parameter θ in its 3σ interval θ 3σ presented in Table 8. The value of χ 2 min for Model 1.1 is displayed with a black solid curve (when acting as a test model), whereas it is marked for Model 2.1 (i.e., the test model) by a black dotted curve.
of the leptonic mixing angles, see Fig. 2. For our analysis with ESSnuSB, we have selected only those models that are allowed by the current data within 3σ, which are Models 1.1-1.5 and 2.1-2.5. Among the five allowed one-parameter models, Models 1.3 and 1.5 predict | sin δ CP | = 1 and θ 23 = 45 • , whereas Models 1.1, 1.2, and 1.4 predict sin δ CP = 0 with θ 23 lying in the higher octant. For the two-parameter models, all models (i.e., Models 2.1-2.5) predict the higher octant of θ 23 . Furthermore, the predictions for δ CP of Models 2.3-2.5 are compatible with CP conservation, while Model 2.2 predicts δ CP close to maximal CP violation. On the other hand, Model 2.1 is compatible with neither CP conservation nor maximal CP violation.
Next, we have studied the capability of ESSnuSB to distinguish the different one-and twoparameter models. We have found that the one-parameter models, which predict sin δ CP = 0, i.e., Models 1.1, 1.2, and 1.4, can be separated from the one-parameter models that predict | sin δ CP | = 1, i.e., Models 1.3 and 1.5, by at least 7σ, as can be seen from Fig. 3. However, it is not possible to discriminate among Models 1.1, 1.2, and 1.4 and also not between Models 1.3 and 1.5. For the two-parameter models, we have found that ESSnuSB suggests Models 2.1-2.3 to be the most preferred models if they are assumed to be the true models for most of the parameter space, see Figs. 4-8. In general, ESSnuSB is helpful to disfavor the other models. In addition, we have studied the capability of ESSnuSB to discriminate between Models 1.1 and 2.1, which are the models that lead to the best fit to the current data for the one-and two-parameter models, respectively. We have found that Models 1.1 and 2.1 can be separated by 6σ or more (cf. Figs. 9 and 10), depending on the true parameter values assumed. Finally, while studying the capability of ESSnuSB to exclude models, we have found that if the best-fit point (after the running of ESSnuSB) remains close to the current best-fit point from global neutrino oscillation data, then Models 1.1, 1.2, and 1.4 are excluded by ESSnuSB at more than 5σ. However, Models 1.3, 1.5, 2.1, and 2.2 are consistent with the data at 5σ and Models 2.3-2.5 are consistent with the data at 3σ, see Figs. 11 and 12.  Figure 10: Comparison between the best two-parameter model (Model 2.1) and the best one-parameter model (Model 1.1) for ESSnuSB using priors. The upper (lower) panels show the quantity χ 2 min , given Model 2.1 as the true model, as function of the model parameter θ (φ) in its 3σ interval θ 3σ (φ 3σ ) presented in Table 8. In the panels, χ 2 min for Model 2.1 is displayed with black solid curves (when acting as a test model), whereas it is marked for Model 1.1 (i.e., the test model) by black dotted curves. Test Model-1.5 ≤ 1σ 1σ -3σ 3σ -5σ Figure 11: Compatibility of the five allowed one-parameter models with any potentially true values of sin 2 θ 23 and δ CP for ESSnuSB. For the other neutrino oscillation parameters, we assume the following true values sin 2 θ 12 = 0.310, sin 2 θ 13 = 0.02237, ∆m 2 21 = 7.39 · 10 −5 eV 2 , and ∆m 2 31 = 2.528 · 10 −3 eV 2 . A star (" ") indicates the present best-fit values of sin 2 θ 23 and δ CP from global neutrino oscillation data. The contours correspond to the indicated number of sigmas for 3 d.o.f. Test Model-2.5 ≤ 1σ 1σ -3σ 3σ -5σ Figure 12: Compatibility of the five allowed two-parameter models with any potentially true values of sin 2 θ 23 and δ CP for ESSnuSB. For the other neutrino oscillation parameters, we assume the following true values sin 2 θ 12 = 0.310, sin 2 θ 13 = 0.02237, ∆m 2 21 = 7.39 · 10 −5 eV 2 , and ∆m 2 31 = 2.528 · 10 −3 eV 2 . A star (" ") indicates the present best-fit values of sin 2 θ 23 and δ CP from global neutrino oscillation data. The contours correspond to the indicated number of sigmas for 2 d.o.f.